搜索
bottom↓
回复: 152

一个快速浮点开方程序,据说能比float sqrt(x)快4倍

  [复制链接]

出0入0汤圆

发表于 2006-6-26 23:30:51 | 显示全部楼层 |阅读模式
据说能比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编辑过

阿莫论坛20周年了!感谢大家的支持与爱护!!

月入3000的是反美的。收入3万是亲美的。收入30万是移民美国的。收入300万是取得绿卡后回国,教唆那些3000来反美的!

出0入4汤圆

发表于 2006-6-28 14:33:35 | 显示全部楼层
太NB了。按也作过牛顿迭代开方,不过一般是10次。

出0入0汤圆

发表于 2009-8-25 13:42:02 | 显示全部楼层
太谢谢楼主了。我找到的资料都是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

出0入0汤圆

发表于 2009-8-25 13:46:05 | 显示全部楼层
刚才又在网上找到了一篇:就是楼主给出的函数!

[直译]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;

}


牛顿对根的近似值
上面的代码执行了众所周知的牛顿对根的近似值[3],像绝大多数其它迭代求近似值的计算一样,牛顿近似值假定是迭代的;每一次迭代都增强了它的准确度直至达到需要的准确度。

在牛顿近似值中的一般想法是我们我们猜测一个数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" [2] ,普度大学的数学家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

出0入42汤圆

发表于 2009-8-25 14:18:38 | 显示全部楼层
不可思议的神奇

出0入0汤圆

发表于 2009-8-25 14:24:44 | 显示全部楼层

出0入0汤圆

发表于 2009-8-25 14:42:04 | 显示全部楼层
强悍!

出0入0汤圆

发表于 2009-8-25 14:55:47 | 显示全部楼层
mark

出0入0汤圆

发表于 2009-8-25 16:45:29 | 显示全部楼层
mark mark

出0入0汤圆

发表于 2009-8-25 16:47:29 | 显示全部楼层
MARK

出0入0汤圆

发表于 2009-8-25 16:55:37 | 显示全部楼层
好多学问

出0入0汤圆

发表于 2009-8-25 18:52:59 | 显示全部楼层
mark

出0入0汤圆

发表于 2009-8-25 19:03:00 | 显示全部楼层
数学的力量……

出0入0汤圆

发表于 2009-8-25 19:08:18 | 显示全部楼层
mark

出0入0汤圆

发表于 2009-8-25 19:08:44 | 显示全部楼层
按个爪印
哪位高人验算一下
用这个算法开方再乘回去看看误差多少

出0入0汤圆

发表于 2009-8-25 19:48:38 | 显示全部楼层
mark

出0入0汤圆

发表于 2009-8-25 23:36:44 | 显示全部楼层

出0入0汤圆

发表于 2009-8-26 00:22:18 | 显示全部楼层
谢谢LZ分享,马克!

出0入0汤圆

发表于 2009-8-26 00:39:45 | 显示全部楼层
记号

出0入0汤圆

发表于 2009-8-26 00:56:49 | 显示全部楼层
好!

出0入22汤圆

发表于 2009-8-26 01:38:41 | 显示全部楼层

出0入0汤圆

发表于 2009-8-26 08:29:18 | 显示全部楼层
谢谢

出0入0汤圆

发表于 2009-8-26 08:40:53 | 显示全部楼层
mark

出0入10汤圆

发表于 2009-8-26 08:54:15 | 显示全部楼层
牛人呀

出0入0汤圆

发表于 2009-8-26 09:07:47 | 显示全部楼层
数学是非常有用的,顶!!!

出0入0汤圆

发表于 2009-8-26 09:35:11 | 显示全部楼层
学习了

出0入0汤圆

发表于 2009-8-26 09:42:39 | 显示全部楼层
it's amazing

出0入4汤圆

发表于 2009-8-26 09:55:08 | 显示全部楼层
you rock!
make mark.

出0入0汤圆

发表于 2009-8-26 09:55:24 | 显示全部楼层
神秘

出0入0汤圆

发表于 2009-8-26 09:56:34 | 显示全部楼层
先记号

出0入0汤圆

发表于 2009-11-13 11:03:44 | 显示全部楼层
挖到宝!!!

出0入0汤圆

发表于 2009-11-13 12:14:59 | 显示全部楼层
好东西,留名,以后一定会用得着。

出0入0汤圆

