对数学 ( 尤其是大学高等数学 ) 有恐惧心理的请略过吧…

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);
}