The Art and Science of C 阅读笔记 II

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

Algorithms (算法) 这一章有个运用泰勒公式(Taylor Series)求平方根(square root)的程序。虽然C语言的math.h库里已经有个sqrt的function了，不过这个sqrt function本身怎么写呢? 大学微积分里的泰勒公式终于派上用途了…我正好忘的一干二净了:)

先来学习/温习下泰勒公式吧:

后面的Rn(x) 是拉格朗日余项，当|x-x₀|无穷小时，可以去掉拉格朗日余项，近似认为f(x) 等于前面的n项之和。

要求的是f(x) = √x，根据求导公式f'(X) = c Xc-1，得出f'(x0) = 1/2 x-1/2，二阶，三阶…N阶导以此类推。假设x₀ = 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)。

不过这个程序还有个问题，就是当假设x₀ = 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);
}