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

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

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

[复制链接]
无心人
 楼主| 发表于 2009-2-12 15:55:41 | 显示全部楼层
你39#代码存在问题
首行应该是or ecx, ecx吧
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

liangbch
发表于 2009-2-12 16:03:46 | 显示全部楼层

回复 41# 无心人 的帖子

是，不小心写错了，已修改

下面贴出在无心人的基础上，经我改写的3个FPU的版本。


  3. double b32[] = {0.0,4294967296.0};
  4. //实验证明zero5不能大于0.49999999999636
  5. double zero5=        0.49999999999636;;

  7. __declspec(naked)
  8. DWORD __fastcall iSqrt_FPU1(DWORD n)
  9. {
  10. __asm
  11. {
  12. push ecx

  14. shr ecx, 31
  15. fld qword ptr [b32 + ecx  * 8]
  16. fild dword ptr [esp]
  17. faddp  st(1),st
  18. fsqrt
  19. fsub  qword ptr [zero5]
  20. fistp  dword ptr [esp]
  21. pop eax
  22. ret
  23. }
  24. }

  26. __declspec(naked)
  27. DWORD __fastcall iSqrt_FPU2(DWORD n)
  28. {
  29. __asm
  30. {
  31. or ecx,ecx
  32. jnz next00
  33. mov eax,0
  34. ret
  35. next00:
  36. push ecx

  38. shr ecx, 31
  39. fld qword ptr [b32 + ecx * 8]
  40. fild dword ptr [esp]
  41. faddp  st(1),st
  42. fsqrt
  43. fstp  qword ptr [esp-8]

  45. fwait
  46. mov  ecx,dword ptr [esp-4]
  47. mov  edx,0xfff00000
  48. mov  eax,0xfffff

  50. and  edx,ecx //阶码
  51. and  eax,ecx //尾数
  52. shr  edx,20 //得到阶段码
  53. add  eax,0x100000 //设置位数最高有效位
  54. mov  ecx,1043
  55. sub  ecx,edx
  56. shr  eax,cl //得到整数

  58. pop ecx
  59. ret
  60. }
  61. }

  63. __declspec(naked)
  64. DWORD __fastcall iSqrt_FPU3(DWORD n)
  65. {
  66. __asm
  67. {
  68. push ecx
  69. shr ecx, 31
  70. fld qword ptr [b32 + ecx * 8]
  71. fild dword ptr [esp]
  72. faddp  st(1),st
  73. fsqrt

  75. //设置FPU的取整方式  为了直接使用fistp浮点指令
  76. FNSTCW  [esp-4]            // 保存协处理器控制字,用来恢复
  77. FNSTCW  [esp-6]            // 保存协处理器控制字,用来修改
  78. FWAIT
  79. OR  word ptr [esp-6], 0x0F00    // 改为 RC=11  使FPU向零取整
  80. FLDCW   [esp-6]            // 载入协处理器控制字,RC场已经修改

  82. fistp   dword ptr [esp]

  84. //恢复FPU的取整方式
  85. FWAIT
  86. FLDCW   [esp-4]

  88. pop eax
  89. ret
  90. }
  91. }

复制代码

[ 本帖最后由 liangbch 于 2009-2-13 10:28 编辑 ]
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

gxqcn
发表于 2009-2-12 16:26:35 | 显示全部楼层
我也写了个FPU的版本，有兴趣的朋友可把相关函数添加进来一并测试一下：
  1. #include < stdio.h >
  2. #include < stdlib.h >
  3. #include < time.h >
  4. #include < float.h >

  6. typedef unsigned int UINT32;
  7. typedef UINT32 __fastcall func_t( UINT32 );

  9. _declspec(naked)
  10. UINT32 __fastcall Sqrt_43F( UINT32 n )
  11. {
  12.     __asm
  13.     {
  14.         xor     eax, eax;
  15.         push    eax;
  16.         push    ecx;

  18.         fild    qword ptr[esp];
  19.         fsqrt;
  20.         pop     eax;

  22.         fistp   dword ptr[esp];
  23.         pop     eax;

  25.         ret;
  26.     }
  27. }

  29. #define FN(ID)      Sqrt_##ID##F /* ex.: Sqrt_43F */

  31. void test_fun( void )
  32. {
  33.     UINT32 i, j;
  34.     UINT32 fID[] = { 43, /* ... */ };
  35.     func_t * funcTable[] = { FN(43), /* FN(..),... */ };
  36.     func_t * pfn;

  38. #ifdef _DEBUG
  39.     UINT32 m = 256;
  40. #else
  41.     UINT32 m = 1UL<<28;
  42.     time_t t;
  43. #endif

  45.     printf( "\n\nn: 0-0x%x\n", m );
  46.     for ( i=0; i<sizeof(fID)/sizeof(fID[0]); ++i )
  47.     {
  48.         pfn = funcTable[i];
  49. #ifdef _DEBUG
  50.         // check result
  51.         for ( j=0; j<=m; ++j )
  52.         {
  53.             printf( "Sqrt_%uF(%u)=%u\n", fID[i], j, pfn(j) );
  54.         }
  55.         system( "pause" );
  56. #else
  57.         t = time(NULL);
  58.         for ( j=m; 0!=j; --j )
  59.         {
  60.             pfn(j);
  61.         }
  62.         t = time(NULL)-t;

  64.         printf( "%u#: %u s\n", fID[i], t );
  65. #endif
  66.     }
  67. }

  69. int main(int argc, char* argv[])
  70. {
  71.     _control87( RC_CHOP | PC_64, MCW_RC | MCW_PC );
  72.     test_fun();

  74.     return 0;
  75. }
