数学研发网设为首页收藏本站
切换到窄版

 找回密码
 欢迎注册
数学研发论坛»论坛 › 【算法研发】 › 编程擂台 › 二进制32位整数快速平方根
12345678910... 15下一页
返回列表 发新帖
楼主: 无心人

[讨论] 二进制32位整数快速平方根

[复制链接]
liangbch
发表于 2009-2-10 21:20:10 | 显示全部楼层
下面给出实现这个功能的代码和测试代码。 顺便说一下，在2楼长提到的各个算法的选择都是在实验的基础上得出的，既能保证正确性，又是的速度不致太慢。
  1. #include <math.h>
  3. typedef unsigned long DWORD;
  4. typedef unsigned char BYTE;
  6. inline DWORD log2(DWORD n)
  7. {
  8. _asm
  9. {
  10. mov ecx,n
  11. bsr eax,ecx
  12. }
  13. }
  15. extern "C" unsigned char sqrtTab[]=
  16. {
  17. 0, 1, 1, 1, 2, 2, 2, 2,
  18. 2, 3, 3, 3, 3, 3, 3, 3,
  19. 4, 4, 4, 4, 4, 4, 4, 4,
  20. 4, 5, 5, 5, 5, 5, 5, 5,
  21. 5, 5, 5, 5, 6, 6, 6, 6,
  22. 6, 6, 6, 6, 6, 6, 6, 6,
  23. 6, 7, 7, 7, 7, 7, 7, 7,
  24. 7, 7, 7, 7, 7, 7, 7, 7,
  25. 128,128,129,130,131,132,133,134,
  26. 135,136,137,138,139,140,141,142,
  27. 143,144,144,145,146,147,148,149,
  28. 150,150,151,152,153,154,155,155,
  29. 156,157,158,159,160,160,161,162,
  30. 163,163,164,165,166,167,167,168,
  31. 169,170,170,171,172,173,173,174,
  32. 175,176,176,177,178,178,179,180,
  33. 181,181,182,183,183,184,185,185,
  34. 186,187,187,188,189,189,190,191,
  35. 192,192,193,193,194,195,195,196,
  36. 197,197,198,199,199,200,201,201,
  37. 202,203,203,204,204,205,206,206,
  38. 207,208,208,209,209,210,211,211,
  39. 212,212,213,214,214,215,215,216,
  40. 217,217,218,218,219,219,220,221,
  41. 221,222,222,223,224,224,225,225,
  42. 226,226,227,227,228,229,229,230,
  43. 230,231,231,232,232,233,234,234,
  44. 235,235,236,236,237,237,238,238,
  45. 239,240,240,241,241,242,242,243,
  46. 243,244,244,245,245,246,246,247,
  47. 247,248,248,249,249,250,250,251,
  48. 251,252,252,253,253,254,254,255
  49. };
  51. extern "C"
  52. DWORD __fastcall UintSqrt(DWORD n)
  53. {
  54. DWORD bc;
  55. DWORD r;
  56. DWORD r1;
  58. if (n<64)
  59. return sqrtTab[n];
  61. bc=log2(n)/2;
  62. switch (bc)
  63. {
  64. case 3:
  65. r=(sqrtTab[n]>>4);
  66. r1=r+1;
  67. if (r1*r1<=n)
  68. r++;
  69. break;
  70. case 4:
  71. r=(sqrtTab[n>>2]>>3);
  72. r1=r+1;
  73. if (r1*r1<=n)
  74. r++;
  75. break;
  77. case 5:
  78. r=(sqrtTab[n>>4]>>2);
  79. r1=r+1;
  80. if (r1*r1<=n)
  81. r++;
  82. break;
  84. case 6:
  85. r=(sqrtTab[n>>6]>>1);
  86. r1=r+1;
  87. if (r1*r1<=n)
  88. r++;
  90. break;
  92. case 7:
  93. r=(sqrtTab[n>>8]);
  94. r1=r+1;
  95. if (r1*r1<=n)
  96. r++;
  97. break;
  99. case 8:
  100. r=(sqrtTab[n>>10]<<1);
  101. r= (n/r + r)/2;
  102. if (r*r>n) r--;
  103. break;
  105. case 9:
  106. r=(sqrtTab[n>>12]<<2);
  107. r= (n/r + r)/2;
  108. if (r*r>n) r--;
  109. break;
  111. case 10:
  112. r=(sqrtTab[n>>14]<<3);
  113. r= (n/r + r)/2;
  114. if (r*r>n) r--;
  115. break;
  117. case 11:
  118. r=(sqrtTab[n>>16]<<4);
  119. r= (n/r + r)/2;
  120. if (r*r>n) r--;
  121. break;
  123. case 12:
  124. r=(sqrtTab[n>>18]<<5);
  125. r= (n/r + r)/2;
  126. if (r*r>n) r--;
  127. break;
  129. case 13:
  130. r=(sqrtTab[n>>20]<<6);
  131. r= (n/r + r)/2;
  132. if (r*r>n) r--;
  133. break;
  135. case 14:
  136. r=(sqrtTab[n>>22]<<7);
  137. r= (n/r + r)/2; //2次牛顿迭代
  138. r= (n/r + r)/2;
  139. if (r*r>n) r--;//调整
  140. break;
  142. case 15:
  143. r=(sqrtTab[n>>24]<<8);
  144. r= (n/r + r)/2; //2次牛顿迭代
  145. r= (n/r + r)/2;
  146. if (r*r>n) r--;//调整
  147. break;
  148. }
  149. return r;
  150. }
