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

比二分更快的方法

---

如果要求一個高次方程的根，我們可以用二分法來做，這是最基礎的方法了。但是有沒有更好更快的方法呢？

 

我們先來考察一個方程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)