Vous recherchez un algorithme de racine carrée en nombre entier efficace pour ARM Thumb2

Je cherche un algorithme rapide, entier seulement pour trouver la racine carrée (partie entière de celui-ci) d’un entier non signé. Le code doit avoir d’excellentes performances sur les processeurs ARM Thumb 2. Il pourrait s’agir d’un langage d’assemblage ou d’un code C.

Tous les conseils sont les bienvenus.

Integer Square Roots de Jack W. Crenshaw pourrait être une autre référence.

L’ archive C Snippets a également une implémentation de racine carrée entière . Celui-ci va au-delà du résultat entier et calcule des bits fractionnaires supplémentaires (en virgule fixe) de la réponse. (Mise à jour: malheureusement, l’archive des extraits de code C est désormais obsolète. Le lien pointe vers l’archive Web de la page.) Voici le code de l’archive C Snippets:

#define BITSPERLONG 32 #define TOP2BITS(x) ((x & (3L << (BITSPERLONG-2))) >> (BITSPERLONG-2)) struct int_sqrt { unsigned sqrt, frac; }; /* usqrt: ENTRY x: unsigned long EXIT returns floor(sqrt(x) * pow(2, BITSPERLONG/2)) Since the square root never uses more than half the bits of the input, we use the other half of the bits to contain extra bits of precision after the binary point. EXAMPLE suppose BITSPERLONG = 32 then usqrt(144) = 786432 = 12 * 65536 usqrt(32) = 370727 = 5.66 * 65536 NOTES (1) change BITSPERLONG to BITSPERLONG/2 if you do not want the answer scaled. Indeed, if you want n bits of precision after the binary point, use BITSPERLONG/2+n. The code assumes that BITSPERLONG is even. (2) This is really better off being written in assembly. The line marked below is really a "arithmetic shift left" on the double-long value with r in the upper half and x in the lower half. This operation is typically expressible in only one or two assembly instructions. (3) Unrolling this loop is probably not a bad idea. ALGORITHM The calculations are the base-two analogue of the square root algorithm we all learned in grammar school. Since we're in base 2, there is only one nonsortingvial sortingal multiplier. Notice that absolutely no multiplications or divisions are performed. This means it'll be fast on a wide range of processors. */ void usqrt(unsigned long x, struct int_sqrt *q) { unsigned long a = 0L; /* accumulator */ unsigned long r = 0L; /* remainder */ unsigned long e = 0L; /* sortingal product */ int i; for (i = 0; i < BITSPERLONG; i++) /* NOTE 1 */ { r = (r << 2) + TOP2BITS(x); x <<= 2; /* NOTE 2 */ a <<= 1; e = (a << 1) + 1; if (r >= e) { r -= e; a++; } } memcpy(q, &a, sizeof(long)); } 

