Racine de cube entier

Je recherche un code rapide pour les racines de cube 64 bits (non signé). (J’utilise C et comstack avec gcc, mais j’imagine que la plus grande partie du travail requirejs sera agnostique vis-à-vis du langage et du compilateur).

Étant donné une entrée n, je demande que la valeur de retour (intégrale) r soit telle que

r * r * r <= n && n < (r + 1) * (r + 1) * (r + 1) 

C’est-à-dire que je veux la racine cubique de n, arrondie. Code de base comme

 return (ulong)pow(n, 1.0/3); 

est incorrect en raison de l’arrondi vers la fin de la plage. Code non sophistiqué comme

 ulong cuberoot(ulong n) { ulong ret = pow(n + 0.5, 1.0/3); if (n = 18446724184312856125ULL) return 2642245ULL; if (ret * ret * ret > n) { ret--; while (ret * ret * ret > n) ret--; return ret; } while ((ret + 1) * (ret + 1) * (ret + 1) <= n) ret++; return ret; } 

donne le résultat correct, mais est plus lent que nécessaire.

Ce code est destiné à une bibliothèque mathématique et sera appelé plusieurs fois à partir de différentes fonctions. La rapidité est importante, mais vous ne pouvez pas compter sur une mémoire cache chaude (les suggestions telles qu’une recherche binary de 2 642 245 entrées sont donc parfaitement justes).

