二进制32位整数快速平方根

无心人

FWAIT是 Causes the processor to check for and handle pending, unmasked, floating-point exceptions before proceeding.

无心人

在386时代 FWAIT是等待FST FSTP FISTP写入内存完毕的指令 我想，P4时代应该不必了

无心人

重新回到整数算法上来 假设bsr方法得到初始值 因为，初始值是个2^k类型的值 所以首次的牛顿迭代是无需乘法和除法的 现在的问题是，首次迭代的值和精确值的差距是多少 二次迭代是否能得到准确值？？ 如果不能得到 那不能得到精确值的输入是否有范围？ 范围多大？

无心人

网上找到的，似乎有用处 http://tech.ddvip.com/2008/08/121928805856259.html 为工作的需要，要在单片机上实现开根号的操作。目前开平方的方法大部分是用牛顿 　　迭代法。我在查了一些资料以后找到了一个比牛顿迭代法更加快速的方法。 　　1.原理 　　因为排版的原因，用pow(X,Y)表示X的Y次幂，用B[0]，B[1]，...，B[m-1]表示一个序列， 　　其中[x]为下标。 　　假设： 　　B[x],b[x]都是二进制序列,取值0或1。 　　M = B[m-1]*pow(2,m-1) + B[m-2]*pow(2,m-2) + ... + B[1]*pow(2,1) + B[0]*pow 　　(2,0) 　　N = b[n-1]*pow(2,n-1) + b[n-2]*pow(2,n-2) + ... + b[1]*pow(2,1) + n[0]*pow 　　(2,0) 　　pow(N,2) = M 　　(1) N的最高位b[n-1]可以根据M的最高位B[m-1]直接求得。 　　设 m 已知,因为 pow(2, m-1) <= M <= pow(2, m)，所以 pow(2, (m-1)/2) <= N <= 　　pow(2, m/2) 　　如果 m 是奇数，设m=2*k+1, 　　那么 pow(2,k) <= N < pow(2, 1/2+k) < pow(2, k+1), 　　n-1=k, n=k+1=(m+1)/2 　　如果 m 是偶数，设m=2k, 　　那么 pow(2,k) > N >= pow(2, k-1/2) > pow(2, k-1), 　　n-1=k-1,n=k=m/2 　　所以b[n-1]完全由B[m-1]决定。 　　余数 M[1] = M - b[n-1]*pow(2, 2*n-2) 　　(2) N的次高位b[n-2]可以采用试探法来确定。 　　因为b[n-1]=1，假设b[n-2]=1，则 pow(b[n-1]*pow(2,n-1) + b[n-1]*pow(2,n-2), 　　2) = b[n-1]*pow(2,2*n-2) + (b[n-1]*pow(2,2*n-2) + b[n-2]*pow(2,2*n-4)), 　　然后比较余数M[1]是否大于等于 (pow(2,2)*b[n-1] + b[n-2]) * pow(2,2*n-4)。这种 　　比较只须根据B[m-1]、B[m-2]、...、B[2*n-4]便可做出判断，其余低位不做比较。 　　若 M[1] >= (pow(2,2)*b[n-1] + b[n-2]) * pow(2,2*n-4), 则假设有效，b[n-2] = 　　1； 　　余数 M[2] = M[1] - pow(pow(2,n-1)*b[n-1] + pow(2,n-2)*b[n-2], 2) = M[1] - 　　(pow(2,2)+1)*pow(2,2*n-4)； 　　若 M[1] < (pow(2,2)*b[n-1] + b[n-2]) * pow(2,2*n-4), 则假设无效，b[n-2] = 　　0；余数 M[2] = M[1]。 　　(3) 同理，可以从高位到低位逐位求出M的平方根N的各位。 　　使用这种算法计算32位数的平方根时最多只须比较16次，而且每次比较时不必把M的各位逐 　　一比较，尤其是开始时比较的位数很少，所以消耗的时间远低于牛顿迭代法。 　　2. 流程图 　　(制作中，稍候再上) 　　3. 实现代码 　　这里给出实现32位无符号整数开方得到16位无符号整数的C语言代码。

1. /****************************************/
2. /*Function: 开根号处理　　　　　　　　　*/
3. /*入口参数：被开方数，长整型　　　　　　*/
4. /*出口参数：开方结果，整型　　　　　　　*/
5. /****************************************/
6. unsigned int sqrt_16(unsigned long M)
7. {
8. 　　unsigned int N, i;
9. 　　unsigned long tmp, ttp;　// 结果、循环计数
10. 　　if (M == 0)　　　　　　　// 被开方数，开方结果也为0
11. 　　　　return 0;
12. 　　N = 0;
13. 　　tmp = (M >> 30);　　　　　// 获取最高位：B[m-1]
14. 　　M <<= 2;
15. 　　if (tmp > 1)　　　　　　　// 最高位为1
16. 　　{
17. 　　　　N ++;　　　　　　　　// 结果当前位为1，否则为默认的0
18. 　　　　tmp -= N;
19. 　　}
20. 　　for (i=15; i>0; i--)　　　// 求剩余的15位
21. 　　{
22. 　　　　N <<= 1;　　　　　　　// 左移一位
23. 　　　　tmp <<= 2;
24. 　　　　tmp += (M >> 30);　　// 假设
25. 　　　　ttp = N;
26. 　　　　ttp = (ttp<<1)+1;
27. 　　　　M <<= 2;
28. 　　　　if (tmp >= ttp)　　　// 假设成立
29. 　　　　{
30. 　　　　　　tmp -= ttp;
31. 　　　　　　N ++;
32. 　　　　}
33. 　　}
34. 　　return N;
35. }