发表于 2009-11-13 15:24:20 | 显示全部楼层
不知道是否我的计算机有问题,怎么SquareRootFloat()比AVR STUDIO的sqrt()慢好多!

出0入0汤圆

发表于 2009-11-13 16:20:09 | 显示全部楼层
mark

出0入0汤圆

发表于 2009-11-14 11:30:17 | 显示全部楼层
留名

出0入0汤圆

发表于 2009-11-14 12:42:37 | 显示全部楼层
mark

出0入0汤圆

发表于 2009-11-14 12:56:17 | 显示全部楼层
强贴留印

出0入0汤圆

发表于 2009-11-15 13:17:08 | 显示全部楼层
mark sqrt

出0入0汤圆

发表于 2009-11-15 15:54:30 | 显示全部楼层
mark

出0入0汤圆

发表于 2009-11-15 23:03:13 | 显示全部楼层
这个是否落后了?
我用汇编写的可以在GCC下直接使用的fsqrt_qianhng()执行才526个机器周期,是GCC float sqrt()的1/8。嘿嘿

出0入0汤圆

发表于 2009-11-16 10:31:35 | 显示全部楼层
额呵呵。年轻的朋友都喜欢呀。

出0入0汤圆

发表于 2010-1-4 22:48:11 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-1-4 22:52:30 | 显示全部楼层
留个脚印...

出0入0汤圆

发表于 2010-1-5 16:53:16 | 显示全部楼层
顶一下

出0入0汤圆

发表于 2010-1-7 10:17:06 | 显示全部楼层
顶一个

出0入0汤圆

发表于 2010-1-7 10:34:39 | 显示全部楼层
标记

出0入0汤圆

发表于 2010-1-27 21:12:42 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-1-27 21:26:37 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-1-27 21:44:51 | 显示全部楼层
好帖,mark

出0入0汤圆

发表于 2010-1-27 22:51:41 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-1-27 23:06:01 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-1-28 00:18:34 | 显示全部楼层
记号

出0入0汤圆

发表于 2010-2-17 14:20:21 | 显示全部楼层
winavr 的sqrt 大概494个时钟周期

出0入0汤圆

发表于 2010-2-17 16:25:43 | 显示全部楼层
回复【53楼】jingle jingle
winavr 的sqrt 大概494个时钟周期
-----------------------------------------------------------------------

不会吧,哪个版本的winavr,是浮点数的吗?你开方根3.1试试

出0入0汤圆

发表于 2010-2-17 20:01:59 | 显示全部楼层
在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次方内,这个函数可以节约很多时间.

出0入0汤圆

发表于 2010-2-17 20:25:42 | 显示全部楼层
附上e文,给牛人研究研究.说不定以后有更好的用.^_^
点击此处下载 ourdev_534129.pdf(文件大小:148K) (原文件名:InvSqrt.pdf)

出0入0汤圆

发表于 2010-2-17 20:51:56 | 显示全部楼层
记号了

出0入0汤圆

发表于 2010-2-17 21:14:01 | 显示全部楼层
不错!!!!!

出0入0汤圆

发表于 2010-4-23 22:53:50 | 显示全部楼层
快速浮点开方程序,我是GOOGLE搜过来的找到的

出0入42汤圆

发表于 2010-4-24 00:45:10 | 显示全部楼层
高中时老师教过一个徒手开跟的方法,只要用到乘法和加减法,不知是什么原理。

出0入42汤圆

发表于 2010-4-24 00:45:33 | 显示全部楼层

(原文件名:徒手开跟.JPG)

出0入0汤圆

发表于 2010-4-24 08:27:13 | 显示全部楼层
顶一顶.

出0入0汤圆

发表于 2010-4-24 09:03:35 | 显示全部楼层
回复【楼主位】123_zh 多来米
-----------------------------------------------------------------------

dddddddddd

出0入0汤圆

发表于 2010-4-24 09:12:34 | 显示全部楼层
to:【61楼】 wshtyr

顺便问一下有没有徒手算对数的方法?不用级数展开的。

出0入0汤圆

发表于 2010-4-24 09:33:30 | 显示全部楼层
标记。

出0入0汤圆

发表于 2010-4-24 11:55:01 | 显示全部楼层
回复【60楼】wshtyr
高中时老师教过一个徒手开跟的方法,只要用到乘法和加减法,不知是什么原理。
-----------------------------------------------------------------------

回复【61楼】wshtyr

(原文件名:徒手开跟.JPG)

-----------------------------------------------------------------------

