平方根运算的软件与硬件的加速计算

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