常用的平方根算法详解与实现

本文从属于笔者的数据结构与算法系列文章。

SquareRoot

平方根计算一直是计算系统的常用算法，本文列举出几张简单易懂的平方根算法讲解与实现。其中Java版本的代码参考这里

Reference

• 计算平方根的算法

• Wiki-Methods of computing square roots

Babylonian:巴比伦算法/牛顿法

巴比伦算法可能算是最早的用于计算$sqrt{S}$的算法之一，因为其可以用牛顿法导出，因此在很多地方也被成为牛顿法。其核心思想在于为了计算x的平方根，可以从某个任意的猜测值g开始计算。在真实的运算中，我们往往将g直接设置为x，不过也可以选择其他任何的正数值。那么其计算的迭代过程为:

1.如果猜测值g已经足够接近于正确的平方根，算法结束，函数将g作为结果返回。

2.如果猜测值g不够精确，那么使用g和x/g的平均值作为新的猜测值。因为这两个值中的一个小于确切的平方根，另一个则大于确切的平方根，选择平均值有助于你得到一个更接近于正确答案的值。

3.把新的猜测值赋予给变量g，重复第一步的判断。

综上所述，其计算公式可以表述为:

其实现的参考代码地址为SquareRoots:

public double Babylonian() {

    double g = this.value;

    while (isApproximate(g)) {
        g = (g + this.value / g) / 2;
    }

    return g;

}

基于泰勒公式的级数逼近

微积分中的泰勒级数可以表示为:

在这个公式中，符号a表示某个常量，记号f'、f''和f'''表示函数f的一阶、二阶和三阶导数，以此类推，这个公式称为泰勒公式，基于这个公式，我们平方根公式的展开式为:

根据该公式我们可以在一定精度内逼近真实值，不过这个公式仍然存在一个问题，即是公式的收敛问题。在泰勒级数展开中，平方根函数的公式当且仅当参数值位于一个有效范围内时才有效，在该范围内计算趋于收敛。该范围即是收敛半径，当我们对平方根函数用a=1进行计算时，泰勒级数公式希望x处于范围:$0<x<2$之间。如果x在收敛半径之外，则展开式中的项会越来越大，泰勒级数离答案也就越来越远。为了解决该问题，我们可以考虑当待开平方数大于4时以4去除它，最后将得到的数乘以相同次数的2即可。其实现的参考代码地址为SquareRoots:

public double TSqrt() {

    //设置修正系数
    double correction = 1;

    //因为要对原值进行缩小,因此设置临时值
    double tempValue = value;

    while (tempValue >= 2) {
        tempValue = tempValue / 4;
        correction *= 2;
    }

    return this.TSqrtIteration(tempValue) * correction;
}

private double TSqrtIteration(double value) {

    double sum = 0, coffe = 1, factorial = 1, xpower = 1, term = 1;

    int i = 0;

    while (Math.abs(term) > 0.000001) {

        sum += term;

        coffe *= (0.5 - i);

        factorial *= (i + 1);

        xpower *= (value - 1);

        term = coffe * xpower / factorial;

        i++;
    }

    return sum;

}

平方根倒数速算法

首先接收一个32位带符浮点数，然后将之作为一个32位整数看待，以将其向右进行一次逻辑移位的方式将之取半，并用十六进制“魔术数字”0x5f3759df减之，如此即可得对输入的浮点数的平方根倒数的首次近似值；而后重新将其作为浮点数，以牛顿法反复迭代，以求出更精确的近似值，直至求出符合精确度要求的近似值。在计算浮点数的平方根倒数的同一精度的近似值时，此算法比直接使用浮点数除法要快四倍。此算法最早被认为是由约翰·卡马克所发明，但后来的调查显示，该算法在这之前就于计算机图形学的硬件与软件领域有所应用，如SGI和3dfx就曾在产品中应用此算法。而就现在所知，此算法最早由Gary Tarolli在SGI Indigo的开发中使用。虽说随后的相关研究也提出了一些可能的来源，但至今为止仍未能确切知晓此常数的起源。

其实现的参考代码地址为SquareRoots:

public double FastInverseSquareRoot() {

    double tempValue = value;

    double xhalf = 0.5d * tempValue;

    long i = Double.doubleToLongBits(tempValue);

    i = 0x5fe6ec85e7de30daL - (i >> 1);

    tempValue = Double.longBitsToDouble(i);

    tempValue = tempValue * (1.5d - xhalf * tempValue * tempValue);

    tempValue = this.value * tempValue;

    return tempValue;

}

Comparsion:比较

笔者建立了一个专门的单元测试类来比较上述算法的准确度与性能，代码参考SquareRootsTest，首先在准确度与稳定性测试方面，这几种算法都能达到较好地稳定性，其中平方根倒数速算法相对而言是较好。

@Test
public void testBabylonian() {

    for (int i = 0; i < 10000; i++) {
        Assert.assertEquals(2.166795861438391, squareRoots.Babylonian(), 0.000001);

    }

}

@Test
public void testTSqrt() {

    for (int i = 0; i < 10000; i++) {

        Assert.assertEquals(2.166795861438391, squareRoots.TSqrt(), 0.000001);

    }
}

@Test
public void testFastInverseSquareRoot() {

    for (int i = 0; i < 10000; i++) {

        Assert.assertEquals(2.1667948388864198, squareRoots.FastInverseSquareRoot(), 0.000001);

    }
}

而在性能测试方面，级数逼近的性能最差，巴比伦算法次之，平方根倒数速算法最好:

@Test
public void benchMark() {

    //巴比伦算法计时器
    long babylonianTimer = 0;

    //级数逼近算法计时器
    long tSqrtTimer = 0;

    //平方根倒数速算法计时器
    long fastInverseSquareRootTimer = 0;

    //随机数生成器
    Random r = new Random();

    for (int i = 0; i < 100000; i++) {

        double value = r.nextDouble() * 1000;

        SquareRoots squareRoots = new SquareRoots(value);

        long start, stop;

        start = System.currentTimeMillis();

        squareRoots.Babylonian();

        babylonianTimer += (System.currentTimeMillis() - start);

        start = System.currentTimeMillis();

        squareRoots.TSqrt();

        tSqrtTimer += (System.currentTimeMillis() - start);

        start = System.currentTimeMillis();

        squareRoots.FastInverseSquareRoot();

        fastInverseSquareRootTimer += (System.currentTimeMillis() - start);

    }


    System.out.println("巴比伦算法:" + babylonianTimer);

    System.out.println("级数逼近算法:" + tSqrtTimer);

    System.out.println("平方根倒数速算法:" + fastInverseSquareRootTimer);


}


/**
结果为:
巴比伦算法:17
级数逼近算法:34
平方根倒数速算法:7
**/