原理就是(a*10+b)^2=(a^2)*100+20*a*b+b^2

出0入0汤圆

发表于 2010-4-24 12:42:25 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-4-24 12:44:37 | 显示全部楼层
数论这东西,果然够神奇。

出0入0汤圆

发表于 2010-4-24 13:13:58 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-4-25 10:57:45 | 显示全部楼层
数学中的美
尽情体会吧

出0入0汤圆

发表于 2010-4-25 11:41:08 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-10-26 15:53:28 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-10-26 16:17:05 | 显示全部楼层
很好的东西

出0入9汤圆

发表于 2010-10-26 16:29:31 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-10-26 16:57:36 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-10-26 17:01:02 | 显示全部楼层
推荐大家用cordic坐标旋转算法

出0入0汤圆

发表于 2010-10-26 17:44:22 | 显示全部楼层
标记

出0入0汤圆

发表于 2010-10-26 19:48:24 | 显示全部楼层
这个当然要马克啊。

出0入0汤圆

发表于 2010-10-26 20:07:35 | 显示全部楼层
MARK

出0入0汤圆

发表于 2010-11-16 10:47:26 | 显示全部楼层
这个当然要马克啊。YE.

出0入0汤圆

发表于 2010-11-16 11:01:15 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-11-16 11:07:35 | 显示全部楼层
收藏了
非常好!

出0入0汤圆

发表于 2010-11-16 12:42:38 | 显示全部楼层
收藏

出0入0汤圆

发表于 2010-11-16 12:46:13 | 显示全部楼层
收藏

出0入0汤圆

发表于 2010-11-16 13:08:25 | 显示全部楼层
mark 数论忘得差不多了。。。

出0入0汤圆

发表于 2010-11-16 13:16:08 | 显示全部楼层
数学没学好,是个悲剧..

出0入0汤圆

发表于 2010-11-16 13:20:41 | 显示全部楼层
记号

出0入0汤圆

发表于 2010-11-16 13:37:59 | 显示全部楼层
mark,感觉数学就一个字  “变”!

出0入0汤圆

发表于 2010-12-22 14:37:15 | 显示全部楼层
挺好

出0入0汤圆

发表于 2010-12-22 15:01:21 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-12-22 15:05:22 | 显示全部楼层
简直神奇呢!

出0入0汤圆

发表于 2010-12-22 18:45:20 | 显示全部楼层
整型数的适合么?

出0入0汤圆

发表于 2010-12-22 20:08:26 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-12-22 20:49:37 | 显示全部楼层
注意误差接近2%,某些情况不能接受

出0入0汤圆

发表于 2010-12-22 21:07:04 | 显示全部楼层
真神奇

出0入0汤圆

发表于 2010-12-22 21:10:42 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-12-22 22:04:07 | 显示全部楼层
谢谢

出0入0汤圆

发表于 2011-4-7 15:34:41 | 显示全部楼层
mark

出0入0汤圆

发表于 2011-4-8 11:11:34 | 显示全部楼层
mark

出0入0汤圆

发表于 2011-4-8 12:33:17 | 显示全部楼层
回复【1楼】liuqian  刘汧
太nb了。按也作过牛顿迭代开方,不过一般是10次。
-----------------------------------------------------------------------

我做过开4次方的,但迭代4次的精度已经很高了,10次感觉没必要。
回帖提示: 反政府言论将被立即封锁ID 在按“提交”前,请自问一下:我这样表达会给举报吗,会给自己惹麻烦吗? 另外:尽量不要使用Mark、顶等没有意义的回复。不得大量使用大字体和彩色字。【本论坛不允许直接上传手机拍摄图片,浪费大家下载带宽和论坛服务器空间,请压缩后(图片小于1兆)才上传。压缩方法可以在微信里面发给自己(不要勾选“原图),然后下载,就能得到压缩后的图片】。另外,手机版只能上传图片,要上传附件需要切换到电脑版(不需要使用电脑,手机上切换到电脑版就行,页面底部)。
您需要登录后才可以回帖 登录 | 注册

本版积分规则

手机版|Archiver|amobbs.com 阿莫电子技术论坛 ( 粤ICP备2022115958号, 版权所有:东莞阿莫电子贸易商行 创办于2004年 (公安交互式论坛备案:44190002001997 ) )

GMT+8, 2024-4-26 21:47

© Since 2004 www.amobbs.com, 原www.ourdev.cn, 原www.ouravr.com

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