最近看到了这个比较有意思的题目，探究了一下。
当然有最暴力的方法，直接遍历，(0.0, x)区间内所有的数据(也可以是 x/2)，看值是否相等，但该方法太过暴力，在此不讨论。
可以采用下面的 二分法、牛顿法，Quake III源码中的方法。
浮点数的比较方法比较特别，不能使用平常的 x == 0的方法进行比较，下面都用到了equal比较方法，equal的代码：
const double Min_value = 1e-7; //定义最小的浮点数数误差，这里的用法比 宏更好
bool equal(double num1, double num2)
{
if ((num1 - num2 > -Min_value) && (num1 - num2 < Min_value))
return true;
else
return false;
}
类似于二分查找的思路，时间复杂度 O ( l o g n ) O(logn) O(logn)。
double sqrt_binary_search(double square)
{
if(square < Min_value)
return -1;
double low = 0.0;
double high = square;
// while (high - low > Min_value)
while (!equal(high, low))
{
/* code */
double mid = (high + low) / 2;
if (mid * mid > square)
{
high = mid;
}
else
{
low = mid;
}
}
return low;
}
先分析下牛顿法的原理，用的是牛顿迭代方法。
牛顿迭代法又称为牛顿-拉弗森方法，实际上是由牛顿、拉弗森各自独立提出来的。牛顿-拉弗森方法提出来的思路就是利用切线是曲线的线性逼近这个思想。
随便找一个曲线上的A点（为什么随便找，根据切线是切点附近的曲线的近似，应该在根点附近找，但是很显然我们现在还不知道根点在哪里），做一个切线，切线的根（就是和x轴的交点）与曲线的根，还有一定的距离。牛顿、拉弗森们想，没关系，我们从这个切线的根出发，做一根垂线，和曲线相交于B点，继续重复刚才的工作：
经过多次迭代后会越来越接近曲线的根。
上面就是牛顿迭代的几何直觉。只有将上述过程用代数表示出来，才方便我们写代码。
首先， x 0 x_0 x0处的切线方程为： y − y 0 = f ′ ( x ) ( x − x 0 ) y - y_0 = f'(x) (x - x_0) y−y0=f′(x)(x−x0)，即： f ( x ) = y = f ′ ( x ) ( x − x 0 ) + f ( x 0 ) f(x) = y = f'(x) (x - x_0) + f(x_0) f(x)=y=f′(x)(x−x0)+f(x0)。
要 求在 x n x_n xn所在切线方程交于x轴的点 x n + 1 x_{n+1} xn+1，即求解 0 = f ′ ( x ) ( x n + 1 − x n ) + f ( x n ) 0 = f'(x) (x_{n+1} - x_n) + f(x_n) 0=f′(x)(xn+1−xn)+f(xn)，==> x n + 1 = x n − f ( x n ) f ′ ( x n ) x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} xn+1=xn−f′(xn)f(xn)。
通过上面的公式可以得出： x n + 1 = x n − f ( x n ) f ′ ( x n ) x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} xn+1=xn−f′(xn)f(xn)，该公式可用于求任何次方的解。
这里求的是 平方根，即求解 f ( x ) = x 2 = n f(x) = x^2 = n f(x)=x2=n时，x的值是多少。也就是求解： f ( x ) = x 2 − n = 0 f(x) = x^2 - n = 0 f(x)=x2−n=0的解。
由于 f ′ ( x ) = 2 x f'(x) = 2x f′(x)=2x，所以 x n + 1 = x n − x n 2 − n 2 x n x_{n+1} = x_n - \frac{x_n^2 - n}{2x_n} xn+1=xn−2xnxn2−n，即 x n + 1 = x n + n x n 2 x_{n+1} = \frac{x_n + \frac{ n}{x_n}}{2} xn+1=2xn+xnn
实现上诉公式时，一个小改进：计算机不善于做除法，所以这里修改为乘 0.5。
double sqrt_newton(double square)
{
if (square < Min_value)
return -1;
double ret = 1.0;
while (!equal(ret*ret, square))
{
ret = (ret + square / ret) * 0.5;
}
return ret;
}
上面的初始值是随意选择的，选择的好坏一定程度上影响了计算的速度。如，下面的方法选择的初始值是 0x5f3759df。
double InvSqrt(double x)
{
double xhalf = 0.5f * x;
int i = *(int *)&x; // get bits for floating value
i = 0x5f3759df - (i >> 1); // gives initial guess y0
x = *(double *)&i; // convert bits back to float
x = x * (1.5f - xhalf * x * x); // Newton step, repeating increases accuracy
return x;
}
#include <cstdlib>
#include <cstdio>
#include <ctime>
#include <cmath>
using namespace std;
// const double Min_value = 0.0000001;
const double Min_value = 1e-7;
bool equal(double num1, double num2)
{
if ((num1 - num2 > -Min_value) && (num1 - num2 < Min_value))
return true;
else
return false;
}
double sqrt_binary_search(double square)
{
if(square < Min_value)
return -1;
double low = 0.0;
double high = square;
// while (high - low > Min_value)
while (!equal(high, low))
{
/* code */
double mid = (high + low) / 2;
if (mid * mid > square)
{
high = mid;
}
else
{
low = mid;
}
}
return low;
}
double sqrt_newton(double square)
{
if (square < Min_value)
return -1;
// double ret = 0x5f3759df; // 另一种初值
double ret = 1.0;
while (!equal(ret*ret, square))
{
ret = (ret + square / ret) * 0.5;
}
return ret;
}
double InvSqrt(double x)
{
double xhalf = 0.5f * x;
int i = *(int *)&x; // get bits for floating value
i = 0x5f3759df - (i >> 1); // gives initial guess y0
x = *(double *)&i; // convert bits back to float
x = x * (1.5f - xhalf * x * x); // Newton step, repeating increases accuracy
return x;
}
int main()
{
double x;
scanf("%lf", &x);
int start = clock();
printf("sqrt_binary_search:%lf sqrt = %lf, time = %lf\n", x,
sqrt_binary_search(x), (double)(clock() - start) / CLOCKS_PER_SEC);
start = clock();
printf("sqrt_newton:%lf sqrt = %lf, time = %lf\n", x, sqrt_newton(x),
(double)(clock() - start) / CLOCKS_PER_SEC);
start = clock();
printf("InvSqrt:%lf sqrt = %lf, time = %lf\n", x, InvSqrt(x),
(double)(clock() - start) / CLOCKS_PER_SEC);
start = clock();
printf("sqrt:%lf sqrt = %lf, time = %lf\n", x, sqrt(x),
(double)(clock() - start) / CLOCKS_PER_SEC);
}
上述代码经过测试，可以发现 来自于Quake III源码的解法 和 cmath中的 sqrt 的速度差不多。