复制代码

无心人

/****************************************/ /*Function: 开根号处理　　　　　　　　　*/ /*入口参数：被开方数，长整型　　　　　　*/ /*出口参数：开方结果，整型　　　　　　　*/ /****************************************/ #include #include unsigned int sqrt_16(unsigned int M) { unsigned int N, i; unsigned int tmp, ttp; // 结果、循环计数 if (M == 0) // 被开方数，开方结果也为0 return 0; N = 0; tmp = (M >> 30); // 获取最高位：B[m-1] M <<= 2; if (tmp > 1) //最高位为1 { N ++; // 结果当前位为1，否则为默认的0 tmp -= N; } for (i=15; i>0; i--) // 求剩余的15位 { N <<= 1; // 左移一位 tmp <<= 2; tmp += (M >> 30); // 假设 ttp = N; ttp = (ttp<<1)+1; M <<= 2; if (tmp >= ttp) // 假设成立 { tmp -= ttp; N ++; } } return N; } int main(void) { unsigned int M, N; printf("Input M:"); scanf("%u", &M); N = sqrt_16(M); printf("\nSQRT(%u) = %u\n", M, N); return 0; } 写了个测试程序 测试了几个例子 E:\MinGW\MSYS\math>sqrt Input M:1023 SQRT(1023) = 31 E:\MinGW\MSYS\math>sqrt Input M:4000000000 SQRT(4000000000) = 63245 E:\MinGW\MSYS\math>sqrt Input M:0 SQRT(0) = 0 E:\MinGW\MSYS\math>sqrt Input M:1 SQRT(1) = 1 E:\MinGW\MSYS\math>sqrt Input M:23 SQRT(23) = 4 E:\MinGW\MSYS\math>sqrt Input M:400000000 SQRT(400000000) = 20000

gxqcn

很奇怪的现象： 虽说我的43#代码够简洁， 但测试结果速度并不理想，基本上比无心人和liangbch写的函数要慢一倍。 我又测试了HugeCalc内部采用的纯C/C++版本（不含浮点计算）， 虽然语句比较多，还有跳转语句（不见得编译后就有跳转指令）， 但测试起来速度反而比用fsqrt的都快很多很多！ 不信大家可以测试测试：

1. UINT32 Sqrt( UINT32 n )
2. {
3. UINT32 t;
4. UINT32 ret = 0;
6. #define SQRT_CORE(x) \
7. t = ret + (1UL<<(x)); \
8. ret >>= 1; \
9. if ( n >= t ) \
10. { \
11. n -= t; \
12. ret += (1UL<<(x)); \
13. }
15. SQRT_CORE(30);
16. SQRT_CORE(28);
17. SQRT_CORE(26);
18. SQRT_CORE(24);
19. SQRT_CORE(22);
20. SQRT_CORE(20);
21. SQRT_CORE(18);
22. SQRT_CORE(16);
23. SQRT_CORE(14);
24. SQRT_CORE(12);
25. SQRT_CORE(10);
26. SQRT_CORE(8);
27. SQRT_CORE(6);
28. SQRT_CORE(4);
29. SQRT_CORE(2);
30. SQRT_CORE(0);
32. #undef SQRT_CORE
34. return ret;
35. }

复制代码

无心人

http://www.dianyuan.com/bbs/d/50/143623.html 这个也许有点用处

无心人

56#的似乎类似于54#的算法 呵呵 因为不存在乘法和除法 只有逻辑移位和比较 所以比较快 =================== fsqrt的大概58-60个时钟 郭老大的代码最少32个时钟吧 有时间依照54#的逻辑写个汇编 呵呵

liangbch

他这个算法是固定15次循环，我想这样做。 1.做一个256 byte 的表。 2. 如果被开放数小于8bit 直接查表得到结果。 否则如果被开方数是k(<8k<=32)bit，则首先通过一次得到最高结果的最高4bit，然后循环 (k-8)/2次循环得到结果的其他bit. 3.程序仍然采用先前的结果，bsr 加跳表。

无心人

测试了 FPU2， 54#， 56#代码 2100000000 -- 2200000000 FPU2 Version: 3587.139 ms iSqrt_16 Version: 26195.432 ms iSqrt_GxQ Version: 16059.504 ms