Le moyen le plus rapide d’obtenir la partie entière de sqrt (n)?

Comme on le sait si n n’est pas un carré parfait, alors sqrt(n) ne serait pas un entier. Comme je n’ai besoin que de la partie entière, je pense que l’appel de sqrt(n) ne serait pas aussi rapide, car cela prend du temps de calculer la partie fractionnaire.

Donc, ma question est,

Peut-on obtenir uniquement la partie entière de sqrt (n) sans calculer la valeur réelle de sqrt(n) ? L’algorithme devrait être plus rapide que sqrt(n) (défini dans ou )?

Si possible, vous pouvez aussi écrire le code dans le bloc asm .

Je voudrais essayer le truc de racine carrée inverse rapide .

C’est une façon d’obtenir une très bonne approximation de 1/sqrt(n) sans n’importe quelle twig, basée sur des bit-twiddings donc pas portables (notamment entre les plates-formes 32 bits et 64 bits).

Une fois que vous l’avez obtenu, il vous suffit d’inverser le résultat et de prendre la partie entière.

Il peut y avoir des astuces plus rapides, bien sûr, car celui-ci est un peu compliqué.

EDIT : faisons-le!

D’abord un petit aide:

 // benchmark.h #include  template  double benchmark(Func f, size_t iterations) { f(); timeval a, b; gettimeofday(&a, 0); for (; iterations --> 0;) { f(); } gettimeofday(&b, 0); return (b.tv_sec * (unsigned int)1e6 + b.tv_usec) - (a.tv_sec * (unsigned int)1e6 + a.tv_usec); } 

Ensuite, le corps principal:

 #include  #include  #include "benchmark.h" class Sqrt { public: Sqrt(int n): _number(n) {} int operator()() const { double d = _number; return static_cast(std::sqrt(d) + 0.5); } private: int _number; }; // http://www.codecodex.com/wiki/Calculate_an_integer_square_root class IntSqrt { public: IntSqrt(int n): _number(n) {} int operator()() const { int remainder = _number; if (remainder < 0) { return 0; } int place = 1 <<(sizeof(int)*8 -2); while (place > remainder) { place /= 4; } int root = 0; while (place) { if (remainder >= root + place) { remainder -= root + place; root += place*2; } root /= 2; place /= 4; } return root; } private: int _number; }; // http://en.wikipedia.org/wiki/Fast_inverse_square_root class FastSqrt { public: FastSqrt(int n): _number(n) {} int operator()() const { float number = _number; float x2 = number * 0.5F; float y = number; long i = *(long*)&y; //i = (long)0x5fe6ec85e7de30da - (i >> 1); i = 0x5f3759df - (i >> 1); y = *(float*)&i; y = y * (1.5F - (x2*y*y)); y = y * (1.5F - (x2*y*y)); // let's be precise return static_cast(1/y + 0.5f); } private: int _number; }; int main(int argc, char* argv[]) { if (argc != 3) { std::cerr < < "Usage: %prog integer iterations\n"; return 1; } int n = atoi(argv[1]); int it = atoi(argv[2]); assert(Sqrt(n)() == IntSqrt(n)() && Sqrt(n)() == FastSqrt(n)() && "Different Roots!"); std::cout << "sqrt(" << n << ") = " << Sqrt(n)() << "\n"; double time = benchmark(Sqrt(n), it); double intTime = benchmark(IntSqrt(n), it); double fastTime = benchmark(FastSqrt(n), it); std::cout << "Number iterations: " << it << "\n" "Sqrt computation : " << time << "\n" "Int computation : " << intTime << "\n" "Fast computation : " << fastTime << "\n"; return 0; } 

Et les résultats:

 sqrt(82) = 9 Number iterations: 4096 Sqrt computation : 56 Int computation : 217 Fast computation : 119 // Note had to tweak the program here as Int here returns -1 :/ sqrt(2147483647) = 46341 // real answer sqrt(2 147 483 647) = 46 340.95 Number iterations: 4096 Sqrt computation : 57 Int computation : 313 Fast computation : 119 

