牛顿迭代法求高次方程的根

比二分更快的方法

---

如果要求一个高次方程的根，我们可以用二分法来做，这是最基础的方法了。但是有没有更好更快的方法呢？

 

我们先来考察一个方程f(x)的在点a的泰勒展开，展开到一阶就可以了（假设f(x)在点a可以泰勒展开，也就是泰勒展开的那个余项在n趋于无穷时趋于0）

f(x) -> f(a) + (x - a)f'(a)

现在我们令这个一阶展开为0，当f'(a)是非0值，移项一下就有

x = a - f(a)/f'(a)

 

实际上当我们把f(a)改成f(x) - f(a)，这就是一个过了f(a)的关于f(x)的切线方程，如果我们令f(x) = 0就可以得到这条切线在x轴上的交点了。

重复这个过程，我们就可以得到逼近我们所想要的答案的解了。这就是牛顿迭代法的原理。

如上图，当我们求f[x] = x^2 - 2这个方程的根时，我们可以先猜这个解是4，然后在(4,f(4))这个点上做切线交x轴零点9/4，然后我们继续在9/4这个点做切线，继续逼近。

 

当然代码写起来是很简单的了，比如我们要算简单的y = sqrt(b)，变形一下就是y^2 - b == 0的解，比如在leetcode写一下

#include <iostream>
#include <algorithm>

class Solution 
{
public:
    int mySqrt(int x) 
    {
        double result = x;//随便猜一个
        double precision = 1e-6;
        double tmpX;

        if (x == 0)
            return 0;

        do
        {
            tmpX = result;
            result = result / 2.0 + x / 2.0 / result;

        } while (std::abs(tmpX - result) >= precision);

        return result;
    }
};

写了一个matlab来和二分法的速度进行了一下对比：

测试代码：

clear
sum = 100000000;

rb = zeros(sum,2);
rn = zeros(sum,2);
index = 1;

for n = 1:sum
     s =tic;
        Newton(n);
        elasped= toc(s);
        rb(index,1) = n;
        rb(index,2) = elasped*1e4;
        
        s =tic;
        Binary(n);
        elasped= toc(s);
        
        rn(index,1) = n;
        rn(index,2) = elasped*1e4;
        index = index + 1;
end

x1 =zeros(sum,1);
y1 =zeros(sum,1);

x2 = zeros(sum,1);
y2 = zeros(sum,1);

for n = 1: size(rb)
    x1(n) = rb(n,1);
    y1(n) = rb(n,2);
end

for n = 1: size(rn)
    x2(n) = rn(n,1);
    y2(n) = rn(n,2);
end

plot(x1, y1,'*', x2, y2,'+')
set(gca,'YTick')
title('Newton Method (+)和 Binary Search (*) 算sqrt(n)(n->1~100000000)的时间比较')


function result = Newton(x)
    result = x;
    while 1
        tmpX = result;
        result = result / 2.0 + x / 2.0 / result;
        
        if abs(tmpX - result) < 1e-6
            break
        end
    end
end

function result = Binary(x)
    left = 0;
    right = x;
    
    while 1
        result = left + (right - left)/2;
        if result*result - x > 0
            right = result;
        elseif result*result - x < 0
            left = result;
        else
            return
        end
        if abs(right - left) > 1e-6
            break
        end
    end  
end

测试结果：

不过由于是在matlab环境下，速度上感觉其实没什么太大差别。

 

更快速的开方

---

QUAKE-III的源码有一段非常厉害的开方代码：

float Q_rsqrt(float number) {
        long i; float x2, y; const float threehalfs = 1.5F;
        x2 = number * 0.5F;
        y = number;
        i = *(long *)&y; // evil floating point bit level hacking 
        i = 0x5f3759df - (i >> 1); // what the fuck? 
        y = *(float *)&i;
        y = y * (threehalfs - (x2 * y * y)); // 1st iteration 
                                             // y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed
        return y; 
    }

这段代码惊为天人，开方连迭代都不用，直接出结果，那个神奇的0x5f3759df是线性拟合出来的，如果想结果更准，把y = y * ( threehalfs - ( x2 * y * y ) );这句话多写几下就好了。

 

不过，这一段代码是过不了leetcode的某个case，不过这已经不重要了。

 

Reference: 

牛顿迭代法：介绍、原理与运用

牛顿迭代法(Newton's Method)