SakuraTears的博客

越努力越幸运

编程与数学 - 求平方根

发表于 更新于 分类于 阅读次数： Valine：
本文字数： 6.5k 阅读时长 ≈ 6 分钟

前言

前几天看到一道题目，是关于“对整数 n 开平方，不使用Math.sqrt实现”，感觉蛮有意思的，其解法用到了牛顿迭代法（Newton’s Method）。

就顺便研究了一下该解法和其他一些解法，特来分享记录一下。

正文

题目是非常容易理解的，我们直接来看相关解法吧。

牛顿迭代法

因为对没有求根公式的函数，求解它的零点是非常困难的，因此就发明了 牛顿迭代法（Newton‘s Method） 来逼近该函数的零点。具体方法如下图所示:

upload successful

设 $r$ 是 $f(x)=0$ 的根，选取 $x_0$ 作为 $r$ 的初始近似值，过点 $(x_0,f(x_0))$ 做曲线 $y=f(x)$ 的切线 $L$，$L:y=f(x_0)+f’(x_0)(x-x_0)$ ，则 $L$ 与 $x$ 轴交点的横坐标 $x_1=x_0-\frac{f(x_0)}{f’(x_0)}$，称 $x_1$ 为 $r$ 的一次近似值。过点 $(x_1,f(x_1))$ 做曲线 $y=f(x)$ 的切线，并求该切线与 $x$ 轴交点的横坐标 $x_2=x_1-\frac{f(x_1)}{f’(x_1)}$，称 $x_2$ 为 $r$ 的二次近似值。重复以上过程，得 $r$ 的近似值序列，其中，$x_{n+1}=x_n-\frac{f(x_n)}{f’(x_n)}$ 称为 $r$ 的 $n+1$ 次近似值，上式称为牛顿迭代公式。

用牛顿迭代法解非线性方程，是把非线性方程 $f(x)=0$ 线性化的一种近似方法。把 $f(x)$ 在点 $x_0$ 的某邻域内展开成泰勒级数
$$f(x)=f(x_0)+f’(x_0)(x-x_0)+\frac{f’’(x_0)(x-x_0)^2}{2!}+…+\frac{f^{(n)}(x_0)(x-x_0)^n}{n!}+R_n(x)$$
，取其线性部分（即泰勒展开的前两项），并令其等于 $0$，即 $f(x_0)+f’(x_0)(x-x_0)=0$，以此作为非线性方程 $f(x)=0$ 的近似方程，若 $f’(x)\not ={0}$，则其解为 $x_1=x_0-\frac{f(x_0)}{f’(x_0)}$， 这样，得到牛顿迭代法的一个迭代关系式： $x_{n+1}=x_n-\frac{f(x_n)}{f’(x_n)}$。

已经证明，如果是连续的，并且待求的零点是孤立的，那么在零点周围存在一个区域，只要初始值位于这个邻近区域内，那么牛顿法必定收敛。 并且，如果不为0, 那么牛顿法将具有平方收敛的性能。 粗略的说，这意味着每迭代一次，牛顿法结果的有效数字将增加一倍。

说了这些，那牛顿迭代法为什么会跟 “对整数n开平方” 有关呢？

若我们另 $f(x)=x^2-n$，则 $f(x)$ 的零点即为 $\sqrt{n}$ ，此时 $f’(x)=2x$，则迭代公式如下：

$$x_{n+1}=x_n-\frac{f(x_n)}{f’(x_n)}=x_n-\frac{x_n^2-n}{2x_n}=\frac{x_n^2+n}{2x_n}=\frac{1}{2}(x_n+\frac{n}{x_n})$$

相关代码实现如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
/**
 * 牛顿迭代法
 * @param n 要开方的数字，需要大于等于0
 * @return
 */
public static double sqrt(double n) {
    if (n < 0) {
        return Double.NaN;
    }
    //精度
    double e = 1e-15;
    double x = n;
    double y = (x + n / x) / 2;
    //当x_n 与 x_n+1 精度小于e，认为逼近0点，返回符合要求的数据
    while (Math.abs(x - y) > e) {
        x = y;
        y = (x + n / x) / 2;
    }
    return x;
}
public static void main(String[] args) {
    System.out.println(sqrt(4));
    System.out.println(sqrt(3));
}

输出：

1
2
2.0
1.7320508075688772

二分查找法

相比于快速的牛顿迭代法，二分查找法也是可以实现开方需求的。不过相比牛顿迭代法其速度较慢。

这种方法十分好理解，就是上界初始化为数字本身，下界初始化为0.0，这样用二分，判断中间数字的平方和目标数字比较，再修改上界和下界，直到小于一定的阈值。需要注意结束条件和精度判断。

