找回密码
 欢迎注册
楼主: 无心人

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

[复制链接]
发表于 2009-2-10 21:20:10 | 显示全部楼层
下面给出实现这个功能的代码和测试代码。
顺便说一下,在2楼长提到的各个算法的选择都是在实验的基础上得出的,既能保证正确性,又是的速度不致太慢。
  1. #include <math.h>

  2. typedef unsigned long DWORD;
  3. typedef unsigned char BYTE;

  4. inline DWORD log2(DWORD n)
  5. {
  6.         _asm
  7.         {
  8.                 mov ecx,n
  9.                 bsr eax,ecx
  10.         }
  11. }

  12. extern "C" unsigned char sqrtTab[]=
  13. {
  14.           0,  1,  1,  1,  2,  2,  2,  2,
  15.           2,  3,  3,  3,  3,  3,  3,  3,
  16.           4,  4,  4,  4,  4,  4,  4,  4,
  17.           4,  5,  5,  5,  5,  5,  5,  5,
  18.           5,  5,  5,  5,  6,  6,  6,  6,
  19.           6,  6,  6,  6,  6,  6,  6,  6,
  20.           6,  7,  7,  7,  7,  7,  7,  7,
  21.           7,  7,  7,  7,  7,  7,  7,  7,
  22.         128,128,129,130,131,132,133,134,
  23.         135,136,137,138,139,140,141,142,
  24.         143,144,144,145,146,147,148,149,
  25.         150,150,151,152,153,154,155,155,
  26.         156,157,158,159,160,160,161,162,
  27.         163,163,164,165,166,167,167,168,
  28.         169,170,170,171,172,173,173,174,
  29.         175,176,176,177,178,178,179,180,
  30.         181,181,182,183,183,184,185,185,
  31.         186,187,187,188,189,189,190,191,
  32.         192,192,193,193,194,195,195,196,
  33.         197,197,198,199,199,200,201,201,
  34.         202,203,203,204,204,205,206,206,
  35.         207,208,208,209,209,210,211,211,
  36.         212,212,213,214,214,215,215,216,
  37.         217,217,218,218,219,219,220,221,
  38.         221,222,222,223,224,224,225,225,
  39.         226,226,227,227,228,229,229,230,
  40.         230,231,231,232,232,233,234,234,
  41.         235,235,236,236,237,237,238,238,
  42.         239,240,240,241,241,242,242,243,
  43.         243,244,244,245,245,246,246,247,
  44.         247,248,248,249,249,250,250,251,
  45.         251,252,252,253,253,254,254,255
  46. };

  47. extern "C"
  48. DWORD __fastcall UintSqrt(DWORD n)
  49. {
  50.         DWORD bc;
  51.         DWORD r;
  52.         DWORD r1;
  53.        
  54.         if (n<64)
  55.                 return sqrtTab[n];

  56.         bc=log2(n)/2;
  57.         switch (bc)
  58.         {
  59.         case  3:
  60.                 r=(sqrtTab[n]>>4);
  61.                 r1=r+1;
  62.                 if (r1*r1<=n)
  63.                         r++;
  64.                 break;
  65.         case  4:
  66.                 r=(sqrtTab[n>>2]>>3);
  67.                 r1=r+1;
  68.                 if (r1*r1<=n)
  69.                         r++;
  70.                 break;

  71.         case 5:
  72.                 r=(sqrtTab[n>>4]>>2);
  73.                 r1=r+1;
  74.                 if (r1*r1<=n)
  75.                         r++;
  76.                 break;

  77.         case 6:
  78.                 r=(sqrtTab[n>>6]>>1);
  79.                 r1=r+1;
  80.                 if (r1*r1<=n)
  81.                         r++;
  82.                
  83.                 break;

  84.         case 7:
  85.                 r=(sqrtTab[n>>8]);
  86.                 r1=r+1;
  87.                 if (r1*r1<=n)
  88.                         r++;
  89.                 break;

  90.         case 8:
  91.                 r=(sqrtTab[n>>10]<<1);
  92.                 r= (n/r + r)/2;
  93.                 if (r*r>n) r--;
  94.                 break;

  95.         case 9:
  96.                 r=(sqrtTab[n>>12]<<2);
  97.                 r= (n/r + r)/2;
  98.                 if (r*r>n) r--;
  99.                 break;

  100.         case 10:
  101.                 r=(sqrtTab[n>>14]<<3);
  102.                 r= (n/r + r)/2;
  103.                 if (r*r>n) r--;
  104.                 break;

  105.         case 11:
  106.                 r=(sqrtTab[n>>16]<<4);
  107.                 r= (n/r + r)/2;
  108.                 if (r*r>n) r--;
  109.                 break;

  110.         case 12:
  111.                 r=(sqrtTab[n>>18]<<5);
  112.                 r= (n/r + r)/2;
  113.                 if (r*r>n) r--;
  114.                 break;

  115.         case 13:
  116.                 r=(sqrtTab[n>>20]<<6);
  117.                 r= (n/r + r)/2;
  118.                 if (r*r>n) r--;
  119.                 break;

  120.         case 14:
  121.                 r=(sqrtTab[n>>22]<<7);
  122.                 r= (n/r + r)/2;  //2次牛顿迭代
  123.                 r= (n/r + r)/2;
  124.                 if (r*r>n) r--;//调整
  125.                 break;

  126.         case 15:
  127.                 r=(sqrtTab[n>>24]<<8);
  128.                 r= (n/r + r)/2;  //2次牛顿迭代
  129.                 r= (n/r + r)/2;
  130.                 if (r*r>n) r--;//调整
  131.                 break;
  132.         }
  133.         return r;
  134. }

复制代码

[ 本帖最后由 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
得到的结果是错误的

需要修正
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 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以内的平方根时,至多使用一次牛顿迭代。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 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)<n) r++;
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 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;

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

  2. __declspec(naked)
  3. DWORD __fastcall iSqrt_SSE2(DWORD n)
  4. {
  5. __asm
  6. {
  7.     cvtsi2sd xmm0, ecx
  8. cmp ecx, 0x80000000
  9. jb FixOver
  10. movsd xmm1, qword ptr [b32]
  11. addsd xmm0, xmm1
  12. FixOver:
  13. sqrtsd xmm0, xmm0
  14. cvttsd2si eax, xmm0
  15. ret
  16. }
  17. }
复制代码
SSE2版本, 同样修正符号整数问题
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2009-2-11 17:00:55 | 显示全部楼层
循环,从1到2^32-1
时间比较:
SSE3 Version: 155727.475 ms
SSE2 Version: 172307.976 ms
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

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

GMT+8, 2024-3-29 08:08 , Processed in 0.044299 second(s), 15 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

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