Où, comme prévu, le calcul rapide fonctionne beaucoup mieux que le calcul Int .

Oh, et au fait, sqrt est plus rapide 🙂

Edit: cette réponse est stupide – utilisez (int) sqrt(i)

Après le profilage avec les parameters appropriés ( -march=native -m64 -O3 ), ce qui précède était beaucoup plus rapide.


Bon, une question un peu ancienne, mais la réponse “la plus rapide” n’a pas encore été donnée. Le plus rapide (je pense) est l’algorithme Binary Square Root, expliqué en détail dans cet article d’Embedded.com .

Cela revient essentiellement à ceci:

 unsigned short isqrt(unsigned long a) { unsigned long rem = 0; int root = 0; int i; for (i = 0; i < 16; i++) { root <<= 1; rem <<= 2; rem += a >> 30; a < <= 2; if (root < rem) { root++; rem -= root; root++; } } return (unsigned short) (root >> 1); } 

Sur ma machine (Q6600, Ubuntu 10.10), j’ai profilé en prenant la racine carrée des nombres 1-100000000. En utilisant iqsrt(i) fallu 2750 ms. L’utilisation de (unsigned short) sqrt((float) i) pris 3600ms. Cela a été fait en utilisant g++ -O3 . En utilisant l’option de compilation -ffast-math , les temps étaient 2100ms et 3100ms respectivement. Notez que ceci est sans utiliser une seule ligne d’assembleur, donc cela pourrait probablement être beaucoup plus rapide.

Le code ci-dessus fonctionne à la fois pour C et C ++ et avec des modifications de syntaxe mineures également pour Java.

Ce qui fonctionne encore mieux pour une plage limitée est une recherche binary. Sur ma machine, la version ci-dessus sort de l’eau d’un facteur 4. Malheureusement, sa scope est très limitée:

 #include  const uint16_t squares[] = { 0, 1, 4, 9, 16, 25, 36, 49, 64, 81, 100, 121, 144, 169, 196, 225, 256, 289, 324, 361, 400, 441, 484, 529, 576, 625, 676, 729, 784, 841, 900, 961, 1024, 1089, 1156, 1225, 1296, 1369, 1444, 1521, 1600, 1681, 1764, 1849, 1936, 2025, 2116, 2209, 2304, 2401, 2500, 2601, 2704, 2809, 2916, 3025, 3136, 3249, 3364, 3481, 3600, 3721, 3844, 3969, 4096, 4225, 4356, 4489, 4624, 4761, 4900, 5041, 5184, 5329, 5476, 5625, 5776, 5929, 6084, 6241, 6400, 6561, 6724, 6889, 7056, 7225, 7396, 7569, 7744, 7921, 8100, 8281, 8464, 8649, 8836, 9025, 9216, 9409, 9604, 9801, 10000, 10201, 10404, 10609, 10816, 11025, 11236, 11449, 11664, 11881, 12100, 12321, 12544, 12769, 12996, 13225, 13456, 13689, 13924, 14161, 14400, 14641, 14884, 15129, 15376, 15625, 15876, 16129, 16384, 16641, 16900, 17161, 17424, 17689, 17956, 18225, 18496, 18769, 19044, 19321, 19600, 19881, 20164, 20449, 20736, 21025, 21316, 21609, 21904, 22201, 22500, 22801, 23104, 23409, 23716, 24025, 24336, 24649, 24964, 25281, 25600, 25921, 26244, 26569, 26896, 27225, 27556, 27889, 28224, 28561, 28900, 29241, 29584, 29929, 30276, 30625, 30976, 31329, 31684, 32041, 32400, 32761, 33124, 33489, 33856, 34225, 34596, 34969, 35344, 35721, 36100, 36481, 36864, 37249, 37636, 38025, 38416, 38809, 39204, 39601, 40000, 40401, 40804, 41209, 41616, 42025, 42436, 42849, 43264, 43681, 44100, 44521, 44944, 45369, 45796, 46225, 46656, 47089, 47524, 47961, 48400, 48841, 49284, 49729, 50176, 50625, 51076, 51529, 51984, 52441, 52900, 53361, 53824, 54289, 54756, 55225, 55696, 56169, 56644, 57121, 57600, 58081, 58564, 59049, 59536, 60025, 60516, 61009, 61504, 62001, 62500, 63001, 63504, 64009, 64516, 65025 }; inline int isqrt(uint16_t x) { const uint16_t *p = squares; if (p[128] < = x) p += 128; if (p[ 64] <= x) p += 64; if (p[ 32] <= x) p += 32; if (p[ 16] <= x) p += 16; if (p[ 8] <= x) p += 8; if (p[ 4] <= x) p += 4; if (p[ 2] <= x) p += 2; if (p[ 1] <= x) p += 1; return p - squares; } 