相关代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
/**
 * 二分开方法
 * @param n
 * @return
 */
public static double sqrt1(double n) {
    double left = 0.0;
    double hight = n;
    // 此处为精度，当满足该精度时返回近似值
    double p = 1e-15;
    double mid = (left + hight)/2.0;
    // 精度比较
    while (Math.abs(mid * mid - n) > p){
        if(mid * mid > n) {
            hight = mid;
        }else if (mid * mid < n) {
            left = mid;
        }else {
            return mid;
        }
        mid = (left + hight)/2.0;
    }
    return mid;
}


public static void main(String[] args) {
    System.out.println(sqrt1(4));
    System.out.println(sqrt1(3));
}

输出：

1
2
2.0
1.7320508075688772

Java源码中的开方

开始一直以为Java源码中Math.sqrt方法使用的是牛顿迭代法来实现的。最近研究源码，发现并不是这样。特地研究记录一下。

我们跟踪Math.sqrt源码，会发现它其实调用的StrictMath的sqrt方法，此方法为native方法。

1
public static native double sqrt(double a);

因此我们需要查找它的具体实现了，这就需要找到openjdk源码了，其实现位于openjdk\jdk\src\share\native\java\lang\fdlibm\src\e_sqrt.c路径下。

当然我们也可以在线查看，文件如下 e_sqrt.c。

我们可以看到对于开方的操作，其源码实际是使用了一种叫Bit by bit method，我这儿称为逐位法。

upload successful

根据上图，我们来看下该方法的优势及特点吧。

Bit by bit method using integer arithmetic. (Slow, but portable)

源码中提到该方法虽然“慢”但合适，其相关原理如下。

归一化

在 $[1,4)$ 中以 $2$ 的偶数次幂缩放 $x$ 到 $y$ :

求一个整数 $k$，使 $1 \leq (y=x * 2^{2k}) < 4$ ，

即 $\sqrt{x} = 2^k * \sqrt{y}$

逐位计算

设 $q_i =\sqrt{y}$ 在二进制点($q_0 = 1$)后截断到 $i$ 位，

$s_i=2q_i$ 并且 $y_i=2^{i+1}(y-q_i^2)$.       (1)

要从 $q_{i+1}$ 计算 $q_i$，首先要检查是否

$(q_i+2^{-(i+1)})^2 \leq y$.       (2)

如果 (2) 式结果为 false，就有 $q_{i+1}=q_i$，否则 $q_{i+1}=q_i+2^{-(i+1)}$.

通过一些代数运算，不难看出 (2) 式等价于

$s_i+2^{-(i+1)} \leq y_i$.       (3)

变为 (3) 式的优点是，$s_i$ 和 $y_i$ 可以用递归式计算:

如果 (3) 式为false

$s_{i+1}=s_i$，$y_{i+1}=y_i$；       (4)

否则

$s_{i+1}=s_i+2^{-i}$，$y_{i+1}=y_i-s_i-2^{-(i+1)}$；       (4)

用归纳法可以很容易地证明 (4) 和 (5)。因为 (3) 的左边只包含 $i+2$ 位，所以 (3) 中不需要进行完整的(53-bit)比较。

最终

在生成53位结果后，我们再计算一个位。连同余数，我们可以确定结果是正确的，大于1/2ulp，还是小于1/2ulp(它永远不会等于1/2ulp)。

四舍五入可以通过检查 huge + tiny 是否等于 huge，以及对于某个浮点数“huge”和“tiny”，huge - tiny 是否等于 huge 来检测。

上面的算法我们可以通过一个简单例子来理解。

假设现在我们要求 $\sqrt{36}$ 的根。

根据第一步得到，$1 \leq y = 36*2^{-4}=2.25_{2(10.01)}<4$。

迭代：

$q_0=1_{2(1)}$，初始化

$q_1=1.5_{2(1.1)}$，$(1+0.5)^2 \leq 2.25$

$q_2=1.5_{2(1.1)}$，$(1.5+0.25)^2 > 2.25$

$q_3=1.5_{2(1.1)}$，$(1.5+0.125)^2 > 2.25$

……

最终 $\sqrt{y} = 1.5$， $\sqrt{36}=1.5*2^2= 6$。

可以看出这种方法不同于牛顿迭代法，它将带求解的数映射于 $[1,4)$ 范围内，通过逐位计算，逐步缩小解的精度，逼近结果。

因为以上迭代过程涉及到平方的操作，为了优化这一点，在逐位计算这一步，使用归纳法消除了平方操作。

