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

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[n]*x[n])$。这个方法的一个关键地方还在于起始值的选取，算法选取一个“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，其平方是$100a^2+20ab+b^2$。我们可以暂时忽略$b^2$，而产生一个“试商”b。即$b = (745 - 100a^2) / 20a = (745 - 100*2*2) / (20 * 2) = 8.625 ≈ 8$(向下取整)。但是，我们发现$28^2 = 784>745$，于是这个试商需要减少为7。然而，当a=0时，上述求试商的方法不在适用，但我们可以直接取下一段的两位数开方。如√45≈6。求出试商后，用$745-ab^2$得到新的“第一段”的数。
具体过程：
取出第一段的数mn，开方得到a，然后接上第二段的数pq，用$mnpq-100a^2$得到“余数”x，$x/20a$得到试商b，然后调整b（当$20*a*b+b*b>x$时b需要减少），调整后，将$x-20a*b-b*b$作为新的余数x’，ab就是第二轮开平方的结果。由于前面已经将$100a^2$减掉，所以后面每次都不用再减去$100a^2$。重复步骤，直到开方完毕，或达到要求的精度为止。最后得到的a就是平方根。
如求$745836942$的平方根：
$7|45|83|69|42$
$①a=sqrt(7)≈2，b=(745-400)/40≈8，28^2=784>745 ==> b=7 ==> 27^2=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，其平方是$a^2<<2+ab<<2+b^2$。我们可以暂时忽略$b^2$，而产生一个“试商”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