J’ai choisi le code suivant. C’est essentiellement de l’article de Wikipedia sur les méthodes de calcul de la racine carrée . Mais il a été modifié pour utiliser les types stdint.h uint32_t etc. Ssortingctement parlant, le type de retour pourrait être changé en uint16_t .

 /** * \brief Fast Square root algorithm * * Fractional parts of the answer are discarded. That is: * - SquareRoot(3) --> 1 * - SquareRoot(4) --> 2 * - SquareRoot(5) --> 2 * - SquareRoot(8) --> 2 * - SquareRoot(9) --> 3 * * \param[in] a_nInput - unsigned integer for which to find the square root * * \return Integer square root of the input value. */ uint32_t SquareRoot(uint32_t a_nInput) { uint32_t op = a_nInput; uint32_t res = 0; uint32_t one = 1uL << 30; // The second-to-top bit is set: use 1u << 14 for uint16_t type; use 1uL<<30 for uint32_t type // "one" starts at the highest power of four <= than the argument. while (one > op) { one >>= 2; } while (one != 0) { if (op >= res + one) { op = op - (res + one); res = res + 2 * one; } res >>= 1; one >>= 2; } return res; } 

La bonne chose, j’ai découvert, est qu’une modification assez simple peut retourner la réponse “arrondie”. J’ai trouvé cela utile dans une application donnée pour plus de précision. Notez que dans ce cas, le type de retour doit être uint32_t car la racine carrée arrondie de 2 32 – 1 est 2 16 .

 /** * \brief Fast Square root algorithm, with rounding * * This does arithmetic rounding of the result. That is, if the real answer * would have a fractional part of 0.5 or greater, the result is rounded up to * the next integer. * - SquareRootRounded(2) --> 1 * - SquareRootRounded(3) --> 2 * - SquareRootRounded(4) --> 2 * - SquareRootRounded(6) --> 2 * - SquareRootRounded(7) --> 3 * - SquareRootRounded(8) --> 3 * - SquareRootRounded(9) --> 3 * * \param[in] a_nInput - unsigned integer for which to find the square root * * \return Integer square root of the input value. */ uint32_t SquareRootRounded(uint32_t a_nInput) { uint32_t op = a_nInput; uint32_t res = 0; uint32_t one = 1uL << 30; // The second-to-top bit is set: use 1u << 14 for uint16_t type; use 1uL<<30 for uint32_t type // "one" starts at the highest power of four <= than the argument. while (one > op) { one >>= 2; } while (one != 0) { if (op >= res + one) { op = op - (res + one); res = res + 2 * one; } res >>= 1; one >>= 2; } /* Do arithmetic rounding to nearest integer */ if (op > res) { res++; } return res; } 

Si la précision exacte n’est pas requirejse, j’ai une approximation rapide pour vous, qui utilise 260 octets de ram (vous pouvez réduire de moitié cela, mais pas).

 int ftbl[33]={0,1,1,2,2,4,5,8,11,16,22,32,45,64,90,128,181,256,362,512,724,1024,1448,2048,2896,4096,5792,8192,11585,16384,23170,32768,46340}; int ftbl2[32]={ 32768,33276,33776,34269,34755,35235,35708,36174,36635,37090,37540,37984,38423,38858,39287,39712,40132,40548,40960,41367,41771,42170,42566,42959,43347,43733,44115,44493,44869,45241,45611,45977}; int fisqrt(int val) { int cnt=0; int t=val; while (t) {cnt++;t>>=1;} if (6>=cnt) t=(val<<(6-cnt)); else t=(val>>(cnt-6)); return (ftbl[cnt]*ftbl2[t&31])>>15; } 

Voici le code pour générer les tables:

 ftbl[0]=0; for (int i=0;i<32;i++) ftbl[i+1]=sqrt(pow(2.0,i)); printf("int ftbl[33]={0"); for (int i=0;i<32;i++) printf(",%d",ftbl[i+1]); printf("};\n"); for (int i=0;i<32;i++) ftbl2[i]=sqrt(1.0+i/32.0)*32768; printf("int ftbl2[32]={"); for (int i=0;i<32;i++) printf("%c%d",(i)?',':' ',ftbl2[i]); printf("};\n"); 

Sur l'intervalle 1-> 2 ^ 20, l'erreur maximale est de 11 et sur l'intervalle 1-> 2 ^ 30, elle est d'environ 256. Vous pouvez utiliser des tableaux plus grands et minimiser cela. Il convient de mentionner que l'erreur sera toujours négative, c'est-à-dire que lorsque la valeur est incorrecte, la valeur sera INFERIEURE à la valeur correcte.

Vous feriez bien de suivre cela avec une étape de raffinage.

L'idée est assez simple: (ab) ^ 0.5 = a ^ 0.b * b ^ 0.5.

Donc, nous prenons l'entrée X = A * B où A = 2 ^ N et 1 <= B <2 Nous avons alors une recherche pour sqrt (2 ^ N), et une recherche pour sqrt (1 <= B <2) . Nous stockons le lookuptable pour sqrt (2 ^ N) en entier, ce qui peut être une erreur (les tests ne montrent aucun effet négatif), et nous stockons la recherche de sqrt (1 <= B <2) à 15 bits de point fixe.

Nous soaps que 1 <= sqrt (2 ^ N) <65536, donc 16bit, et nous savons que nous ne pouvons vraiment multiplier que 16bitx15bit sur un ARM, sans crainte de représailles, alors c'est ce que nous faisons.

En termes d'implémentation, le while (t) {cnt ++; t >> = 1;} est en fait une instruction de décompte de bits (CLB), donc si votre version du chipset a cela, vous gagnez! En outre, l’instruction de décalage serait facile à mettre en œuvre avec un levier de vitesses bidirectionnel, si vous en avez un? Il y a un algorithme Lg [N] pour compter le bit le plus élevé ici.

En termes de nombres magiques, pour changer la taille des tables, le nombre magique de ftbl2 est de 32, mais notez que 6 (Lg [32] +1) est utilisé pour le décalage.

Une approche commune est la bissection.

 hi = number lo = 0 mid = ( hi + lo ) / 2 mid2 = mid*mid while( lo < hi-1 and mid2 != number ) { if( mid2 < number ) { lo = mid else hi = mid mid = ( hi + lo ) / 2 mid2 = mid*mid 

Quelque chose comme ça devrait fonctionner raisonnablement bien. Il fait des tests log2 (nombre), en faisant log2 (nombre) se multiplie et se divise. Comme la division est une division par 2, vous pouvez la remplacer par un >> .

La condition de terminaison peut ne pas être exacte, alors assurez-vous de tester divers entiers pour vous assurer que la division par 2 n'oscille pas de manière incorrecte entre deux valeurs paires. ils différeraient de plus de 1.

Ce n’est pas rapide mais c’est petit et simple:

 int isqrt(int n) { int b = 0; while(n >= 0) { n = n - b; b = b + 1; n = n - b; } return b - 1; } 

Je trouve que la plupart des algorithmes sont basés sur des idées simples, mais sont implémentés de manière plus compliquée que nécessaire. J’ai pris l’idée ici: http://ww1.microchip.com/downloads/en/AppNotes/91040a.pdf (par Ross M. Fosler) et en ai fait une fonction C très courte:

 uint16_t int_sqrt32(uint32_t x) { uint16_t res=0; uint16_t add= 0x8000; int i; for(i=0;i<16;i++) { uint16_t temp=res | add; uint32_t g2=temp*temp; if (x>=g2) { res=temp; } add>>=1; } return res; } 

Ceci comstack à 5 cycles / bit sur mon blackfin. Je pense que votre code compilé sera en général plus rapide si vous utilisez des boucles plutôt que des boucles while et vous bénéficiez de l’avantage du temps déterministe (bien que cela dépend dans une certaine mesure de la façon dont votre compilateur optimise l’instruction if).

Cela dépend de l’utilisation de la fonction sqrt. J’utilise souvent des approx pour faire des versions rapides. Par exemple, lorsque je dois calculer le module de vecteur:

 Module = SQRT( x^2 + y^2) 

J’utilise :

 Module = MAX( x,y) + Min(x,y)/2 

Qui peut être codé en 3 ou 4 instructions comme:

 If (x > y ) Module = x + y >> 1; Else Module = y + x >> 1; 

J’ai opté pour quelque chose de similaire à l’algorithme binary digit-by-digit décrit dans cet article de Wikipedia .

Voici une solution en Java qui combine un entier log_2 et la méthode de Newton pour créer un algorithme sans boucle. En contrepartie, il faut une division. Les lignes commentées doivent être converties en un algorithme 64 bits.

 private static final int debruijn= 0x07C4ACDD; //private static final long debruijn= ( ~0x0218A392CD3D5DBFL)>>>6; static { for(int x= 0; x<32; ++x) { final long v= ~( -2L<<(x)); DeBruijnArray[(int)((v*debruijn)>>>27)]= x; //>>>58 } for(int x= 0; x<32; ++x) SQRT[x]= (int) (Math.sqrt((1L<>>1; // first round up to one less than a power of 2 v|= v>>>2; v|= v>>>4; v|= v>>>8; v|= v>>>16; //v|= v>>>32; y= SQRT[(v*debruijn)>>>27]; //>>>58 } //y= (y+num/y)>>>1; y= (y+num/y)>>>1; y= (y+num/y)>>>1; y= (y+num/y)>>>1; return y*y>num?y-1:y; } 

Comment cela fonctionne: La première partie produit une racine carrée précise à environ trois bits. La ligne [y = (y + num / y) >> 1;] double la précision en bits. La dernière ligne élimine les racines du toit pouvant être générées.

Cette méthode est similaire à la division longue: vous construisez une hypothèse pour le chiffre suivant de la racine, effectuez une soustraction et entrez le chiffre si la différence répond à certains critères. Avec une version binary, votre seul choix pour le chiffre suivant est 0 ou 1, donc vous devinez toujours 1, faites la soustraction et entrez un 1 sauf si la différence est négative.

http://www.realitypixels.com/turk/opensource/index.html#FractSqrt

J’ai implémenté la suggestion de Warren et la méthode de Newton en C # pour les entiers 64 bits. Isqrt utilise la méthode de Newton, tandis que Isqrt utilise la méthode de Warren. Voici le code source:

 using System; namespace Cluster { public static class IntegerMath { ///  /// Compute the integer square root, the largest whole number less than or equal to the true square root of N. /// /// This uses the integer version of Newton's method. ///  public static long Isqrt(this long n) { if (n < 0) throw new ArgumentOutOfRangeException("n", "Square root of negative numbers is not defined."); if (n <= 1) return n; var xPrev2 = -2L; var xPrev1 = -1L; var x = 2L; // From Wikipedia: if N + 1 is a perfect square, then the algorithm enters a two-value cycle, so we have to compare // to two previous values to test for convergence. while (x != xPrev2 && x != xPrev1) { xPrev2 = xPrev1; xPrev1 = x; x = (x + n/x)/2; } // The two values x and xPrev1 will be above and below the true square root. Choose the lower one. return x < xPrev1 ? x : xPrev1; } #region Sqrt using Bit-shifting and magic numbers. // From http://stackoverflow.com/questions/1100090/looking-for-an-efficient-integer-square-root-algorithm-for-arm-thumb2 // Converted to C#. private static readonly ulong debruijn= ( ~0x0218A392CD3D5DBFUL)>>6; private static readonly ulong[] SQRT = new ulong[64]; private static readonly int[] DeBruijnArray = new int[64]; static IntegerMath() { for(int x= 0; x<64; ++x) { ulong v= (ulong) ~( -2L<<(x)); DeBruijnArray[(v*debruijn)>>58]= x; } for(int x= 0; x<64; ++x) SQRT[x]= (ulong) (Math.Sqrt((1L<>1; // first round up to one less than a power of 2 v|= v>>2; v|= v>>4; v|= v>>8; v|= v>>16; v|= v>>32; y= SQRT[(v*debruijn)>>58]; } y= (y+num/y)>>1; y= (y+num/y)>>1; y= (y+num/y)>>1; y= (y+num/y)>>1; // Make sure that our answer is rounded down, not up. return (long) (y*y>num?y-1:y); } #endregion } } 

J’ai utilisé ce qui suit pour évaluer le code:

 using System; using System.Diagnostics; using Cluster; using Microsoft.VisualStudio.TestTools.UnitTesting; namespace ClusterTests { [TestClass] public class IntegerMathTests { [TestMethod] public void Isqrt_Accuracy() { for (var n = 0L; n <= 100000L; n++) { var expectedRoot = (long) Math.Sqrt(n); var actualRoot = n.Isqrt(); Assert.AreEqual(expectedRoot, actualRoot, String.Format("Square root is wrong for N = {0}.", n)); } } [TestMethod] public void Isqrt2_Accuracy() { for (var n = 0L; n <= 100000L; n++) { var expectedRoot = (long)Math.Sqrt(n); var actualRoot = n.Isqrt2(); Assert.AreEqual(expectedRoot, actualRoot, String.Format("Square root is wrong for N = {0}.", n)); } } [TestMethod] public void Isqrt_Speed() { var integerTimer = new Stopwatch(); var libraryTimer = new Stopwatch(); integerTimer.Start(); var total = 0L; for (var n = 0L; n <= 300000L; n++) { var root = n.Isqrt(); total += root; } integerTimer.Stop(); libraryTimer.Start(); total = 0L; for (var n = 0L; n <= 300000L; n++) { var root = (long)Math.Sqrt(n); total += root; } libraryTimer.Stop(); var isqrtMilliseconds = integerTimer.ElapsedMilliseconds; var libraryMilliseconds = libraryTimer.ElapsedMilliseconds; var msg = String.Format("Isqrt: {0} ms versus library: {1} ms", isqrtMilliseconds, libraryMilliseconds); Debug.WriteLine(msg); Assert.IsTrue(libraryMilliseconds > isqrtMilliseconds, "Isqrt2 should be faster than Math.Sqrt! " + msg); } [TestMethod] public void Isqrt2_Speed() { var integerTimer = new Stopwatch(); var libraryTimer = new Stopwatch(); var warmup = (10L).Isqrt2(); integerTimer.Start(); var total = 0L; for (var n = 0L; n <= 300000L; n++) { var root = n.Isqrt2(); total += root; } integerTimer.Stop(); libraryTimer.Start(); total = 0L; for (var n = 0L; n <= 300000L; n++) { var root = (long)Math.Sqrt(n); total += root; } libraryTimer.Stop(); var isqrtMilliseconds = integerTimer.ElapsedMilliseconds; var libraryMilliseconds = libraryTimer.ElapsedMilliseconds; var msg = String.Format("isqrt2: {0} ms versus library: {1} ms", isqrtMilliseconds, libraryMilliseconds); Debug.WriteLine(msg); Assert.IsTrue(libraryMilliseconds > isqrtMilliseconds, "Isqrt2 should be faster than Math.Sqrt! " + msg); } } } 

Mes résultats sur un Dell Latitude E6540 en mode Libération, Visual Studio 2012 étaient que l’appel de bibliothèque Math.Sqrt est plus rapide.

  • 59 ms – Newton (Isqrt)
  • 12 ms – Décalage de bits (Isqrt2)
  • 5 ms – Math.Sqrt

Je ne suis pas habile avec les directives du compilateur, il est donc possible d’ajuster le compilateur pour obtenir plus rapidement le calcul de l’entier. Clairement, l’approche de changement de bit est très proche de la bibliothèque. Sur un système sans coprocesseur mathématique, ce serait très rapide.

J’ai récemment rencontré la même tâche sur ARM Cortex-M3 (STM32F103CBT6) et après avoir effectué une recherche sur Internet, j’ai trouvé la solution suivante. Ce n’est pas la comparaison la plus rapide avec les solutions proposées ici, mais elle a une bonne précision (erreur maximale 1, c.-à-d. LSB sur toute la plage d’entrée UI32) et relativement bonne vitesse (environ 1,3 million de racines carrées par seconde). M3 ou environ 55 cycles par racine unique, y compris l’appel de fonction).

 // FastIntSqrt is based on Wikipedia article: // https://en.wikipedia.org/wiki/Methods_of_computing_square_roots // Which involves Newton's method which gives the following iterative formula: // // X(n+1) = (X(n) + S/X(n))/2 // // Thanks to ARM CLZ instruction (which counts how many bits in a number are // zeros starting from the most significant one) we can very successfully // choose the starting value, so just three iterations are enough to achieve // maximum possible error of 1. The algorithm uses division, but fortunately // it is fast enough here, so square root computation takes only about 50-55 // cycles with maximum comstackr optimization. uint32_t FastIntSqrt (uint32_t value) { if (!value) return 0; uint32_t xn = 1 << ((32 - __CLZ (value))/2); xn = (xn + value/xn)/2; xn = (xn + value/xn)/2; xn = (xn + value/xn)/2; return xn; } 

J'utilise IAR et il produit le code assembleur suivant:

  SECTION `.text`:CODE:NOROOT(1) THUMB _Z11FastIntSqrtj: MOVS R1,R0 BNE.N ??FastIntSqrt_0 MOVS R0,#+0 BX LR ??FastIntSqrt_0: CLZ R0,R1 RSB R0,R0,#+32 MOVS R2,#+1 LSRS R0,R0,#+1 LSL R0,R2,R0 UDIV R3,R1,R0 ADDS R0,R3,R0 LSRS R0,R0,#+1 UDIV R2,R1,R0 ADDS R0,R2,R0 LSRS R0,R0,#+1 UDIV R1,R1,R0 ADDS R0,R1,R0 LSRS R0,R0,#+1 BX LR ;; return