我正在寻找一种快速、仅整数的算法来查找无符号整数的平方根（其整数部分）。代码必须在 ARM Thumb 2 处理器上具有出色的性能。它可以是汇编语言或 C 代码。
欢迎任何提示。
问问题
64793 次
14 回答
36
Jack W. Crenshaw 的Integer Square Roots可以作为另一个参考。
C Snippets Archive也有一个整数平方根实现。这不仅仅是整数结果，还计算答案的额外小数（定点）位。（更新：不幸的是，C 片段存档现已失效。链接指向页面的 Web 存档。）这是来自 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;
}
我发现，好消息是一个相当简单的修改可以返回“四舍五入”的答案。我发现这在某些应用程序中很有用，可以提高准确性。请注意，在这种情况下，返回类型必须是因为 2 32 - 1uint32_t的四舍五入的平方根是 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;
}
于 2009-07-09T00:12:53.897 回答
16
如果不需要精确的精度，我有一个快速的近似值，它使用 260 字节的 RAM（你可以减半，但不要）。
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 < 2
然后我们有一个用于 sqrt(2 N ) 的查找表和一个用于 sqrt(1 ≤ B < 2) 的查找表。我们将 sqrt(2 N ) 的查找表存储为整数，这可能是一个错误（测试显示没有不良影响），我们将 sqrt(1 ≤ B < 2) 的查找表存储为 15 位定点。
我们知道 1 ≤ sqrt(2 N ) < 65536，所以这是 16 位，而且我们知道我们实际上只能在 ARM 上乘以 16 位 × 15 位，而不必担心遭到报复，所以这就是我们所做的。
在实现方面，while(t) {cnt++;t>>=1;}它实际上是一个计数前导位指令 (CLB)，因此，如果您的芯片组版本具有该指令，那么您就赢了！此外，如果你有一个双向移位器，移位指令很容易实现？
这里有一个用于计算最高设置位的 Lg[N] 算法。
就幻数而言，对于更改表大小，幻数ftbl2是 32，但请注意，6 (Lg[32]+1) 用于移位。
于 2009-07-08T21:19:36.240 回答
9
一种常见的方法是二分法。
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(number) 测试，进行 log2(number) 乘法和除法。由于除法是除以 2，因此您可以将其替换为>>.
终止条件可能不准确，所以一定要测试各种整数，以确保除以 2 不会在两个偶数之间错误地振荡；它们的差异将超过 1。
于 2009-07-08T19:35:02.363 回答
8
I find that most algorithms are based on simple ideas, but are implemented in a way more complicated manner than necessary. I've taken the idea from here: http://ww1.microchip.com/downloads/en/AppNotes/91040a.pdf (by Ross M. Fosler) and made it into a very short C-function:
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;
}
This compiles to 5 cycles/bit on my blackfin. I believe your compiled code will in general be faster if you use for loops instead of while loops, and you get the added benefit of deterministic time (although that to some extent depends on how your compiler optimizes the if statement.)
于 2012-04-26T09:45:13.140 回答
8
它并不快，但它又小又简单：
int isqrt(int n)
{
int b = 0;
while(n >= 0)
{
n = n - b;
b = b + 1;
n = n - b;
}
return b - 1;
}
于 2009-08-27T20:54:19.487 回答
7
这取决于 sqrt 函数的用法。我经常使用一些近似来制作快速版本。例如，当我需要计算 vector 的模块时：
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;
于 2012-07-04T13:37:38.077 回答
4
我已经习惯了类似于这篇维基百科文章中描述的二进制逐位算法。
于 2009-07-08T23:32:54.553 回答
2
ARM 的最巧妙编码的按位整数平方根实现实现了每个结果位 3 个周期，对于 32 位无符号整数的平方根，其下限为 50 个周期。一个示例显示在 Andrew N. Sloss、Dominic Symes、Chris Wright，“ARM 系统开发人员指南”，Morgan Kaufman 2004 中。
由于大多数 ARM 处理器还具有非常快速的整数乘法器，并且大多数甚至提供了非常快速的宽乘法指令UMULL实现，因此可以实现大约 35 到 45 个周期的执行时间的另一种方法是通过平方根的倒数 1/ √x 使用定点计算。为此，有必要借助 count-leading-zeros 指令对输入进行规范化，该指令在大多数 ARM 处理器上都可用作指令CLZ。
计算从一个查找表的初始低精度倒数平方根近似值开始，该查找表由归一化参数的一些最高有效位索引。Newton-Raphson 迭代细化具有二次收敛性r的数的倒数平方根为 r n+1 = r n + r n * (1 - a * r n 2 ) / 2。这可以重新排列为代数等价形式方便。在下面的示例 C99 代码中，从 96 项查找表中读取8 位近似值 r 0 。这个近似值精确到大约 7 位。第一次 Newton-Raphson 迭代计算 r 1 = (3 * r 0 - a * r 0a3 ) / 2 以潜在地利用小操作数乘法指令。第二次 Newton-Raphson 迭代计算 r 2 = (r 1 * (3 - r 1 * (r 1 * a))) / 2。
然后通过反向乘法 s 2 = a * r 2计算归一化的平方根，并通过基于原始参数的前导零计数的非归一化来实现最终近似a。重要的是，期望的结果⌊√a⌋被低估了。这通过保证余数 ⌊√a⌋ - s 2 * s 2是正数来简化检查是否已达到预期结果。如果发现最终的近似值太小，则将结果增加一。通过针对“黄金”参考对所有可能的 2 32 个输入进行详尽的测试，可以很容易地证明该算法的正确操作，这只需要几分钟。
可以通过预先计算 3 * r 0和 r 0 3来简化第一次 Newton-Raphson 迭代，以牺牲查找表的额外存储空间为代价来加速此实现。前者需要 10 位存储，后者需要 24 位。为了将每一对组合成一个 32 位的数据项，立方体被四舍五入到 22 位，这在计算中引入了可忽略的误差。这会产生一个 96 * 4 = 384 字节的查找表。
另一种方法是观察到所有起始近似值都具有相同的最高有效位集，因此可以隐式假设并且不必存储。这允许将 9 位近似值 r 0压缩到 8 位数据项中，并在表查找后恢复前导位。使用包含 384 个条目的查找表，所有这些都被低估了，可以达到大约 7.5 位的准确度。将反向乘法与倒数平方根的 Newton-Raphson 迭代相结合，计算 s 0 = a * r 0， s 1 = s 0 + r 0 * (a - s 0 * s 0) / 2. 因为起始逼近的精度对于非常精确的最终平方根逼近来说不够高，所以最多可以偏离三倍，并且必须根据余数下限（sqrt (a)) - s 1 * s 1。
另一种方法的一个优点是所需的乘法次数减半，特别是只需要一次宽乘法UMULL。尤其是宽乘相当慢的处理器，这是一个值得尝试的替代方案。
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <math.h>
#if defined(_MSC_VER) && defined(_WIN64)
#include <intrin.h>
#endif // defined(_MSC_VER) && defined(_WIN64)
#define CLZ_BUILTIN (1) // use compiler's built-in count-leading-zeros
#define CLZ_FPU (2) // emulate count-leading-zeros via FPU
#define CLZ_CPU (3) // emulate count-leading-zeros via CPU
#define ALT_IMPL (0) // alternative implementation with fewer multiplies
#define LARGE_TABLE (0) // ALT_IMPL=0: incorporate 1st NR-iter into table
#define CLZ_IMPL (CLZ_CPU)// choose count-leading-zeros implementation
#define GEN_TAB (0) // generate tables
uint32_t umul32_hi (uint32_t a, uint32_t b);
uint32_t float_as_uint32 (float a);
int clz32 (uint32_t a);
#if ALT_IMPL
uint8_t rsqrt_tab [384] =
{
0xfe, 0xfc, 0xfa, 0xf8, 0xf6, 0xf4, 0xf2, 0xf0, 0xee, 0xed, 0xeb, 0xe9,
0xe7, 0xe6, 0xe4, 0xe2, 0xe1, 0xdf, 0xdd, 0xdc, 0xda, 0xd8, 0xd7, 0xd5,
0xd4, 0xd2, 0xd1, 0xcf, 0xce, 0xcc, 0xcb, 0xc9, 0xc8, 0xc7, 0xc5, 0xc4,
0xc2, 0xc1, 0xc0, 0xbe, 0xbd, 0xbc, 0xba, 0xb9, 0xb8, 0xb7, 0xb5, 0xb4,
0xb3, 0xb2, 0xb0, 0xaf, 0xae, 0xad, 0xac, 0xab, 0xa9, 0xa8, 0xa7, 0xa6,
0xa5, 0xa4, 0xa3, 0xa2, 0xa0, 0x9f, 0x9e, 0x9d, 0x9c, 0x9b, 0x9a, 0x99,
0x98, 0x97, 0x96, 0x95, 0x94, 0x93, 0x92, 0x91, 0x90, 0x8f, 0x8e, 0x8d,
0x8c, 0x8b, 0x8b, 0x8a, 0x89, 0x88, 0x87, 0x86, 0x85, 0x84, 0x83, 0x83,
0x82, 0x81, 0x80, 0x7f, 0x7e, 0x7d, 0x7d, 0x7c, 0x7b, 0x7a, 0x79, 0x79,
0x78, 0x77, 0x76, 0x75, 0x75, 0x74, 0x73, 0x72, 0x72, 0x71, 0x70, 0x6f,
0x6f, 0x6e, 0x6d, 0x6c, 0x6c, 0x6b, 0x6a, 0x6a, 0x69, 0x68, 0x67, 0x67,
0x66, 0x65, 0x65, 0x64, 0x63, 0x63, 0x62, 0x61, 0x61, 0x60, 0x5f, 0x5f,
0x5e, 0x5d, 0x5d, 0x5c, 0x5c, 0x5b, 0x5a, 0x5a, 0x59, 0x58, 0x58, 0x57,
0x57, 0x56, 0x55, 0x55, 0x54, 0x54, 0x53, 0x52, 0x52, 0x51, 0x51, 0x50,
0x50, 0x4f, 0x4e, 0x4e, 0x4d, 0x4d, 0x4c, 0x4c, 0x4b, 0x4b, 0x4a, 0x4a,
0x49, 0x48, 0x48, 0x47, 0x47, 0x46, 0x46, 0x45, 0x45, 0x44, 0x44, 0x43,
0x43, 0x42, 0x42, 0x41, 0x41, 0x40, 0x40, 0x3f, 0x3f, 0x3e, 0x3e, 0x3d,
0x3d, 0x3c, 0x3c, 0x3c, 0x3b, 0x3b, 0x3a, 0x3a, 0x39, 0x39, 0x38, 0x38,
0x37, 0x37, 0x36, 0x36, 0x36, 0x35, 0x35, 0x34, 0x34, 0x33, 0x33, 0x33,
0x32, 0x32, 0x31, 0x31, 0x30, 0x30, 0x30, 0x2f, 0x2f, 0x2e, 0x2e, 0x2d,
0x2d, 0x2d, 0x2c, 0x2c, 0x2b, 0x2b, 0x2b, 0x2a, 0x2a, 0x29, 0x29, 0x29,
0x28, 0x28, 0x27, 0x27, 0x27, 0x26, 0x26, 0x26, 0x25, 0x25, 0x24, 0x24,
0x24, 0x23, 0x23, 0x23, 0x22, 0x22, 0x21, 0x21, 0x21, 0x20, 0x20, 0x20,
0x1f, 0x1f, 0x1f, 0x1e, 0x1e, 0x1e, 0x1d, 0x1d, 0x1d, 0x1c, 0x1c, 0x1c,
0x1b, 0x1b, 0x1a, 0x1a, 0x1a, 0x19, 0x19, 0x19, 0x18, 0x18, 0x18, 0x17,
0x17, 0x17, 0x17, 0x16, 0x16, 0x16, 0x15, 0x15, 0x15, 0x14, 0x14, 0x14,
0x13, 0x13, 0x13, 0x12, 0x12, 0x12, 0x11, 0x11, 0x11, 0x11, 0x10, 0x10,
0x10, 0x0f, 0x0f, 0x0f, 0x0e, 0x0e, 0x0e, 0x0e, 0x0d, 0x0d, 0x0d, 0x0c,
0x0c, 0x0c, 0x0c, 0x0b, 0x0b, 0x0b, 0x0a, 0x0a, 0x0a, 0x0a, 0x09, 0x09,
0x09, 0x08, 0x08, 0x08, 0x08, 0x07, 0x07, 0x07, 0x07, 0x06, 0x06, 0x06,
0x05, 0x05, 0x05, 0x05, 0x04, 0x04, 0x04, 0x04, 0x03, 0x03, 0x03, 0x03,
0x02, 0x02, 0x02, 0x02, 0x01, 0x01, 0x01, 0x01, 0x00, 0x00, 0x00, 0x00,
};
/* compute floor (sqrt (a)) */
uint32_t my_isqrt32 (uint32_t a)
{
uint32_t b, r, s, scal, rem;
if (a == 0) return a;
/* Normalize argument */
scal = clz32 (a) & ~1;
b = a << scal;
/* Compute initial approximation to 1/sqrt(a) */
r = rsqrt_tab [(b >> 23) - 128] | 0x100;
/* Compute initial approximation to sqrt(a) */
s = umul32_hi (b, r << 8);
/* Refine sqrt approximation */
b = b - s * s;
s = s + ((r * (b >> 2)) >> 23);
/* Denormalize result*/
s = s >> (scal >> 1);
/* Ensure we got the floor correct */
rem = a - s * s;
if (rem < (2 * s + 1)) s = s + 0;
else if (rem < (4 * s + 4)) s = s + 1;
else if (rem < (6 * s + 9)) s = s + 2;
else s = s + 3;
return s;
}
#else // ALT_IMPL
#if LARGE_TABLE
uint32_t rsqrt_tab [96] =
{
0xfa0bfafa, 0xee6b2aee, 0xe5f02ae5, 0xdaf26ed9, 0xd2f002d0, 0xc890c2c4,
0xc1037abb, 0xb9a75ab2, 0xb4da42ac, 0xadcea2a3, 0xa6f27a9a, 0xa279c294,
0x9beb4a8b, 0x97a5ca85, 0x9163427c, 0x8d4fca76, 0x89500270, 0x8563ba6a,
0x818ac264, 0x7dc4ea5e, 0x7a120258, 0x7671da52, 0x72e4424c, 0x6f690a46,
0x6db24243, 0x6a52423d, 0x67042637, 0x6563c234, 0x62302a2e, 0x609cea2b,
0x5d836a25, 0x5bfd1a22, 0x58fd421c, 0x5783ae19, 0x560e4a16, 0x53300210,
0x51c7120d, 0x50623a0a, 0x4da4c204, 0x4c4c1601, 0x4af769fe, 0x49a6b9fb,
0x485a01f8, 0x471139f5, 0x45cc59f2, 0x448b5def, 0x4214fde9, 0x40df89e6,
0x3fade1e3, 0x3e8001e0, 0x3d55e1dd, 0x3c2f79da, 0x3c2f79da, 0x3b0cc5d7,
0x39edc1d4, 0x38d265d1, 0x37baa9ce, 0x36a689cb, 0x359601c8, 0x348909c5,
0x348909c5, 0x337f99c2, 0x3279adbf, 0x317741bc, 0x30784db9, 0x30784db9,
0x2f7cc9b6, 0x2e84b1b3, 0x2d9001b0, 0x2d9001b0, 0x2c9eb1ad, 0x2bb0b9aa,
0x2bb0b9aa, 0x2ac615a7, 0x29dec1a4, 0x29dec1a4, 0x28fab5a1, 0x2819e99e,
0x2819e99e, 0x273c599b, 0x273c599b, 0x26620198, 0x258ad995, 0x258ad995,
0x24b6d992, 0x24b6d992, 0x23e5fd8f, 0x2318418c, 0x2318418c, 0x224d9d89,
0x224d9d89, 0x21860986, 0x21860986, 0x20c18183, 0x20c18183, 0x20000180,
};
#else // LARGE_TABLE
uint8_t rsqrt_tab [96] =
{
0xfe, 0xfa, 0xf7, 0xf3, 0xf0, 0xec, 0xe9, 0xe6, 0xe4, 0xe1, 0xde, 0xdc,
0xd9, 0xd7, 0xd4, 0xd2, 0xd0, 0xce, 0xcc, 0xca, 0xc8, 0xc6, 0xc4, 0xc2,
0xc1, 0xbf, 0xbd, 0xbc, 0xba, 0xb9, 0xb7, 0xb6, 0xb4, 0xb3, 0xb2, 0xb0,
0xaf, 0xae, 0xac, 0xab, 0xaa, 0xa9, 0xa8, 0xa7, 0xa6, 0xa5, 0xa3, 0xa2,
0xa1, 0xa0, 0x9f, 0x9e, 0x9e, 0x9d, 0x9c, 0x9b, 0x9a, 0x99, 0x98, 0x97,
0x97, 0x96, 0x95, 0x94, 0x93, 0x93, 0x92, 0x91, 0x90, 0x90, 0x8f, 0x8e,
0x8e, 0x8d, 0x8c, 0x8c, 0x8b, 0x8a, 0x8a, 0x89, 0x89, 0x88, 0x87, 0x87,
0x86, 0x86, 0x85, 0x84, 0x84, 0x83, 0x83, 0x82, 0x82, 0x81, 0x81, 0x80,
};
#endif //LARGE_TABLE
/* compute floor (sqrt (a)) */
uint32_t my_isqrt32 (uint32_t a)
{
uint32_t b, r, s, t, scal, rem;
if (a == 0) return a;
/* Normalize argument */
scal = clz32 (a) & ~1;
b = a << scal;
/* Initial approximation to 1/sqrt(a)*/
t = rsqrt_tab [(b >> 25) - 32];
/* First NR iteration */
#if LARGE_TABLE
r = (t << 22) - umul32_hi (b, t);
#else // LARGE_TABLE
r = ((3 * t) << 22) - umul32_hi (b, (t * t * t) << 8);
#endif // LARGE_TABLE
/* Second NR iteration */
s = umul32_hi (r, b);
s = 0x30000000 - umul32_hi (r, s);
r = umul32_hi (r, s);
/* Compute sqrt(a) = a * 1/sqrt(a). Adjust to ensure it's an underestimate*/
r = umul32_hi (r, b) - 1;
/* Denormalize result */
r = r >> ((scal >> 1) + 11);
/* Make sure we got the floor correct */
rem = a - r * r;
if (rem >= (2 * r + 1)) r++;
return r;
}
#endif // ALT_IMPL
uint32_t umul32_hi (uint32_t a, uint32_t b)
{
return (uint32_t)(((uint64_t)a * b) >> 32);
}
uint32_t float_as_uint32 (float a)
{
uint32_t r;
memcpy (&r, &a, sizeof r);
return r;
}
int clz32 (uint32_t a)
{
#if (CLZ_IMPL == CLZ_FPU)
// Henry S. Warren, Jr, " Hacker's Delight 2nd ed", p. 105
int n = 158 - (float_as_uint32 ((float)(int32_t)(a & ~(a >> 1))+.5f) >> 23);
return (n < 0) ? 0 : n;
#elif (CLZ_IMPL == CLZ_CPU)
static const uint8_t clz_tab[32] = {
31, 22, 30, 21, 18, 10, 29, 2, 20, 17, 15, 13, 9, 6, 28, 1,
23, 19, 11, 3, 16, 14, 7, 24, 12, 4, 8, 25, 5, 26, 27, 0
};
a |= a >> 16;
a |= a >> 8;
a |= a >> 4;
a |= a >> 2;
a |= a >> 1;
return clz_tab [0x07c4acddu * a >> 27] + (!a);
#else // CLZ_IMPL == CLZ_BUILTIN
#if defined(_MSC_VER) && defined(_WIN64)
return __lzcnt (a);
#else // defined(_MSC_VER) && defined(_WIN64)
return __builtin_clz (a);
#endif // defined(_MSC_VER) && defined(_WIN64)
#endif // CLZ_IMPL
}
/* Henry S. Warren, Jr., "Hacker's Delight, 2nd e.d", p. 286 */
uint32_t ref_isqrt32 (uint32_t x)
{
uint32_t m, y, b;
m = 0x40000000U;
y = 0U;
while (m != 0) {
b = y | m;
y = y >> 1;
if (x >= b) {
x = x - b;
y = y | m;
}
m = m >> 2;
}
return y;
}
#if defined(_WIN32)
#if !defined(WIN32_LEAN_AND_MEAN)
#define WIN32_LEAN_AND_MEAN
#endif
#include <windows.h>
double second (void)
{
LARGE_INTEGER t;
static double oofreq;
static int checkedForHighResTimer;
static BOOL hasHighResTimer;
if (!checkedForHighResTimer) {
hasHighResTimer = QueryPerformanceFrequency (&t);
oofreq = 1.0 / (double)t.QuadPart;
checkedForHighResTimer = 1;
}
if (hasHighResTimer) {
QueryPerformanceCounter (&t);
return (double)t.QuadPart * oofreq;
} else {
return (double)GetTickCount() * 1.0e-3;
}
}
#elif defined(__linux__) || defined(__APPLE__)
#include <stddef.h>
#include <sys/time.h>
double second (void)
{
struct timeval tv;
gettimeofday(&tv, NULL);
return (double)tv.tv_sec + (double)tv.tv_usec * 1.0e-6;
}
#else
#error unsupported platform
#endif
int main (void)
{
#if ALT_IMPL
printf ("Alternate integer square root implementation\n");
#else // ALT_IMPL
#if LARGE_TABLE
printf ("Integer square root implementation w/ large table\n");
#else // LARGE_TAB
printf ("Integer square root implementation w/ small table\n");
#endif
#endif // ALT_IMPL
#if GEN_TAB
printf ("Generating lookup table ...\n");
#if ALT_IMPL
for (int i = 0; i < 384; i++) {
double x = 1.0 + (i + 1) * 1.0 / 128;
double y = 1.0 / sqrt (x);
uint8_t val = (uint8_t)((y * 512) - 256);
rsqrt_tab[i] = val;
printf ("0x%02x, ", rsqrt_tab[i]);
if (i % 12 == 11) printf("\n");
}
#else // ALT_IMPL
for (int i = 0; i < 96; i++ ) {
double x1 = 1.0 + i * 1.0 / 32;
double x2 = 1.0 + (i + 1) * 1.0 / 32;
double y = (1.0 / sqrt(x1) + 1.0 / sqrt(x2)) * 0.5;
uint32_t val = (uint32_t)(y * 256 + 0.5);
#if LARGE_TABLE
uint32_t cube = val * val * val;
rsqrt_tab[i] = (((cube + 1) / 4) << 10) + (3 * val);
printf ("0x%08x, ", rsqrt_tab[i]);
if (i % 6 == 5) printf ("\n");
#else // LARGE_TABLE
rsqrt_tab[i] = val;
printf ("0x%02x, ", rsqrt_tab[i]);
if (i % 12 == 11) printf ("\n");
#endif // LARGE_TABLE
}
#endif // ALT_IMPL
#endif // GEN_TAB
printf ("Running exhaustive test ... ");
uint32_t i = 0;
do {
uint32_t ref = ref_isqrt32 (i);
uint32_t res = my_isqrt32 (i);
if (res != ref) {
printf ("error: arg=%08x res=%08x ref=%08x\n", i, res, ref);
return EXIT_FAILURE;
}
i++;
} while (i);
printf ("PASSED\n");
printf ("Running benchmark ...\n");
i = 0;
uint32_t sum[8] = {0, 0, 0, 0, 0, 0, 0, 0};
double start = second();
do {
sum [0] += my_isqrt32 (i + 0);
sum [1] += my_isqrt32 (i + 1);
sum [2] += my_isqrt32 (i + 2);
sum [3] += my_isqrt32 (i + 3);
sum [4] += my_isqrt32 (i + 4);
sum [5] += my_isqrt32 (i + 5);
sum [6] += my_isqrt32 (i + 6);
sum [7] += my_isqrt32 (i + 7);
i += 8;
} while (i);
double stop = second();
printf ("%08x \relapsed=%.5f sec\n",
sum[0]+sum[1]+sum[2]+sum[3]+sum[4]+sum[5]+sum[6]+sum[7],
stop - start);
return EXIT_SUCCESS;
}
于 2021-12-19T19:51:40.023 回答
2
我最近在 ARM Cortex-M3 (STM32F103CBT6) 上遇到了同样的任务，在网上搜索后想出了以下解决方案。与此处提供的解决方案相比，它不是最快的，但它具有良好的精度（最大误差为 1，即整个 UI32 输入范围的 LSB）和相对较好的速度（在 72-MHz ARM Cortex-上大约每秒 1.3M 平方根- M3 或每个单根大约 55 个周期，包括函数调用）。
// 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;
}
我正在使用 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
于 2018-07-29T22:54:27.247 回答
1
这是一个结合整数 log_2 和牛顿方法来创建无循环算法的 Java 解决方案。不利的一面是，它需要划分。注释行需要上转换为 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;位精度加倍。最后一行消除了可以生成的屋顶根部。
于 2014-02-07T11:44:50.947 回答
1
如果您只需要 ARM Thumb 2 处理器，那么 ARM 的 CMSIS DSP 库是您的最佳选择。它是由设计 Thumb 2 处理器的人制作的。还有谁能打败它？
实际上，您甚至不需要算法，而是需要专门的平方根硬件指令，例如VSQRT。这家 ARM 公司通过尝试使用 VSQRT 等硬件来维护针对 Thumb 2 支持的处理器高度优化的数学和 DSP 算法实现。您可以获取源代码：
arm_sqrt_f32()arm_sqrt_q15.c/arm_sqrt_q31.c（q15 和 q31 是 ARM DSP 内核专用的定点数据类型，通常与 Thum 2 兼容处理器一起提供。）
请注意，ARM 还维护 CMSIS DSP 的编译二进制文件，以保证 ARM Thumb 体系结构特定指令的最佳性能。所以你应该在使用库时考虑静态链接它们。您可以在此处获取二进制文件。
于 2019-08-23T13:33:59.070 回答
0
此方法类似于长除法：您对根的下一位数字进行猜测，进行减法，如果差值符合特定标准，则输入该数字。对于二进制版本，下一个数字的唯一选择是 0 或 1，所以你总是猜 1，做减法，然后输入 1，除非差值为负。
http://www.realitypixels.com/turk/opensource/index.html#FractSqrt
于 2011-09-13T18:34:36.503 回答
-1
我在 C# 中为 64 位整数实现了 Warren 的建议和 Newton 方法。Isqrt 使用牛顿法，而 Isqrt 使用 Warren 法。这是源代码：
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
我对编译器指令并不聪明，因此可以调整编译器以更快地获得整数数学。显然，移位方法非常接近库。在没有数学协处理器的系统上，它会非常快。
于 2015-05-09T05:27:18.600 回答
-2
我为 RGB 伽玛压缩设计了一个 16 位 sqrt。它根据高 8 位分派到 3 个不同的表中。缺点：它使用大约一千字节的表，如果不可能精确的 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;
}
我已经将它与基于 log2 的 sqrt 进行了比较，使用 clang __builtin_clz（它应该扩展为单个程序集操作码）和库的sqrtf，称为 using (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%，浪费了大量时间，或者如果你使用一些丑陋的按位 hack，可以快 10 倍。而且仍然sqrtss比自定义函数慢一点，所以如果你真的需要这 5%，你可以得到它们，例如 ARM 没有sqrtss，所以 log2_sqrt 应该不会落后那么糟糕。
在可以使用 FPU 的 x86 上，旧的 Quake hack似乎是计算整数 sqrt 的最快方法。它比这张表或 FPU 的 sqrtss 快 2 倍。
于 2019-07-31T14:46:08.910 回答