Recherche d'un algorithme efficace de racine carrée entière pour ARM Thumb2

Integer Racines Carrées par Jack W. Crenshaw pourrait être utile comme une autre référence.

Le C Extraits de l'Archive a également un integer racine carrée de la mise en œuvre . Celui-ci va au-delà juste le résultat entier, et calcule extra fractionnaire (point fixe) bits de la réponse. (Mise à jour: malheureusement, l'archive c snippets n'existe plus. Le lien renvoie à 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 nontrivial trial 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;                   /* trial 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. Il s'agit essentiellement de L'article de Wikipedia sur les méthodes de calcul à racine carrée . Mais il a été modifié pour utiliser stdint.h types uint32_t etc. À proprement parler, le type de retour peut ê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 le "arrondi" réponse. J'ai trouvé ce utile dans une certaine demande pour plus de précision. Notez que dans ce cas, le type de retour doit être uint32_t parce que la racine carrée arrondie de 2 ³² - 1 is 2 ¹⁶ .

/**
 * \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 requise, j'ai une approximation rapide pour vous, qui utilise 260octets de ram (vous pouvez réduire de moitié cela, mais ne le faites 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 la plage 1 - >2^20, l'erreur maximale est de 11, et sur la plage 1->2^30, elle est d'environ 256. Vous pourriez utiliser des tables plus grandes et minimiser cela. Il est intéressant de mentionner que l'erreur sera toujours négative-c.-à-d. quand il est faux, le la valeur sera inférieure à la valeur correcte.

vous pourriez faire 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 Ensuite, nous avons une table de recherche pour sqrt(2^N), et une table de recherche pour sqrt(1<=B<2). Nous stockons la table de recherche pour sqrt (2^n) Comme entier, ce qui pourrait être une erreur (les tests ne montrent aucun mauvais effets), et nous stockons le table de recherche pour sqrt (1<=B<2) à 15bits de point fixe.

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

en termes d'implémentation, le while (t) {cnt++;t>>=1;} est en fait une instruction count-leading-bits (CLB), donc si votre version du chipset a cela, vous gagnez! Aussi, le passage d'instruction serait facile de implémenter avec un shifter bidirectionnel, si vous en avez un? Il y a un algorithme Lg[N] pour compter le bit de réglage le plus élevé ici.

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

une approche commune est la bisection.

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 divise. Depuis la division est une division par 2, vous pouvez le remplacer par un >> .

la condition finale peut ne pas être spot on, alors assurez-vous de tester une variété d'entiers pour être sûr que la division par 2 ne oscille pas incorrectement entre deux valeurs égales; ils différerait 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 mis en œuvre d'une manière plus compliquée que nécessaire. J'ai repris l'idée d'ici: http://ww1.microchip.com/downloads/en/AppNotes/91040a.pdf (par Ross M. Fosler) et en fit une très courte fonction C:

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;
}

cela se compile à 5 cycles/bit sur mon blackfin. Je crois que votre code compilé sera en général plus rapide si vous utilisez des boucles au lieu de tout les boucles, et vous obtenez l'avantage supplémentaire du temps déterministe (bien que cela dépend dans une certaine mesure de la façon dont votre compilateur optimise l'instruction if).)

Ça dépend de l'utilisation de la fonction sqrt. J'en utilise souvent pour faire des versions rapides. Par exemple, quand 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 réglé à quelque chose de similaire à l'algorithme binaire chiffre par chiffre décrit dans cet article de Wikipedia .

Voici une solution en Java qui combine le log_2 entier et la méthode de Newton pour créer un algorithme sans boucle. En revanche, il a besoin de division. Les lignes commentées sont nécessaires pour passer à 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<<DeBruijnArray[x])*Math.sqrt(2)));
}

public static int sqrt(final int num)
{
  int y;
  if(num==0)
    return num;
  {
    int v= num;
    v|= v>>>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 de toit qui peuvent être produites.

cette méthode est similaire à la division longue: vous construisez une conjecture pour le prochain chiffre de la racine, faites une soustraction, et entrez le chiffre si la différence répond à certains critères. Avec la version binaire, votre seul choix pour le prochain chiffre 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 mis en œuvre la suggestion de Warren et la méthode de Newton en C# pour les entiers 64 bits. L'Isqrt utilise la méthode de Newton, tandis que L'Isqrt utilise la méthode de Warren. Voici le code source:

using System;

namespace Cluster
{
    public static class IntegerMath
    {


        /// <summary>
        /// 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.
        /// </summary>
        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 /q/looking-for-an-efficient-integer-square-root-algorithm-for-arm-thumb2-55830/"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 Release, Visual Studio 2012 ont été que la bibliothèque appelle Les Maths.La Sqrt est plus rapide.

• 59 ms-Newton (Isqrt)
• 12 ms - décalage de Bits (Isqrt2)
• 5 ms - Mathématique.Sqrt

Je ne suis pas intelligent avec les directives de compilateur, il peut donc être possible d'accorder le compilateur pour obtenir le nombre entier math plus rapidement. Clairement, le décalage de bits approche est très proche de la bibliothèque. Sur un système sans coprocesseur mathématique, il serait très rapide.

j'ai récemment rencontré la même tâche sur ARM Cortex-M3 (STM32F103CBT6) et après une recherche sur Internet est venu avec 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 de 1, C'est-à-dire LSB sur toute la gamme d'entrée UI32) et une vitesse relativement bonne (environ 1,3 m de racine carrée par seconde sur un cortex-M3 de 72 MHz ou environ 55 cycles par racine simple incluant 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 compiler 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