复制代码
[ 本帖最后由 liangbch 于 2009-2-13 10:45 编辑 ]
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

无心人
 楼主| 发表于 2009-2-10 21:43:52 | 显示全部楼层
我说的就是这个意思 我可能有时间测试下256表的一次迭代的范围 毕竟256表的范围分解很好做 ======================================= 不行， 在测试16位整数时就有一次迭代误差太大的例子 ======================================== 又及：否定上面的否定 换个思路 明天重新测试
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

无心人
 楼主| 发表于 2009-2-11 10:54:39 | 显示全部楼层
刚才测试FPU及SSE2的整数平方根程序 发现仅能计算有符号数的 无符号数的，比如2^31 得到的结果是错误的 需要修正
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

liangbch
发表于 2009-2-11 11:03:34 | 显示全部楼层
如果扩大查找表，应该可以减少牛顿迭代次数。一个方案是： 1. BYTE 数组。256个元素， 直接得到 0到255的平方根。 2. WORD数组，768个元素，可得到256 to 1023的平方根，精确到16bit。 总共需要256+768*2=1792 Bytes, 在计算2^32以内的平方根时，至多使用一次牛顿迭代。
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

liangbch
发表于 2009-2-11 11:05:49 | 显示全部楼层

回复 11# liangbch 的帖子

11# 中case 3 to case 7中的
if ( n-r*r>r*2) r++;
可改为 if ( (r+1)*(r+1)
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

仙剑魔
发表于 2009-2-11 11:11:35 | 显示全部楼层
据说有个"利用1+3+5+7+...+(2n-1)=n^2 原理,只用到加法和移位"的算法 虽然我还没学会...
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

无心人
 楼主| 发表于 2009-2-11 11:51:21 | 显示全部楼层
楼上的主意很慢的
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

无心人
 楼主| 发表于 2009-2-11 14:09:42 | 显示全部楼层
  1. double b32 = 4294967296.0;
  3. __declspec(naked)
  4. DWORD __fastcall iSqrt_SSE3(DWORD n)
  5. {
  6. __asm
  7. {
  8. push ecx
  9. fld qword ptr [b32]
  10. fldz
  11. cmp ecx, 0x8000000
  12. fcmovnb st, st(1)
  13. ffree st(1)
  14. fild dword ptr [esp]
  15. faddp st(1), st
  16. fsqrt
  17. fisttp dword ptr [esp]
  18. pop eax
  19. ret
  20. }
  21. }
复制代码
SSE3版，已进行修正 可正常处理超过2^31-1的整数！！ ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ 另外，似乎及时释放寄存器 会极大提高速度 见ffree一句 ========================= 可惜了，似乎这个版本，大部分人机器不能运行 还是要重新写个X86的FPU版本
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

无心人
 楼主| 发表于 2009-2-11 14:24:17 | 显示全部楼层
  1. double b32 = 4294967296.0;
  3. __declspec(naked)
  4. DWORD __fastcall iSqrt_SSE2(DWORD n)
  5. {
  6. __asm
  7. {
  8. cvtsi2sd xmm0, ecx
  9. cmp ecx, 0x80000000
  10. jb FixOver
  11. movsd xmm1, qword ptr [b32]
  12. addsd xmm0, xmm1
  13. FixOver:
  14. sqrtsd xmm0, xmm0
  15. cvttsd2si eax, xmm0
  16. ret
  17. }
  18. }
复制代码
SSE2版本，　同样修正符号整数问题
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

无心人
 楼主| 发表于 2009-2-11 17:00:55 | 显示全部楼层
循环，从1到2^32-1 时间比较： SSE3 Version: 155727.475 ms SSE2 Version: 172307.976 ms
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

12345678910... 15下一页
返回列表 发新帖
高级模式
B Color Image Link Quote Code Smilies
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2024-11-21 21:28 , Processed in 0.030555 second(s), 15 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表