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

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

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

后面的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)ifactorical是阶乘项(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时跳出循环。因为sumterm都是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);
}

6 Reponses So Far ↓

  1. Tony:

    我讨厌高数

  2. 8miler:

    牛了。。我数学还挂了,看到那个拉格朗就像杀人。

  3. punkid:

    哈哈,原来高中我还自认为数学很牛X,到了大学才发现那是一坨屎啊…

    还好,这个泰勒公式现在总算知道怎么用了,那个拉格朗日就确实是很抓狂,还不知道学来干嘛的。

  4. druggo:

    我就记得《拉格朗日坟场》

  5. zcc:

    看了about
    你的性格偶喜欢。呵呵!
    college in china – failed system

  6. punkid:

    Not only the college, but the whole education system in China, from kindergarten to post-graduation, it’s all about Br&*&*ain@$#$wash.

    Like a review said:

    教育成了我国社会最大的资源浪费,高考陷入了高强度低效率的竞争

Leave a Reply ↓