平方根运算的软件与硬件的加速计算
2017-09-28 16:13
323 查看
1. 平方根运算软件算法
1.1 二分法利用二分进行开平方的思想很简单:假定中值为最终解。假定下限为0,上限为x,然后求中值;然后比较中值的平方和x的大小,并根据大小修改下限或者上限;重新计算中值,开始新的循环,直到前后两次中值的距离小于给定的精度为止。需要注意的一点是,如果x小于1,我们需要将上限置为1。代码如下:
float SqrtByBisection(float n) { float low,up,mid,last; low=0,up=(n<1?1:n); mid=(low+up)/2; do { if(mid*mid>n) up=mid; else low=mid; last=mid; mid=(up+low)/2; }while(fabsf(mid-last) > eps); //eps为精度,一般为0.000001 return mid; }
这里有一点需要特别注意:在精度判别时不能利用上下限而要利用前后两次mid值,否则可能会陷入死循环!这是因为由于精度问题,在循环过程中可能会产生mid值和up或low中的一个相同。这种情况下,后面的计算都不会再改变mid值,因而在达不到精度内时就陷入死循环。
1.2 牛顿迭代法
将中值替换为切线方程的零根作为最终解,原理可以利用下图解释:
具体代码:
float SqrtByNewton(float x) { int temp = (((*(int *)&x)&0xff7fffff)>>1)+(64<<23); float val=*(float*)&temp; //初始值 float last; do { last = val; val =(val + x/val) / 2; //利用一阶梯度更新结果 }while(fabsf(val-last) > eps); //eps为精度,一般为0.000001 return val; }
注意到代码中初始值的选择,解释如下:
IEEE浮点标准用的形式来表示一个数,将该数存入float类型之后变为:
现在需要对这个浮点数进行开方,我们看看各部分都会大致发生什么变化。指数E肯定会除以2,127保持不变,m需要进行开方。由于指数部分是浮点数的大头,所以对指数的修改最容易使初始值接近精确值。幸运的是,对指数的开平方我们只需要除以2即可,也即右移一位。但是由于E+127可能是奇数,右移一位会修改指数,我们将先将指数的最低位清零,这就是& 0xff7fffff的目的。然后将该转换后的整数右移一位,也即将指数除以2,同时尾数也除以2(其实只是尾数的小数部分除以2)。由于右移也会将127除以2,所以我们还需要补偿一个64,这就是最后还需要加一个(64<<23)的原因。
这里大家可能会有疑问,最后为什么加(64<<23)而不是(63<<23),还有能不能不将指数最后一位清零?答案是都可以,但是速度都没有我上面写的快。这说明我上面的估计更接近精确值。下面简单分析一下原因。首先假设e为偶数,不妨设e=2n,开方之后e则应该变为n,127保持不变,我们看看上述代码会变为啥。e+127是奇数,会清零,这等价于e+126,右移一位变为n+63,加上补偿的64,指数为n+127,正是所需!再假设e为奇数,不妨设e=2n+1,开方之后e应该变为n+1(不精确),127保持不变,我们看看上述代码会变为啥。e+127是偶数等于2n+128,右移一位变为n+64,加上补偿的64,指数为n+1+127,也是所需!这确实说明上述的估计比其他方法更精确一些,因而速度也更快一些。
1.3 卡马克算法—快速平方根倒数算法
此法实质上是进行了1~2次的牛顿迭代法求开方倒数,求平方根的倒数,实际就是求方程1/(x^2)-a=0的解。将该方程按牛顿迭代法的公式展开为:x[n+1]=1/2*x [ n ] *(3-a*x
*x
)$。这个方法的一个关键地方还在于起始值的选取,算法选取一个“magic number”,使算法的初始解非常接近精确解!
float SqrtByCarmack( float number ) { int i; float x2, y; const float threehalfs = 1.5F; x2 = number * 0.5F; y = number; i = * ( long * ) &y; i = 0x5f375a86 - ( i >> 1 ); //得到初始化值 y = * ( float * ) &i; //第一次迭代,可以根据精度要求重复此过程,一般一次迭代精度就够用了 y = y * ( threehalfs - ( x2 * y * y ) ); return number*y; }
2. 整数平方根———硬件加速算法
2.1 手工开平方算法法先以10进制为例,解释手动开平方算法:
首先将数字每两位分成一段。如:745836942。就分成:
7|45|83|69|42,即从个位开始分。共分成五段,第一段是7。
对第一段的数字7开方取整,可得a=2。此时,要在2后面接一个数字b,并在7后面加上下一段的数45,使产生的两位数ab的平方不大于745。
我们知道,数ab实际值表示为10a+b,其平方是100a2+20ab+b2。我们可以暂时忽略b2,而产生一个“试商”b。即b=(745−100a2)/20a=(745−100∗2∗2)/(20∗2)=8.625≈8(向下取整)。但是,我们发现282=784>745,于是这个试商需要减少为7。然而,当a=0时,上述求试商的方法不在适用,但我们可以直接取下一段的两位数开方。如√45≈6。求出试商后,用745−ab2得到新的“第一段”的数。
具体过程:
取出第一段的数mn,开方得到a,然后接上第二段的数pq,用mnpq−100a2得到“余数”x,x/20a得到试商b,然后调整b(当20∗a∗b+b∗b>x时b需要减少),调整后,将x−20a∗b−b∗b作为新的余数x’,ab就是第二轮开平方的结果。由于前面已经将100a2减掉,所以后面每次都不用再减去100a2。重复步骤,直到开方完毕,或达到要求的精度为止。最后得到的a就是平方根。
如求745836942的平方根:
7|45|83|69|42
①a=sqrt(7)≈2,b=(745−400)/40≈8,282=784>745==>b=7==>272=729<745,745−729=16
②a=27,b=1683/540≈3,27∗20∗3+3∗3=1629<1683,1683−1629=54
③a=273,b=5469/5460≈1,273∗20∗1+1∗1=5461<5469==>b=1,5469−5461=8
当运用到二进制时,数ab实际值表示为a<<1+b,其平方是a2<<2+ab<<2+b2。我们可以暂时忽略b2,而产生一个“试商”b,依此类推下去。以下代码为求一个32位数的平方根,算法将原始值两两分组进行计算,注意:根应至多为16位且每次计算的b值只能是0或1:
unsigned int sqrt_16_1(unsigned int M) { unsigned int N,N2, i,j; unsigned int tmp; // 每轮计算的残差 if (M == 0) // 被开方数为0,开方结果也为0 return 0; N = 0; N2 = 0; tmp = 0; for (i=16; i>0; i-=1) //结果为16位 { N <<= 1; // 左移一位 N2 = N << 1; //N2代表a<<2 tmp <<= 2; // tmp储存的是每次计算完剩下的残差 tmp += ((M&0xc0000000) >> 30); //再加上此次计算的2bit的值 M <<= 2; if( tmp >= N2 + 1 ) //试值ab<<2+bi^2与残差的比较,此时设bi=1 { tmp -= N2 + 1; //计算此轮的残差 N++; //确定bi=1 } } return N; }
将原始值分类大小由2变成4,那么每次要计算2位二进制数值,可能情况有00,01,10,11,00情况时由于不更新结果,直接移动两位即可,所以可不进行处理,内部循环因此为1~3!
也可把分类大小变成8,16…只要是2的倍数均可,算法流程类似!
具体代码如下:
unsigned int sqrt_16_2(unsigned int M) { unsigned int N,N2, i,j; if (M == 0) // 被开方数,开方结果也为0 return 0; N = 0; N2 = 0; tmp = 0; for (i=16; i>0; i-=2) // 求剩余的15位,每次2位 { N <<= 2; // 左移两位 N2 = N << 1; // N2代表2a tmp <<= 4; // tmp储存的是每次查表完剩下的残差 tmp += (M >> 28); // 再加上下一个4bit的值 for (j=1;j<4;j++) { ttp[j] = N2*j + j*j; //求2ab+bi^2,bi可为01,10,11 } M <<= 4; for (j=3;j>0;j--) { if (tmp >= ttp[j]) //试值ab<<2+bi^2与残差的比较,bi可为01,10,11 { tmp -= ttp[j]; N+=j; break; //bi只可取一个值 } } } return N; }//32->16
相关文章推荐
- 不是运算容错,而是高温降频率,软件劣化老硬件
- Intel硬件指令加速计算CRC32
- 【并行计算-CUDA开发】CUDA软件架构与Nvidia硬件对应关系
- WPF Rendering 1(硬件加速、软件加速)
- Mathcad 是一种工程计算软件,主要运算功能:代数运算、线性代数、微积分、符号计算、2D和3D图表、动画、函数、程序编写、逻辑运算、变量与单位的定义和计算等。
- Mathcad 是一种工程计算软件,主要运算功能:代数运算、线性代数、微积分、符号计算、2D和3D图表、动画、函数、程序编写、逻辑运算、变量与单位的定义和计算等。
- 云计算加速互联网公司与软件公司的融合
- Haskell作业1(1)|实现分数的常用运算(2)|计算平方根的Newton-Raphson公式
- 3D软件加速和硬件加速
- 软件指挥硬件
- 【图像】软件使用--matlab点运算
- PHP应用加速工具软件
- 用C为密集运算函数加速
- 小学四则运算的“软件”
- 计算景深的ZEMAX宏 光学软件
- chrome如何开启硬件加速?
- 由第三方软件or硬件提供商导致的oracle database 损坏问题
- 类脑计算与神经网络加速
- 理解WebKit和Chromium: Chromium硬件加速合成
- WPF/Silverlight深度解决方案:(十八)GPU硬件加速下Silverlight超性能动画实现(下)