Я ищу быстрый, целочисленный алгоритм, чтобы найти квадратный корень (целочисленную его часть) целого числа без знака. Код должен иметь отличную производительность на процессорах ARM Thumb 2. Это может быть язык ассемблера или код C.
Любые приветствия приветствуются.
Целочисленные квадратные корни от Джека У. Креншоу могут быть полезны в качестве другой ссылки.
В C Snippets Archive также реализована целочисленная реализация с квадратным корнем. Это выходит за пределы только целочисленного результата и вычисляет дополнительные дробные (фиксированные) биты ответа. (Обновление: к сожалению, архив фрагментов C теперь не функционирует. Ссылка указывает на веб-архив страницы.) Вот код из архива фрагментов C:
#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));
}
Я остановился на следующем коде. Это в основном из статьи статьи в области квадратичных методов вычислений. Но было изменено использование 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 32 - 1 равен 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;
}
Если точная точность не требуется, у меня есть быстрое приближение для вас, которое использует 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.
Итак, мы берем вход X = A × B, где A = 2 N и 1 ≤ B & lt; 2
Затем у нас есть таблица поиска для sqrt (2 N) и таблица поиска для sqrt (1 ≤ B & lt; 2). Мы сохраняем таблицу поиска для sqrt (2 N) как целое число, что может быть ошибкой (тестирование не выявило вредных последствий), и мы храним таблицу поиска для sqrt (1 ≤ B & lt; 2) как 15-битное исправленное -точка.
Мы знаем, что 1 ≤ sqrt (2 N) & lt; 65536, так что 16-разрядный, и мы знаем, что мы действительно можем только умножить 16-разрядный × 15-разрядный на ARM, не опасаясь репрессий, так что мы делаем.
С точки зрения реализации, while(t) {cnt++;t>>=1;} фактически является инструкцией подсчета начальных битов (CLB), так что если ваша версия чипсета имеет это, вы выигрываете! Кроме того, инструкцию сдвига было бы легко реализовать с помощью двунаправленного переключателя, если он у вас есть?
Здесь есть алгоритм Lg [N] для подсчета старшего установленного бита здесь.
С точки зрения магических чисел, для изменения размеров таблицы, магическое число для ftbl2 равно 32, хотя обратите внимание, что 6 (Lg [32] +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.
Я считаю, что большинство алгоритмов основаны на простых идеях, но реализованы более сложным способом, чем необходимо. Я взял идею отсюда: 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. Я считаю, что ваш скомпилированный код будет в целом быстрее, если вы будете использовать циклы for вместо циклов while, и вы получите дополнительное преимущество детерминированного времени (хотя это в некоторой степени зависит от того, как ваш компилятор оптимизирует оператор if).
Это зависит от использования функции 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;
Это не быстро, но это маленький и простой:
int isqrt(int n)
{
int b = 0;
while(n >= 0)
{
n = n - b;
b = b + 1;
n = n - b;
}
return b - 1;
}
Я решил что-то похожее на двоичный алгоритм с цифрой по цифре, описанный в этой статье в Википедии.
Вот решение на 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; удваивает точность в битах. Последняя строка исключает корни крыши, которые могут быть получены.
Этот метод похож на длинное деление: вы строите предположение для следующей цифры корня, выполняете вычитание и вводите цифру, если разница соответствует определенным критериям. С бинарной версией ваш единственный выбор для следующей цифры - 0 или 1, поэтому вы всегда угадываете 1, выполняете вычитание и вводите 1, если разница не отрицательная.
http://www.realitypixels.com/turk/opensource/index.html#FractSqrt
Недавно я столкнулся с той же задачей на ARM Cortex-M3 (STM32F103CBT6) и после поиска в интернете нашел следующее решение. Это не самое быстрое сравнение с предлагаемыми здесь решениями, но он имеет хорошую точность (максимальная погрешность равна 1, то есть LSB во всем входном диапазоне UI32) и относительно хорошую скорость (около 1,3 М квадратных корней в секунду на ARM Cortex с частотой 72 МГц). M3 или около 55 циклов на один корень, включая вызов функции).
// FastIntSqrt is based on Wikipedia article:
// https://en.wikipedia.org/wiki/Methods_of_computing_square_roots
// Which involves Newton 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;
}
Я использую IAR, и он производит следующий код на ассемблере:
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
Если вам это нужно только для процессоров ARM Thumb 2, библиотека CMSIS DSP от ARM - лучший выбор для вас. Она сделана людьми, разработавшими процессоры Thumb 2. Кто еще может победить это?
На самом деле вам даже не нужен алгоритм, а специальные аппаратные инструкции с квадратным корнем, такие как VSQRT. Компания ARM поддерживает реализации алгоритмов математики и DSP, оптимизированных для процессоров с поддержкой Thumb 2, пытаясь использовать свое аппаратное обеспечение, такое как VSQRT. Вы можете получить исходный код:
arm_sqrt_f32()arm_sqrt_q15.c/arm_sqrt_q31.c (q15 и q31 - это типы данных с фиксированной запятой, специализированные для ядра ARM DSP, которое часто поставляется с Thum 2-совместимыми процессорами.)Обратите внимание, что ARM также поддерживает скомпилированные двоичные файлы CMSIS DSP, которые гарантируют наилучшую возможную производительность для инструкций, специфичных для архитектуры ARM Thumb. Таким образом, вы должны рассмотреть возможность статического связывания их при использовании библиотеки. Вы можете получить двоичные файлы здесь.
Я реализовал предложение Уоррена и метод Ньютона в С# для 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 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 происходит быстрее.
Я не сообразителен с директивами компилятора, так что может быть возможно настроить компилятор, чтобы быстрее получить целочисленную математику. Очевидно, что метод сдвига битов очень близок к библиотеке. В системе без математического сопроцессора это будет очень быстро.
Я разработал 16-битный sqrt для сжатия RGB-гаммы. Распределяется в 3 разные таблицы, основанные на старших 8 битах. Недостатки: он использует около килобайта для таблиц, раунды непредсказуемы, если точный sqrt невозможен, и, в худшем случае, использует одиночное умножение (но только для очень немногих входных значений).
static uint8_t sqrt_50_256[] = {
114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,
133,134,135,136,137,138,139,140,141,142,143,143,144,145,146,147,148,149,150,
150,151,152,153,154,155,155,156,157,158,159,159,160,161,162,163,163,164,165,
166,167,167,168,169,170,170,171,172,173,173,174,175,175,176,177,178,178,179,
180,181,181,182,183,183,184,185,185,186,187,187,188,189,189,190,191,191,192,
193,193,194,195,195,196,197,197,198,199,199,200,201,201,202,203,203,204,204,
205,206,206,207,207,208,209,209,210,211,211,212,212,213,214,214,215,215,216,
217,217,218,218,219,219,220,221,221,222,222,223,223,224,225,225,226,226,227,
227,228,229,229,230,230,231,231,232,232,233,234,234,235,235,236,236,237,237,
238,238,239,239,240,241,241,242,242,243,243,244,244,245,245,246,246,247,247,
248,248,249,249,250,250,251,251,252,252,253,253,254,254,255,255
};
static uint8_t sqrt_0_10[] = {
1,2,3,3,4,4,5,5,5,6,6,6,7,7,7,7,8,8,8,8,9,9,9,9,9,10,10,10,10,10,11,11,11,
11,11,11,12,12,12,12,12,12,13,13,13,13,13,13,13,14,14,14,14,14,14,14,15,15,
15,15,15,15,15,15,16,16,16,16,16,16,16,16,17,17,17,17,17,17,17,17,17,18,18,
18,18,18,18,18,18,18,19,19,19,19,19,19,19,19,19,19,20,20,20,20,20,20,20,20,
20,20,21,21,21,21,21,21,21,21,21,21,21,22,22,22,22,22,22,22,22,22,22,22,23,
23,23,23,23,23,23,23,23,23,23,23,24,24,24,24,24,24,24,24,24,24,24,24,25,25,
25,25,25,25,25,25,25,25,25,25,25,26,26,26,26,26,26,26,26,26,26,26,26,26,27,
27,27,27,27,27,27,27,27,27,27,27,27,27,28,28,28,28,28,28,28,28,28,28,28,28,
28,28,29,29,29,29,29,29,29,29,29,29,29,29,29,29,29,30,30,30,30,30,30,30,30,
30,30,30,30,30,30,30,31,31,31,31,31,31,31,31,31,31,31,31,31,31,31,31,32,32,
32,32,32,32,32,32,32,32,32,32,32,32,32,32,33,33,33,33,33,33,33,33,33,33,33,
33,33,33,33,33,33,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,35,35,
35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,36,36,36,36,36,36,36,36,36,
36,36,36,36,36,36,36,36,36,37,37,37,37,37,37,37,37,37,37,37,37,37,37,37,37,
37,37,37,38,38,38,38,38,38,38,38,38,38,38,38,38,38,38,38,38,38,38,39,39,39,
39,39,39,39,39,39,39,39,39,39,39,39,39,39,39,39,39,40,40,40,40,40,40,40,40,
40,40,40,40,40,40,40,40,40,40,40,40,41,41,41,41,41,41,41,41,41,41,41,41,41,
41,41,41,41,41,41,41,41,42,42,42,42,42,42,42,42,42,42,42,42,42,42,42,42,42,
42,42,42,42,43,43,43,43,43,43,43,43,43,43,43,43,43,43,43,43,43,43,43,43,43,
43,44,44,44,44,44,44,44,44,44,44,44,44,44,44,44,44,44,44,44,44,44,44,45,45,
45,45,45,45,45,45,45,45,45,45,45,45,45,45,45,45,45,45,45,45,45,46,46,46,46,
46,46,46,46,46,46,46,46,46,46,46,46,46,46,46,46,46,46,46,47,47,47,47,47,47,
47,47,47,47,47,47,47,47,47,47,47,47,47,47,47,47,47,47,48,48,48,48,48,48,48,
48,48,48,48,48,48,48,48,48,48,48,48,48,48,48,48,48,49,49,49,49,49,49,49,49,
49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,50,50,50,50,50,50,50,50,
50,50,50,50,50,50,50,50,50,50,50,50,50,50,50,50,50,51,51,51,51,51,51,51,51,
51,51,51,51,51,51,51,51,51,51,51,51,51,51,51,51,51,51,52,52,52,52,52,52,52,
52,52,52,52,52,52,52,52,52,52,52,52,52,52,52,52,52,52,52,53,53
};
static uint8_t sqrt_11_49[] = {
54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,0,76,77,78,
0,79,80,81,82,83,0,84,85,86,0,87,88,89,0,90,0,91,92,93,0,94,0,95,96,97,0,98,0,
99,100,101,0,102,0,103,0,104,105,106,0,107,0,108,0,109,0,110,0,111,112,113
};
uint16_t isqrt16(uint16_t v) {
uint16_t a, b;
uint16_t h = v>>8;
if (h <= 10) return v ? sqrt_0_10[v>>2] : 0;
if (h >= 50) return sqrt_50_256[h-50];
h = (h-11)<<1;
a = sqrt_11_49[h];
b = sqrt_11_49[h+1];
if (!a) return b;
return b*b > v ? a : b;
}
Я сравнил его с sqrt на основе log2, используя clang __builtin_clz (который должен быть расширен до кода операции одной сборки), и библиотеку sqrtf, называемую используя (int)sqrtf((float)i). И получил довольно странные результаты:
$ gcc -O3 test.c -o test && ./test
isqrt16: 6.498955
sqrtf: 6.981861
log2_sqrt: 61.755873
Clang скомпилировал вызов sqrtf к инструкции sqrtss, что почти так же быстро, как эта таблица sqrt. Извлеченный урок: на x86 компилятор может обеспечить достаточно быстрый sqrt, который менее чем на 10% медленнее, чем вы сами можете придумать, тратить много времени или может быть в 10 раз быстрее, если вы используете некрасивое побитовое хаки. И все же sqrtss немного медленнее, чем пользовательская функция, поэтому, если вам действительно нужны эти 5%, вы можете получить их, а ARM, например, не имеет sqrtss, поэтому log2_sqrt не должен сильно отставать.
На x86, где доступен FPU, старый взлом Quake представляется самым быстрым способом вычисления целочисленных значений sqrt. Это в 2 раза быстрее, чем эта таблица или FPU sqrtss.