Ищете эффективный алгоритм квадратного корня целочисленного значения для ARM Thumb2

Question

Я ищу быстрый, целочисленный алгоритм для поиска квадратного корня (целочисленной его части) целых чисел без знака. Код должен иметь отличную производительность на процессорах ARM Thumb 2. Это может быть язык ассемблера или код C.

Любые подсказки приветствуются.

Answer 1

Один общий подход - деление пополам.

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 

Нечто подобное должно работать достаточно хорошо. Он выполняет тесты log2 (число), делая log2 (число) умножает и делит. Поскольку разделение является делением на 2, вы можете заменить его на >>.

Условие завершения может не быть точным, поэтому обязательно проверяйте множество целых чисел, чтобы убедиться, что деление на 2 не колеблется некорректно между двумя четными значениями; они будут отличаться более чем на 1.

Answer 2

Если точная точность не требуется, у меня есть быстрое приближение к вам, которое использует 260 байт памяти (вы можете вдвое уменьшить это, но не делать этого).

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

Вот код для создания таблицы:

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"); 

в диапазоне 1-> 2^20, максимальная погрешность составляет 11, а в диапазоне 1-> 2^30, то около 256. Вы можете использовать большие таблицы и минимизировать это. Стоит отметить, что ошибка всегда будет отрицательной - т. Е. Когда это будет неправильно, значение будет МЕНЬШЕ, чем правильное значение.

Вы можете преуспеть, следуя этому этапу очистки.

Идея достаточно проста: (ab)^0.5 = a^0.b * b^0.5.

Итак, мы берем вход Х = A * B, где А = 2^N и 1 < = В Тогда мы имеем LookupTable для SQRT (2^N), и LookupTable для SQRT (1 < = B < 2). Мы храним поисковик для sqrt (2^N) как целое, что может быть ошибкой (тестирование не показывает никаких негативных последствий), и мы храним поисковую таблицу для sqrt (1 < = B < 2) в 15 бит фиксированной точки.

Мы знаем, что 1 < = SQRT (2^N) < 65536, так что 16bit, и мы знаем, что мы можем на самом деле только умножать 16bitx15bit на ARM, не опасаясь репрессий, так это то, что мы делаем.

С точки зрения реализации, while (t) {cnt ++; t >> = 1;} - это команда с начальным битом (CLB), поэтому, если ваша версия чипсета имеет это, вы победы! Кроме того, инструкция сдвига будет легко реализована с помощью двунаправленного переключателя, если у вас его есть? Там в [N], алгоритм Lg для подсчета высокого набора битого here.

С точкой зрения магических чисел, для изменения размеров таблицы, магическое число для ftbl2 32, хотя отмечает, что (Lg [32] + 1) используется для смещения.

Answer 3

Я решил что-то похожее на двоичный цифровой алгоритм, описанный в this Wikipedia article.

Answer 4

Integer Square Roots от Jack W. Crenshaw может быть полезным в качестве другой ссылки.

C Snippets Archive также имеет integer square root implementation. Это выходит за пределы только целочисленного результата и вычисляет дополнительные дробные (фиксированные) биты ответа. (Update: к сожалению, архивировать отрывки C теперь несуществующей ссылка указывает на веб-архив страницы.). Вот код из C Snippets 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)); 
} 

я остановился на следующем коде. Это по существу от Wikipedia article on square-root computing methods. Но было изменено использование stdint.h типов uint32_t и т. Д. Строго говоря, тип возврата можно было бы изменить на 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; 
} 

Приятная вещь, я обнаружил, заключается в том, что довольно простая модификация может вернуть «округленный» ответ. Я нашел это полезным в определенном приложении для большей точности. Обратите внимание, что в этом случае возвращаемый тип должен быть uint32_t, потому что закругленный квадратный корень из 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; 
} 

Answer 5

Это не быстро, но это маленькое и простые:

int isqrt(int n) 
{ 
    int b = 0; 

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

    return b - 1; 
} 

Answer 6

Этот метод похож на столбик: вы строите предположение для следующей цифры корня, сделать вычитание, и введите если разница соответствует определенным критериям. С бинарной версией ваш единственный выбор для следующей цифры - 0 или 1, поэтому вы всегда угадываете 1, делаете вычитание и вводите 1, если разница не является отрицательной.

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

Answer 7

Я считаю, что большинство алгоритмов основаны на простых идеях, но реализованы таким образом, более сложным образом, чем это необходимо. Я принял эту идею здесь: http://ww1.microchip.com/downloads/en/AppNotes/91040a.pdf (Росс М. Фослер) и сделал это в очень короткий 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; 
} 

Это составляет до 5 циклов/бит на моем Blackfin. Я считаю, что ваш скомпилированный код будет быстрее, если вы будете использовать для циклов вместо циклов while, и вы получите дополнительное преимущество детерминированного времени (хотя это в некоторой степени зависит от того, как ваш компилятор оптимизирует оператор if.)

Answer 8

Это зависит от использования функции sqrt. Я часто использую некоторые приближенные, чтобы создавать быстрые версии. Например, когда мне нужно вычислить модуль вектора:

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

Я использую:

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

который может быть закодирован в 3-х или 4 инструкции как:

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

Answer 9

Вот решение в Java, которое объединяет целочисленные log_2 и метод Ньютона для создания алгоритма, свободного от цикла. Как недостаток, он нуждается в делении. Прокомментированные строки должны быть преобразованы в 64-битный алгоритм.

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

Как это работает: первая часть производит квадратный корень с точностью до трех бит. Линия [y = (y + num/y) >> 1;] удваивает точность в битах. Последняя линия устраняет корни крыши, которые могут быть сгенерированы.

Answer 10

Я реализовал предложение Уоррена и метод Ньютона в C# для 64-битных целых чисел. Isqrt использует метод Ньютона, а Isqrt использует метод Уоррена. Вот исходный код:

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 

    } 
} 

Я использовал следующий для сравнения код:

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

    } 
} 

Мои результаты на Dell Latitude E6540 в режиме выпуска, Visual Studio 2012 были что Math вызов библиотеки .Sqrt быстрее.

• 59 мс - Ньютон (Isqrt)
• 12 мс - Бит передач (Isqrt2)
• 5 мс - Math.Sqrt

Я не умный с директивами компилятора, так что она может можно настроить компилятор, чтобы быстрее получить математику с целым числом. Ясно, что подход с битовым сдвигом очень близок к библиотеке. В системе без математического сопроцессора это было бы очень быстро.