对数学(尤其是大学高等数学)有恐惧心理的请略过吧…
Algorithms (算法) 这一章有个运用泰勒公式(Taylor Series)求平方根(square root)的程序。虽然C语言的math.h库里已经有个sqrt的function了，不过这个sqrt function本身怎么写呢? 大学微积分里的泰勒公式终于派上用途了…我正好忘的一干二净了:)
先来学习/温习下泰勒公式吧:
后面的Rn(x) 是拉格朗日余项，当|x-x0|无穷小时，可以去掉拉格朗日余项，近似认为f(x) 等于前面的n项之和。
要求的是f(x) = √x，根据求导公式f'(X) = c Xc-1，得出f'(x0) = 1/2 x-1/2，二阶，三阶…N阶导以此类推。假设x0 = 1，则f(1) = 1, f'(1) = 1/2, f''(1) = - 1/4, f'''(1) = 3/8 ...。 将每项进行分拆，基本可以分为3个部分
xpower
coeff ------------
factorical
其中的coeff是导数项fi(x0=1)，xpower是(x-1)i，factorical是阶乘项(i+1)! (i = 0,1,2…)
推到可知:
coeff *= (0.5 - i);
xpower *= (x - 1);
factorial *= (i + 1);
所以用泰勒公式求平方根就可以可以编写C程序了
double TSqrt(double x)
{
double sum, factorial, coeff, term, xpower;
int i;
factorial = coeff = xpower = 1;
sum = 0;
term = 1;
for (i=0; sum != sum + term; i++) {
sum += term;
coeff *= (0.5 - i);
xpower *= (x -1);
factorial *= (i + 1);
term = coeff * xpower / factorial;
}
return (sum);
}
泰勒公式自然是要求i值越大，结果越精确，但是这就陷入了C语言中的无限循环，程序永远不会停止也得不到结果。为了限定for循环的运算次数，可以设定条件为sum != sum + term，即当sum = sum + term时跳出循环。因为sum和term都是double型，其值是不精确的，所以随着循环次数的不断进行，后面的项数值term会越来越趋近于无穷小，例如当sum = 23，而term = 0.000000000001时，计算机会认为 23 (sum) = 23.000000000001 (sum + term)。
不过这个程序还有个问题，就是当假设x0 = 1时，泰勒公式要求x的取值范围是:
0 < x <2否则随着x的增大，运算结果会越来越失真。
解决办法的思路是: 因为 √4x = √4 √x = 2√x，我们可以把大数通过不断的除4最终化为(0,2)区间的小数。只要记得除了多少次4，最后在结果上乘上这个次数的2就行了。
double Sqrt(double x)
{
double correction;
if (x == 0) return (0);
if (x < 0) Error("Sqrt called with negative argument %g", x);
correction = 1;
while (x >=2) {
x /= 4;
correction *= 2;
}
return (TSqrt(x) * correction);
}
特别要注意的是，x是double型，才能保证除以4后的数为带小数位的较精确值，而不是去尾的整数。想象如果x是int型，那么16和17的运算结果会是一样的。
OK，整个运用泰勒公式求平方根的程序就出来了:
/*
* Function : Sqrt
* Usage: root = Sqrt(x);
* -------------------------------
* Description: Returns the square root of x.
* /
double Sqrt(double x)
{
double correction;
if (x == 0) return (0);
if (x < 0) Error("Sqrt called with negative argument %g", x);
correction = 1;
while (x >=2) {
x /= 4;
correction *= 2;
}
return (TSqrt(x) * correction);
}
/*
* Function : TSqrt
* Usage: root = TSqrt(x);
* ---------------------------------
* Description: Returns the square root of x in the range 0 < x < 2.
* /
double TSqrt(double x)
{
double sum, factorial, coeff, term, xpower;
int i;
factorial = coeff = xpower = 1;
sum = 0;
term = 1;
for (i=0; sum != sum + term; i++) {
sum += term;
coeff *= (0.5 - i);
xpower *= (x -1);
factorial *= (i + 1);
term = coeff * xpower / factorial;
}
return (sum);
}