À des fins de comparaison, voici le code qui calcule correctement la racine carrée entière.

 ulong squareroot(ulong a) { ulong x = (ulong)sqrt((double)a); if (x > 0xFFFFFFFF || x*x > a) x--; return x; } 

    Le livre “Hacker’s Delight” contient des algorithmes pour résoudre ce problème et de nombreux autres. Le code est en ligne ici . EDIT : Ce code ne fonctionne pas correctement avec les ints 64 bits, et les instructions du manuel expliquant comment le réparer pour 64 bits sont quelque peu déroutantes. Une implémentation 64 bits appropriée (incluant le scénario de test) est en ligne ici .

    Je doute que votre fonction squareroot fonctionne “correctement” – cela devrait être ulong a pour l’argument, pas n 🙂 (mais la même approche fonctionnerait avec cbrt au lieu de sqrt , bien que toutes les bibliothèques de mathématiques C n’aient pas de fonctions de racine de cube).

    Vous pouvez essayer une solution de Newton pour corriger vos erreurs d’arrondi:

     ulong r = (ulong)pow(n, 1.0/3); if(r==0) return r; /* avoid divide by 0 later on */ ulong r3 = r*r*r; ulong slope = 3*r*r; ulong r1 = r+1; ulong r13 = r1*r1*r1; /* making sure to handle unsigned arithmetic correctly */ if(n >= r13) r+= (n - r3)/slope; if(n < r3) r-= (r3 - n)/slope; 

    Un seul pas de Newton devrait suffire, mais vous pouvez avoir des erreurs inédites (ou peut-être plus?). Vous pouvez vérifier / corriger ceux-ci en utilisant une étape finale de vérification et d’incrémentation, comme dans votre QO:

     while(r*r*r > n) --r; while((r+1)*(r+1)*(r+1) <= n) ++r; 

    ou quelque chose comme ça.

    (J'avoue que je suis paresseux; la bonne façon de le faire est de vérifier avec soin pour déterminer lequel des éléments de contrôle et d'incrément est réellement nécessaire ...)

    Si pow est trop cher, vous pouvez utiliser une instruction nombre-zéros-tête pour obtenir une approximation du résultat, puis utiliser une table de recherche, puis quelques étapes de Newton pour le terminer.

     int k = __builtin_clz(n); // counts # of leading zeros (often a single assembly insn) int b = 64 - k; // # of bits in n int top8 = n >> (b - 8); // top 8 bits of n (top bit is always 1) int approx = table[b][top8 & 0x7f]; 

    Étant donné b et top8 , vous pouvez utiliser une table de recherche (dans mon code, 8 000 entrées) pour trouver une bonne approximation de cuberoot(n) . Utilisez quelques étapes de Newton (voir la réponse de comingstorm) pour la terminer.

    Je voudrais rechercher comment le faire à la main , puis traduire cela en un algorithme informatique, travaillant en base 2 plutôt qu’en base 10.

    On se retrouve avec un algorithme du genre (pseudocode):

     Find the largest n such that (1 << 3n) < input. result = 1 << n. For i in (n-1)..0: if ((result | 1 << i)**3) < input: result |= 1 << i. 

    Nous pouvons optimiser le calcul de (result | 1 << i)**3 en observant que le bit à bit ou est équivalent à l'addition, le refactoring au result**3 + 3 * i * result ** 2 + 3 * i ** 2 * result + i ** 3 , en mettant en cache les valeurs du result**3 et du result**2 entre les itérations, et en utilisant des décalages au lieu de la multiplication.

     // On my pc: Math.Sqrt 35 ns, cbrt64 <70ns, cbrt32 <25 ns, (cbrt12 < 10ns) // cbrt64(ulong x) is a C# version of: // http://www.hackersdelight.org/hdcodetxt/acbrt.c.txt (acbrt1) // cbrt32(uint x) is a C# version of: // http://www.hackersdelight.org/hdcodetxt/icbrt.c.txt (icbrt1) // Union in C#: // http://www.hanselman.com/blog/UnionsOrAnEquivalentInCSairamasTipOfTheDay.aspx using System.Runtime.InteropServices; [StructLayout(LayoutKind.Explicit)] public struct fu_32 // float <==> uint { [FieldOffset(0)] public float f; [FieldOffset(0)] public uint u; } private static uint cbrt64(ulong x) { if (x >= 18446724184312856125) return 2642245; float fx = (float)x; fu_32 fu32 = new fu_32(); fu32.f = fx; uint uy = fu32.u / 4; uy += uy / 4; uy += uy / 16; uy += uy / 256; uy += 0x2a5137a0; fu32.u = uy; float fy = fu32.f; fy = 0.33333333f * (fx / (fy * fy) + 2.0f * fy); int y0 = (int) (0.33333333f * (fx / (fy * fy) + 2.0f * fy)); uint y1 = (uint)y0; ulong y2, y3; if (y1 >= 2642245) { y1 = 2642245; y2 = 6981458640025; y3 = 18446724184312856125; } else { y2 = (ulong)y1 * y1; y3 = y2 * y1; } if (y3 > x) { y1 -= 1; y2 -= 2 * y1 + 1; y3 -= 3 * y2 + 3 * y1 + 1; while (y3 > x) { y1 -= 1; y2 -= 2 * y1 + 1; y3 -= 3 * y2 + 3 * y1 + 1; } return y1; } do { y3 += 3 * y2 + 3 * y1 + 1; y2 += 2 * y1 + 1; y1 += 1; } while (y3 <= x); return y1 - 1; } private static uint cbrt32(uint x) { uint y = 0, z = 0, b = 0; int s = x < 1u << 24 ? x < 1u << 12 ? x < 1u << 06 ? x < 1u << 03 ? 00 : 03 : x < 1u << 09 ? 06 : 09 : x < 1u << 18 ? x < 1u << 15 ? 12 : 15 : x < 1u << 21 ? 18 : 21 : x >= 1u << 30 ? 30 : x < 1u << 27 ? 24 : 27; do { y *= 2; z *= 4; b = 3 * y + 3 * z + 1 << s; if (x >= b) { x -= b; z += 2 * y + 1; y += 1; } s -= 3; } while (s >= 0); return y; } private static uint cbrt12(uint x) // x < ~255 { uint y = 0, a = 0, b = 1, c = 0; while (a < x) { y++; b += c; a += b; c += 6; } if (a != x) y--; return y; }