Vous cherchez un efficace entier racine carrée de l'algorithme pour les BRAS Thumb2

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

La C Extraits De L'Archive a aussi un entier racine carrée de la mise en œuvre. Celui-ci va au-delà de tout juste le résultat sous forme d'entier, et calcule supplémentaire fractionnaire (point fixe) bits de la réponse. (Mise à jour: malheureusement, le C des extraits d'archives est désormais obsolète. Le lien pointe vers un site web des archives de la page.) Voici le code à partir de la C des Extraits de l'Archive:

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

Que j'ai réglé sur le code suivant. C'est essentiellement à partir de la Article de Wikipedia sur la racine carrée méthodes de calcul. Mais il a été changé à utiliser stdint.h types uint32_t etc. Strictement parlant, le type de retour peut être changée pour 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 assez simple modification 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 car l'arrondi de la racine carrée de 2³² - 1 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;
}

Une approche commune est non-bloquante.

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 log2(nombre), les tests, faire
log2(nombre) se multiplie et divise. Depuis la division est une division par 2, vous pouvez le remplacer par un >>.

La condition d'arrêt peut ne pas être exactement les mêmes, alors assurez-vous de tester une variété de nombres entiers pour être sûr que la division par 2 n'a pas tort osciller entre deux valeurs; ils s'en distinguent par plus de 1.

Je trouve que la plupart des algorithmes sont basés sur des idées simples, mais sont mis en œuvre de manière plus compliqué que nécessaire. J'ai pris l'idée à partir d'ici: http://ww1.microchip.com/downloads/en/AppNotes/91040a.pdf (par Ross M. Fosler) et en a fait un très bref C-fonction:

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

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

J'ai réglé quelque chose de similaire à la binary digit par digit de l'algorithme décrit dans cet article de Wikipédia.

J'ai récemment rencontré la même tâche sur le processeur ARM Cortex-M3 (STM32F103CBT6) et après des recherches sur Internet est venu avec la solution suivante. Ce n'est pas le plus rapide en comparaison avec les solutions proposées ici, mais il a une bonne précision (de l'erreur maximale est de 1, c'est à dire LSB sur l'ensemble de la UI32 d'entrée de gamme) et de relativement bonne vitesse (environ 1,3 M racines carrées par seconde sur une de 72 MHz ARM Cortex-M3, soit environ 55 cycles par racine unique, y compris l'appel de la 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;
}

Je suis en utilisant IAR et elle produit les suivants du code assembleur:

        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

J'ai mis en place Warren de la suggestion et de la méthode de Newton en C# pour les entiers 64 bits. Isqrt utilise la méthode de Newton, tandis que Isqrt utilise Warren méthode. 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 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<<DeBruijnArray[x])*Math.Sqrt(2)));
        }

        public static long Isqrt2(this long n)
        {
          ulong num = (ulong) n; 
          ulong y;
          if(num==0)
            return (long)num;
          {
            ulong 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)>>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é les méthodes suivantes pour référence 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 Release, Visual Studio 2012 ont été
que l'appel de la Bibliothèque de Mathématiques.Sqrt est plus rapide.

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

Je ne suis pas à l'aise avec les directives du compilateur, de sorte qu'il peut être possible de régler le compilateur pour obtenir le nombre entier maths plus vite. 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.