一个快速浮点开方程序,据说能比float sqrt(x)快4倍
据说能比float sqrt(x)快4倍，在TC2.0下编译通过，计算结果准确度高。float SquareRootFloat(float number)
{
long i;
float x, y;
const float f = 1.5F;
x = number * 0.5F;
y= number;
i= * ( long * ) &y;
i= 0x5f3759df - ( i >> 1 );//卡马克
//i= 0x5f375a86 - ( i >> 1 );//Lomont
y= * ( float * ) &i;
y= y * ( f - ( x * y * y ) );
y= y * ( f - ( x * y * y ) );
return number * y;
}
有关资料
中文
http://blog.donews.com/snailact/archive/2006/04/01/806368.aspx
英文
http://www.codemaestro.com/reviews/review00000105.html
http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf
-----此内容被123_zh于2006-06-26,23:34:30编辑过 太NB了。按也作过牛顿迭代开方，不过一般是10次。 太谢谢楼主了。我找到的资料都是1/sqrt(x)的功能，就是求出的结果是平方根的倒数。
把return y; 换为return number * y; 就和楼主的一样，求出的是平方根了。
我的测试环境：ICC AVR 7.20 ＋AVR Studio 4.14 + M16。
编译后占(8M的flash)大小被开方数 结果 占用时钟周期数（Cycle counter）
自带的函数sqrtf(x); 15% 65536 256 7981
SquareRootFloat(float number): 8% 65536 255.9989 3475 刚才又在网上找到了一篇：就是楼主给出的函数！
[直译]Quake III中不可思议的求解平方根实现方法
Quake III中不可思议的求解平方根实现方法
任何一个3D引擎都是通过其内部的数学模型和实现工具来展现它的力量与速度的，and trust John Carmack of ID software for using really good hacks. 结果，Quake III中使用了一个非常有意思的技巧来计算平方根倒数（inverse square root）
前言
ID software最近发布了它的带有Gpl许可证的Quake III引擎源代码，在这篇文章中我们将会看到Carmark是怎样用他的black magic来极其迅速地计算一个浮点数的平方根的。
Carmack's 不寻常平方根倒数
对文件game/code/q_math.c的快速一瞥就显示出了许多有趣的performance hacks。
第一个跳出来的便是对函数Q_rsqrt中对0x5f3759df的使用，这个数计算了一个浮点数的inverse square root，但是为什么这个函数有这样的功能呢？
观察q_math.c原本的函数：
float Q_rsqrt( float number )
{
long i;
float x2, y;
const float threehalfs = 1.5F;
x2 = number * 0.5F;
y= number;
i= * ( long * ) &y;// evil floating point bit level hacking
i= 0x5f3759df - ( i >> 1 ); // what the (敏感词0386)?
y= * ( float * ) &i;
y= y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
// y= y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed
#ifndef Q3_VM
#ifdef __linux__
assert( !isnan(y) ); // bk010122 - FPE?
#endif
#endif
return y;
}
br/>它不仅有效，甚至在某些CPU上，Carmack的Q_rsqrt 比(float)(1.0/sqrt(x)的计算快4倍，尽管sqrt()通常使用的是FSQRT的汇编指令！
在另一个文件code/common/cm_trace.c 中，我们发现了更简洁的对同样HACK的实现。这一次，它被用来计算一个float - sqrt(x)的平方根。注意，其中的唯一不同是在返回值上－－用返回*y取代了返回y。
/*
================
SquareRootFloat
================
*/
float SquareRootFloat(float number) {
long i;
float x, y;
const float f = 1.5F;
x = number * 0.5F;
y= number;
i= * ( long * ) &y;
i= 0x5f3759df - ( i >> 1 );
y= * ( float * ) &i;
y= y * ( f - ( x * y * y ) );
y= y * ( f - ( x * y * y ) );
return number * y;
}
牛顿对根的近似值
上面的代码执行了众所周知的牛顿对根的近似值,像绝大多数其它迭代求近似值的计算一样，牛顿近似值假定是迭代的；每一次迭代都增强了它的准确度直至达到需要的准确度。
在牛顿近似值中的一般想法是我们我们猜测一个数x的平方根值y，我们可能通过一个简单的操作用x/y来拉平y来取得更好的猜测，使其更接近实际的平方根，例如，我们像下面这样计算2的平方根，我们假定初始的猜测是1：
2/1 = 2 ;(2 + 1) / 2 = 1.5
2/1.5 = 1.3333; ( 1.5 + 1.3333 ) / 2 = 1.4167
2/1.4167 = 1.4117;( 1.4167 + 1.4117 ) / 2 = 1.4142
And so on...
如前面所提到的，牛顿的近似值是一个大家所熟知的用以快速计算平方根的方法。但是，Carmack在初始的猜测中就选取的不寻常的值，它彻底加强了准确度并且将Quake III中计算所要的值的迭代次数降到了1次！
魔数
这个函数中真正有意思的方面是神奇的常量0x5f3759df，用来计算初始猜测的，在
i= 0x5f3759df - ( i >> 1 );
因此，把输入除以2并从神奇常量中减去。这个常数工作起来几乎是完美的－－对于一个 low relative error of 10^-3来说只要一次牛顿近似值迭代就够了。如评论中第二次迭代中展示的，这个近似值对Quake III引擎来说已经足够了。
结果，这个神奇的常数0x5f3759df成了一个迷了，在文章"Fast Inverse Square Root" ，普度大学的数学家Chris Lomont研究了这个常数，用了几种精细的技术，Lomont想自己用数学方法求出这个常数来，结果令人惊奇－－Lomont用数学方法计算出来的最佳常数(0x5f37642f)有一点点不同，并且除了理论上强一些之外，它产生的结果并没有源代码中使用的原始常数好！确实，John Carmack 一定用了天才般的黑盒来找到这个常数。
只在仅仅从数字上来找的方法中，Lomont找到了一个更好的常数，这个数比原始的那个强了那么一点点。然而，实践中两个常数产生了大概相同的结果，Lomont提出这个使用了更好的常数的函数：
float InvSqrt(float x)
{
float xhalf = 0.5f*x;
int i = *(int*)&x; // get bits for floating value
i = 0x5f375a86- (i>>1); // gives initial guess y0
x = *(float*)&i; // convert bits back to float
x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
return x;
}
Trackback: http://tb.donews.net/TrackBack.aspx?PostId=806368 不可思议的神奇 顶 强悍！ mark mark mark MARK 好多学问 mark 数学的力量…… mark 按个爪印
哪位高人验算一下
用这个算法开方再乘回去看看误差多少 mark 顶 谢谢LZ分享，马克！ 记号 好！ 顶 谢谢 mark 牛人呀 数学是非常有用的，顶！！！ 学习了 it's amazing you rock!
make mark. 神秘 先记号 挖到宝!!! 好东西，留名，以后一定会用得着。 不知道是否我的计算机有问题，怎么SquareRootFloat()比AVR STUDIO的sqrt()慢好多！ mark 留名 mark 强贴留印 mark sqrt mark 这个是否落后了？
我用汇编写的可以在GCC下直接使用的fsqrt_qianhng()执行才526个机器周期，是GCC float sqrt()的1/8。嘿嘿 额呵呵。年轻的朋友都喜欢呀。 mark 留个脚印... 顶一下 顶一个 标记 mark mark 好帖，mark mark mark 记号 winavr 的sqrt 大概494个时钟周期 回复【53楼】jingle jingle
winavr 的sqrt 大概494个时钟周期
-----------------------------------------------------------------------
不会吧，哪个版本的winavr,是浮点数的吗？你开方根3.1试试 在keil下仿真:
//STM32@72M
float InvSqrt(float x)
{
float xhalf = 0.5f*x;
int i = *(int*)&x; // get bits for floating value
i = 0x5f375a86- (i>>1); // gives initial guess y0
x = *(float*)&i; // convert bits back to float
x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
return x;
}
int main(void)
{
u8 temp;
float testnum=3.1;
float res,res1;
Stm32_Clock_Init(9);//系统时钟设置 12M*6=72M
delay_init(72); //延时初始化
uart_init(72,57600);//串口1初始化
res=InvSqrt(testnum); //时间:0.00095323
res1=1/sqrt(testnum); //时间:0.00095551
printf("res:%f",res); //时间:0.00096559
printf("res1:%f",res1);
得到的结果:res:0.567654
res1:0.5679619
windows计算结果:0.5679618
可以看出InvSqrt的计算只花了2.28us,而第二个res1使用sqrt函数,使用时间为10.08us
后者是前者的4.42倍!!
有时候还是很有用的,这个函数.
只能感叹数学真的很神奇!
接着楼主的意思,楼主得到的是1/sqrt,我们很多时候只要SQRT就够了,并不需要倒数
所以网上又搜到以快速sqrt的函数:
float CarmSqrt(float x)
{
union
{
int intPart;
float floatPart;
} convertor;
union
{
int intPart;
float floatPart;
} convertor2;
convertor.floatPart = x;
convertor2.floatPart = x;
convertor.intPart = 0x1FBCF800 + (convertor.intPart >> 1);
convertor2.intPart = 0x5f3759df - (convertor2.intPart >> 1);
return 0.5f*(convertor.floatPart + (x * convertor2.floatPart));
}
再次来验证速度:
float InvSqrt(float x)
{
float xhalf = 0.5f*x;
int i = *(int*)&x; // get bits for floating value
i = 0x5f3759df - (i>>1); // gives initial guess y0
x = *(float*)&i; // convert bits back to float
x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
return x;
}
float CarmSqrt(float x)
{
union
{
int intPart;
float floatPart;
} convertor;
union
{
int intPart;
float floatPart;
} convertor2;
convertor.floatPart = x;
convertor2.floatPart = x;
convertor.intPart = 0x1FBCF800 + (convertor.intPart >> 1);
convertor2.intPart = 0x5f3759df - (convertor2.intPart >> 1);
return 0.5f*(convertor.floatPart + (x * convertor2.floatPart));
}
//STM32@72M
int main(void)
{
u8 temp;
float testnum=3.1;
float res,res1,res2;
Stm32_Clock_Init(9);//系统时钟设置 12M*6=72M
delay_init(72); //延时初始化
uart_init(72,57600);//串口1初始化
res=CarmSqrt(testnum); //0.00095335
res2=1/InvSqrt(testnum); //0.00095492
res1=sqrt(testnum); //0.00095791
printf("res:%f",res); //0.00096533
printf("res2:%f",res2);
printf("res1:%f",res1);
通过上面比较,CarmSqrt函数计算时间为1.57us,1/InvSqrt计算时间为:2.99us,系统sqrt函数的计算时间为:7.42us.
看出来CarmSqrt函数的性能在开平方的时候,最好了.
以后,可以试试这个,精度在10的-3次方内,这个函数可以节约很多时间. 附上e文,给牛人研究研究.说不定以后有更好的用.^_^
点击此处下载 ourdev_534129.pdf(文件大小:148K) (原文件名:InvSqrt.pdf) 记号了 不错！！！！！ 快速浮点开方程序，我是GOOGLE搜过来的找到的 高中时老师教过一个徒手开跟的方法，只要用到乘法和加减法，不知是什么原理。 http://cache.amobbs.com/bbs_upload782111/files_28/ourdev_548808.JPG
(原文件名:徒手开跟.JPG) 顶一顶. 回复【楼主位】123_zh 多来米
-----------------------------------------------------------------------
dddddddddd to：【61楼】 wshtyr
顺便问一下有没有徒手算对数的方法？不用级数展开的。 标记。 回复【60楼】wshtyr
高中时老师教过一个徒手开跟的方法，只要用到乘法和加减法，不知是什么原理。
-----------------------------------------------------------------------
回复【61楼】wshtyr
http://cache.amobbs.com/bbs_upload782111/files_28/ourdev_548808.JPG
(原文件名:徒手开跟.JPG)
-----------------------------------------------------------------------
原理就是(a*10+b)^2=(a^2)*100+20*a*b+b^2 mark 数论这东西，果然够神奇。 mark 数学中的美
尽情体会吧 mark mark 很好的东西 mark mark 推荐大家用cordic坐标旋转算法 标记 这个当然要马克啊。 MARK 这个当然要马克啊。YE. mark 收藏了
非常好！ 收藏 收藏 mark 数论忘得差不多了。。。 数学没学好，是个悲剧.. 记号 mark,感觉数学就一个字“变”！ 挺好 mark 简直神奇呢！ 整型数的适合么？ mark 注意误差接近2%,某些情况不能接受 真神奇
！ mark 谢谢 mark mark 回复【1楼】liuqian刘汧
太nb了。按也作过牛顿迭代开方，不过一般是10次。
-----------------------------------------------------------------------
我做过开4次方的，但迭代4次的精度已经很高了，10次感觉没必要。
页:
[1]
2