Une version 32 bits peut être téléchargée ici: https://gist.github.com/3481770

Bien que je soupçonne que vous pouvez trouver beaucoup d’options en recherchant “la racine carrée de nombre entier rapide”, voici quelques idées potentiellement nouvelles qui pourraient bien fonctionner (chaque indépendant, ou peut-être vous pouvez les combiner):

  1. Créez un tableau static const de tous les carrés parfaits du domaine que vous souhaitez prendre en charge et effectuez une recherche binary rapide sans twig sur celui-ci. L’index résultant dans le tableau est la racine carrée.
  2. Convertissez le nombre en virgule flottante et divisez-le en mantisse et exposant. Réduisez de moitié l’exposant et multipliez la mantisse par un facteur magique (votre travail pour le trouver). Cela devrait pouvoir vous donner une approximation très proche. Incluez une dernière étape pour l’ajuster si elle n’est pas exacte (ou utilisez-la comme sharepoint départ pour la recherche binary ci-dessus).

Je pense que la Google search fournit de bons articles comme Calculate an integer square root qui traitait de trop de manières possibles de calculer rapidement et qu’il existe de bons articles de référence. it), mais si vous les lisez et qu’il y a une ambiguïté avec eux, il se peut que nous puissions vous aider.

Si cela ne vous dérange pas une approximation, que diriez-vous de cette fonction de sqrt entier que j’ai bricolé.

 int sqrti(int x) { union { float f; int x; } v; // convert to float vf = (float)x; // fast aprox sqrt // assumes float is in IEEE 754 single precision format // assumes int is 32 bits // b = exponent bias // m = number of mantissa bits vx -= 1 < < 23; // subtract 2^m vx >>= 1; // divide by 2 vx += 1 < < 29; // add ((b + 1) / 2) * 2^m // convert to int return (int)vf; } 

Il utilise l'algorithme décrit dans cet article Wikipedia . Sur ma machine, il est presque deux fois plus rapide que sqrt 🙂

Pour faire un nombre entier de sqrt, vous pouvez utiliser cette spécialisation de la méthode newtons:

 Def isqrt(N): a = 1 b = N while |ab| > 1 b = N / a a = (a + b) / 2 return a 

Pour tout x, le sqrt se situe dans la plage (x … N / x), nous divisons donc cet intervalle à chaque boucle pour la nouvelle estimation. Un peu comme la recherche binary mais elle converge plus vite.

Cela converge dans O (loglog (N)) qui est très rapide. Il n’utilise pas non plus de virgule flottante, et il fonctionnera aussi bien pour des entiers de précision arbitraires.

Pourquoi personne ne suggère la méthode la plus rapide?

Si:

  1. la gamme de nombres est limitée
  2. la consommation de mémoire n’est pas cruciale
  3. le temps de lancement de l’application n’est pas critique

puis créez int[MAX_X] rempli (au lancement) avec sqrt(x) (vous n’avez pas besoin d’utiliser la fonction sqrt() pour cela).

Toutes ces conditions correspondent bien à mon programme. En particulier, un tableau int[10000000] va consumr 40MB

Que pensez-vous de cela?

Dans de nombreux cas, même la valeur exacte de sqrt n’est pas nécessaire, en ayant une bonne approximation. (Par exemple, cela arrive souvent dans l’optimisation DSP, lorsque le signal 32 bits doit être compressé en 16 bits ou 16 bits en 8 bits, sans perdre beaucoup de précision autour de zéro).

J’ai trouvé cette équation utile:

 k = ceil(MSB(n)/2); - MSB(n) is the most significant bit of "n" 

 sqrt(n) ~= 2^(k-2)+(2^(k-1))*n/(2^(2*k))); - all multiplications and divisions here are very DSP-friendly, as they are only 2^k. 

