牛顿迭代法求高次方程的根
比二分更快的方法
如果要求一个高次方程的根，我们可以用二分法来做，这是最基础的方法了。但是有没有更好更快的方法呢？
我们先来考察一个方程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:
