二进制32位整数快速平方根 - 第2页 - 编程擂台 - 数学研发论坛

liangbch 发表于 2009-2-10 21:20:10

下面给出实现这个功能的代码和测试代码。
顺便说一下，在2楼长提到的各个算法的选择都是在实验的基础上得出的，既能保证正确性，又是的速度不致太慢。#include <math.h>

typedef unsigned long DWORD;
typedef unsigned char BYTE;

inline DWORD log2(DWORD n)
{
_asm
{
mov ecx,n
bsr eax,ecx
}
}

extern "C" unsigned char sqrtTab[]=
{
      0,1,1,1,2,2,2,2,
      2,3,3,3,3,3,3,3,
      4,4,4,4,4,4,4,4,
      4,5,5,5,5,5,5,5,
      5,5,5,5,6,6,6,6,
      6,6,6,6,6,6,6,6,
      6,7,7,7,7,7,7,7,
      7,7,7,7,7,7,7,7,
   128,128,129,130,131,132,133,134,
   135,136,137,138,139,140,141,142,
   143,144,144,145,146,147,148,149,
   150,150,151,152,153,154,155,155,
   156,157,158,159,160,160,161,162,
   163,163,164,165,166,167,167,168,
   169,170,170,171,172,173,173,174,
   175,176,176,177,178,178,179,180,
   181,181,182,183,183,184,185,185,
   186,187,187,188,189,189,190,191,
   192,192,193,193,194,195,195,196,
   197,197,198,199,199,200,201,201,
   202,203,203,204,204,205,206,206,
   207,208,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,224,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,240,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
};

extern "C"
DWORD __fastcall UintSqrt(DWORD n)
{
DWORD bc;
DWORD r;
DWORD r1;

if (n<64)
return sqrtTab;

bc=log2(n)/2;
switch (bc)
{
case3:
r=(sqrtTab>>4);
r1=r+1;
if (r1*r1<=n)
r++;
break;
case4:
r=(sqrtTab>>3);
r1=r+1;
if (r1*r1<=n)
r++;
break;

case 5:
r=(sqrtTab>>2);
r1=r+1;
if (r1*r1<=n)
r++;
break;

case 6:
r=(sqrtTab>>1);
r1=r+1;
if (r1*r1<=n)
r++;

break;

case 7:
r=(sqrtTab);
r1=r+1;
if (r1*r1<=n)
r++;
break;

case 8:
r=(sqrtTab<<1);
r= (n/r + r)/2;
if (r*r>n) r--;
break;

case 9:
r=(sqrtTab<<2);
r= (n/r + r)/2;
if (r*r>n) r--;
break;

case 10:
r=(sqrtTab<<3);
r= (n/r + r)/2;
if (r*r>n) r--;
break;

case 11:
r=(sqrtTab<<4);
r= (n/r + r)/2;
if (r*r>n) r--;
break;

case 12:
r=(sqrtTab<<5);
r= (n/r + r)/2;
if (r*r>n) r--;
break;

case 13:
r=(sqrtTab<<6);
r= (n/r + r)/2;
if (r*r>n) r--;
break;

case 14:
r=(sqrtTab<<7);
r= (n/r + r)/2;//2次牛顿迭代
r= (n/r + r)/2;
if (r*r>n) r--;//调整
break;

case 15:
r=(sqrtTab<<8);
r= (n/r + r)/2;//2次牛顿迭代
r= (n/r + r)/2;
if (r*r>n) r--;//调整
break;
}
return r;
}

[ 本帖最后由 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)<n) r++;

仙剑魔 发表于 2009-2-11 11:11:35

据说有个"利用1+3+5+7+...+(2n-1)=n^2 原理,只用到加法和移位"的算法
虽然我还没学会...:L

无心人 发表于 2009-2-11 11:51:21

:)

楼上的主意很慢的

无心人 发表于 2009-2-11 14:09:42

double b32 = 4294967296.0;

__declspec(naked)
DWORD __fastcall iSqrt_SSE3(DWORD n)
{
__asm
{
push ecx
fld qword ptr
fldz
cmp ecx, 0x8000000
fcmovnb st, st(1)
ffree st(1)
fild dword ptr
faddp st(1), st
fsqrt
fisttp dword ptr
pop eax
ret
}
}SSE3版，已进行修正
可正常处理超过2^31-1的整数！！
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
另外，似乎及时释放寄存器
会极大提高速度
见ffree一句
=========================
可惜了，似乎这个版本，大部分人机器不能运行
还是要重新写个X86的FPU版本

无心人 发表于 2009-2-11 14:24:17

double b32 = 4294967296.0;

__declspec(naked)
DWORD __fastcall iSqrt_SSE2(DWORD n)
{
__asm
{
cvtsi2sd xmm0, ecx
cmp ecx, 0x80000000
jb FixOver
movsd xmm1, qword ptr
addsd xmm0, xmm1
FixOver:
sqrtsd xmm0, xmm0
cvttsd2si eax, xmm0
ret
}
}
SSE2版本，　同样修正符号整数问题

无心人 发表于 2009-2-11 17:00:55

循环，从1到2^32-1
时间比较：
SSE3 Version: 155727.475 ms
SSE2 Version: 172307.976 ms

页: 1 [2] 3 4 5 6 7 8 9 10 11

数学研发论坛's Archiver

回复 11# liangbch 的帖子