Cette équation génère une courbe lisse (n, sqrt (n)), ses valeurs ne sont pas très différentes de sqrt réel (n) et peuvent donc être utiles lorsque la précision approximative est suffisante.

Si vous avez besoin de performances sur le calcul de la racine carrée, je suppose que vous en calculerez beaucoup. Alors pourquoi ne pas mettre la réponse en cache? Je ne connais pas la plage pour N dans votre cas, ni si vous calculerez plusieurs fois la racine carrée du même entier, mais si oui, alors vous pouvez mettre le résultat en cache à chaque fois que votre méthode est appelée (dans un tableau serait le plus efficace sinon trop grand).

C’est tellement court que cela correspond à 99%.

 static inline int sqrtn(int num) { int i; __asm__ ( "pxor %%xmm0, %%xmm0\n\t" // clean xmm0 for cvtsi2ss "cvtsi2ss %1, %%xmm0\n\t" // convert num to float, put it to xmm0 "sqrtss %%xmm0, %%xmm0\n\t" // square root xmm0 "cvttss2si %%xmm0, %0" // float to int :"=r"(i):"r"(num):"%xmm0"); // i: result, num: input, xmm0: scratch register return i; } 

Pourquoi nettoyer xmm0 ? Documentation de cvtsi2ss

L’opérande de destination est un registre XMM. Le résultat est stocké dans le double mot bas de l’opérande de destination et les trois mots doubles supérieurs restnt inchangés.

Version insortingnsèque de GCC (s’exécute uniquement sur GCC):

 #include  int sqrtn2(int num) { register __v4sf xmm0 = {0, 0, 0, 0}; xmm0 = __builtin_ia32_cvtsi2ss(xmm0, num); xmm0 = __builtin_ia32_sqrtss(xmm0); return __builtin_ia32_cvttss2si(xmm0); } 

Version Intel Insortingnsic (testée sur GCC, Clang, ICC):

 #include  int sqrtn2(int num) { register __m128 xmm0 = _mm_setzero_ps(); xmm0 = _mm_cvt_si2ss(xmm0, num); xmm0 = _mm_sqrt_ss(xmm0); return _mm_cvt_ss2si(xmm0); } 

^^^^ Tous nécessitent SSE 1. (même pas SSE 2)

Sur mon ordinateur avec gcc, avec -ffast-math, convertir un entier 32 bits en float et en utilisant sqrtf prend 1,2 s par 10 ^ 9 opérations (sans -ffast-math cela prend 3,54 s).

L’algorithme suivant utilise 0,87 s par 10 ^ 9 aux dépens d’une certaine précision: les erreurs peuvent aller jusqu’à -7 ou +1 bien que l’erreur RMS ne soit que de 0,79:

 uint16_t SQRTTAB[65536]; inline uint16_t approxsqrt(uint32_t x) { const uint32_t m1 = 0xff000000; const uint32_t m2 = 0x00ff0000; if (x&m1) { return SQRTTAB[x>>16]; } else if (x&m2) { return SQRTTAB[x>>8]>>4; } else { return SQRTTAB[x]>>8; } } 

La table est construite en utilisant:

 void maketable() { for (int x=0; x<65536; x++) { double v = x/65535.0; v = sqrt(v); int y = int(v*65535.0+0.999); SQRTTAB[x] = y; } } 

J'ai trouvé que le fait d'affiner la bissection en utilisant davantage d'instructions améliore la précision, mais cela ralentit aussi les choses au point que sqrtf est plus rapide, du moins avec -ffast-math.