可以知道这种逐位计算的方法求解收敛速度某些时候或许比不上牛顿迭代法，但避免了许多乘法和除法操作，所以鲁棒性很好。

源码处理的double数据有52位有效位，在处理时将其分成了高位和低位分开处理，涉及到许多位运算，我们这儿不做详细讨论。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
#include "fdlibm.h"

#ifdef __STDC__
static  const double    one     = 1.0, tiny=1.0e-300;
#else
static  double  one     = 1.0, tiny=1.0e-300;
#endif

#ifdef __STDC__
        double __ieee754_sqrt(double x)
#else
        double __ieee754_sqrt(x)
        double x;
#endif
{
        double z;
        int     sign = (int)0x80000000;
        unsigned r,t1,s1,ix1,q1;
        int ix0,s0,q,m,t,i;

        ix0 = __HI(x);                  /* high word of x */
        ix1 = __LO(x);          /* low word of x */

    /* take care of Inf and NaN */
        if((ix0&0x7ff00000)==0x7ff00000) {
            return x*x+x;               /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
                                           sqrt(-inf)=sNaN */
        }
    /* take care of zero */
        if(ix0<=0) {
            if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */
            else if(ix0<0)
                return (x-x)/(x-x);             /* sqrt(-ve) = sNaN */
        }
    /* normalize x */
        m = (ix0>>20);
        if(m==0) {                              /* subnormal x */
            while(ix0==0) {
                m -= 21;
                ix0 |= (ix1>>11); ix1 <<= 21;
            }
            for(i=0;(ix0&0x00100000)==0;i++) ix0<<=1;
            m -= i-1;
            ix0 |= (ix1>>(32-i));
            ix1 <<= i;
        }
        m -= 1023;      /* unbias exponent */
        ix0 = (ix0&0x000fffff)|0x00100000;
        if(m&1){        /* odd m, double x to make it even */
            ix0 += ix0 + ((ix1&sign)>>31);
            ix1 += ix1;
        }
        m >>= 1;        /* m = [m/2] */

    /* generate sqrt(x) bit by bit */
        ix0 += ix0 + ((ix1&sign)>>31);
        ix1 += ix1;
        q = q1 = s0 = s1 = 0;   /* [q,q1] = sqrt(x) */
        r = 0x00200000;         /* r = moving bit from right to left */

        while(r!=0) {
            t = s0+r;
            if(t<=ix0) {
                s0   = t+r;
                ix0 -= t;
                q   += r;
            }
            ix0 += ix0 + ((ix1&sign)>>31);
            ix1 += ix1;
            r>>=1;
        }

        r = sign;
        while(r!=0) {
            t1 = s1+r;
            t  = s0;
            if((t<ix0)||((t==ix0)&&(t1<=ix1))) {
                s1  = t1+r;
                if(((t1&sign)==sign)&&(s1&sign)==0) s0 += 1;
                ix0 -= t;
                if (ix1 < t1) ix0 -= 1;
                ix1 -= t1;
                q1  += r;
            }
            ix0 += ix0 + ((ix1&sign)>>31);
            ix1 += ix1;
            r>>=1;
        }

    /* use floating add to find out rounding direction */
        if((ix0|ix1)!=0) {
            z = one-tiny; /* trigger inexact flag */
            if (z>=one) {
                z = one+tiny;
                if (q1==(unsigned)0xffffffff) { q1=0; q += 1;}
                else if (z>one) {
                    if (q1==(unsigned)0xfffffffe) q+=1;
                    q1+=2;
                } else
                    q1 += (q1&1);
            }
        }
        ix0 = (q>>1)+0x3fe00000;
        ix1 =  q1>>1;
        if ((q&1)==1) ix1 |= sign;
        ix0 += (m <<20);
        __HI(z) = ix0;
        __LO(z) = ix1;
        return z;
}

其他开方算法

openjdk 源码中除了上述方法外，在注释中还提供了两种开方方法，第一种部分使用了牛顿迭代，涉及四个部分。第二种方法使用reciproot迭代来避免除法，但是需要更多的乘法。

感兴趣的同学可以查看

IEEE754.PDF

reciprt.pdf

等论文进行了解。

总结

本篇文章我们了解了一些开方方法，并分析了一些源码，可以发现数学与编程及算法的巧妙之处。对我们今后的工作学习都是有较大帮助的。

参考资料

e_sqrt.c

Berkeley | EECS




-------------文章结束啦 ~\(≧▽≦)/~ 感谢您的阅读-------------

您的支持就是我创作的动力！

欢迎关注我的其它发布渠道