复制代码
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

无心人
 楼主| 发表于 2009-2-12 16:42:16 | 显示全部楼层
宝宝：
  减0.5不是个好办法，或许因为别的函数改变了RC， 就反而得不到正确的数值了

锅老大：
你这算偷懒的
在函数外
控制RC获得的性能提升有限的
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

liangbch
发表于 2009-2-12 16:44:55 | 显示全部楼层
gxq 的这个版本最简单，但是在程序开头改变了FPU控制字，不知到在运行期间，是不是会影响的一些浮点数学函数。
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

无心人
 楼主| 发表于 2009-2-12 17:07:17 | 显示全部楼层

  2. double b32[] = {0.0,4294967296.0};

  4. __declspec(naked)
  5. DWORD __fastcall iSqrt_FPU(DWORD n)
  6. {
  7. __asm
  8. {
  9.    push ecx
  10.    sub esp, 4
  11.    fnstcw word ptr  [esp]
  12.    mov edx, dword ptr [esp]
  13.    or dword ptr [esp], 0x0C00
  14.    mov eax, ecx   
  15.    and eax, 0x80000000
  16.    shr eax, 31
  17.    fld qword ptr [b32 + eax * 8]
  18.    fild dword ptr [esp + 4]
  19.    faddp st(1), st
  20.    fsqrt
  21.    fldcw word ptr [esp]
  22.    fistp dword ptr [esp + 4]
  23.    mov dword ptr [esp], edx
  24.    fldcw word ptr [esp]
  25.    add esp, 4
  26.    pop eax
  27.    ret
  28.    }
  29. }

  31. __declspec(naked)
  32. DWORD iSqrt_FPU1(DWORD n)
  33. {
  34. __asm
  35. {
  36.    sub esp, 4
  37.    fnstcw word ptr  [esp]
  38.    mov edx, dword ptr [esp]
  39.    or dword ptr [esp], 0x0C00
  40.    mov eax, [esp + 8]   
  41.    and eax, 0x80000000
  42.    shr eax, 31
  43.    fld qword ptr [b32 + eax * 8]
  44.    fild dword ptr [esp + 8]
  45.    faddp st(1), st
  46.    fsqrt
  47.    fldcw word ptr [esp]
  48.    fistp dword ptr [esp + 8]
  49.    mov dword ptr [esp], edx
  50.    fldcw word ptr [esp]
  51.    add esp, 4
  52.    mov eax, [esp + 4]
  53.    ret
  54.    }
  55. }

  57. __declspec(naked)
  58. DWORD __fastcall iSqrt_FPU2(DWORD n)
  59. {
  60. __asm
  61. {
  62.    push ecx
  63.    mov eax, ecx   
  64.    and eax, 0x80000000
  65.    shr eax, 31
  66.    fld qword ptr [b32 + eax * 8]
  67.    fild dword ptr [esp]
  68.    faddp st(1), st
  69.    fsqrt
  70.    sub esp, 8
  71.    fstp qword ptr [esp]
  72.    mov edx, dword ptr [esp + 4]
  73.    mov eax, edx
  74.    and edx,0x7ff00000
  75.    and eax,0xfffff
  76.    shr edx, 20
  77.    or eax, 0x100000
  78.    xchg ecx, edx
  79.    sub ecx, 1043
  80.    neg ecx
  81.    shr eax, cl
  82.    xchg edx, ecx
  83.    add esp, 12
  84.    or ecx, ecx
  85.    cmove eax, ecx
  86.    ret
  87.    }
  88. }

复制代码
比较时间
FPU2 < FPU < FPU1

另外，程序开始就改变RC的
会造成函数结果错误
输入小的时候测试不到
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

无心人
 楼主| 发表于 2009-2-12 17:19:22 | 显示全部楼层
下面就是寻求执行时间少于60周期的整数算法
考虑乘法总延迟16周期
除法总延迟是50周期

所以
1、不能出现除法
2、乘法最多2个
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

liangbch
发表于 2009-2-12 17:26:29 | 显示全部楼层
你的FPU2的版本写成
sub edx, 1022
   mov ecx, 21
   sub ecx, edx


我的则做了优化，直接写成以下形式，实际上是等价的
  mov  ecx,1043
  sub  ecx,edx
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

无心人
 楼主| 发表于 2009-2-12 17:28:50 | 显示全部楼层


我不理解你的做法
所以改了啊
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

无心人
 楼主| 发表于 2009-2-12 17:52:24 | 显示全部楼层
呵呵

    在持续改进FPU2函数，增加了非跳转处理0的代码

   另外，FWAIT似乎并不必须吧
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2024-4-26 18:41 , Processed in 0.100902 second(s), 15 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

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