使用C语言浮点数的精度极限比较精确的求近似值

今天继续学习《C语言的科学和艺术》一书，今晚终于结束了第六章算法部分的课后习题，同时我也非常高兴自己已经完成了本书第一部分的学习，很有一种收获的感觉。

今天做的题中有一个知识点，是关于浮点数精度问题的。要知道，在C语言中浮点数类型如同float，double等都是不精确的，万万不可以使用他们做科学计算，否则结果很可能与答案有很大的偏差。算法这一章中介绍了一个浮点数不精确的例子就是如果for循环中的变量是double类型的，如for(double i = 1.0; i <= 2.0; i += 0.1 )。看上去i的值会从1.0，1.1，1.2...一直取到2.0，但是在大多数电脑上是无法取值到2.0的。因为计算机中i最终的值很有可能是1.9999999...我在自己的虚拟机上验证了一下作者所举的例子，证明其正确。

今天我想说的很我第二段中所讲的有些联系，是利用浮点数的精度限制做计算（当然不可能是精确的计算了）。我举个例子吧，比如使用计算机模拟圆周率PI的计算，我们通过一个公式计算，一直计算下去，直到计算机的浮点精度达到极限。这样计算虽然不能够精确的计算出我们想要的答案，但是这个值是计算机所能计算的最近似这个精确值的值。下面我就通过一道练习题说明，这是《C语言的科学和艺术》第六章中最后一道习题。

其实这道题算法非常容易想出，因为规律非常明显。第一项为1，后面的项中分子为1分母为第n-1项的阶乘。
于是，小使了一下逐步求精的设计原则，在主函数中编写了for循环来进行项的相加，并且在其中调用Factorial函数（自定义的求阶乘函数）。请先不要看代码 。 for循环里面不就是一个sum += 1 / Factorial(n)么？结束了main函数的编写后，开始声明、定义求阶乘的函数Factorial。开始的时候，我将Factorial定义为long类型，参数也定义为long类型，但是最后程序居然报错说浮点数被0所除。我很是不解，最后我将程序改动，将原来的long类型全部修改为double型。现在想想当时的问题，难道是long转换为double型的时候会损失精度？（我曾经测试1.0/Factorial(2)，结果居然为0）(Floating point error: Divided by 0) 精度的问题很是令人烦恼。今天的重点不在这个，虽然你可能觉得我连这个都没有搞清楚而不应该往下讲下去，但是下面的内容我是清楚的。我所讲的就是关于如何将计算结果达到浮点数的精度限制的问题，请看代码18行，for循环的循环条件为sum != sum + term。这两个变量都是double性浮点变量，其中term的值是一个越来越小的变量（term永远大于0，因为term = 1 / f, f又大于0 ）因此，必然会有一个时刻term将等于0，那是sum = sum + term，循环停止。

收获：学会使用浮点数的精度限制来“比较精确”的计算近似值，使用这一方法在许多计算中非常重要。

代码如下：

<!--<br /> <br /> Code highlighting produced by Actipro CodeHighlighter (freeware)<br /> http://www.CodeHighlighter.com/<br /> <br /> -->  1  /*
 2      The Art and Sience of C.
 3     Issue 6, project 10.
 4     Title: e
 5  */
 6  
 7  #include < stdio.h >
 8  #include < math.h >
 9  #include " simpio.h "
10  double  Factorial(  double  n );
11  
12  void  main()
13  {
14       double  e  =   1.0 ;
15       double  term;
16      int  i;
17  
18       for  ( i  =   1 ; e  !=  e  +  term; i ++  )
19      {
20           double  f  =  ( double )Factorial( i );
21          term  =   1.0   /  f;
22          e  +=  term;
23     }
24  
25      printf( " %.10lf\n " ,e);
26  
27      getch();
28      printf( " \n " );
29  }
30  
31  double  Factorial(  double  n )
32  {
33       if ( n  <=   0.0  )
34     {
35        printf( " Illegal parameter!\n " );
36         return   - 1 ;
37     }
38       if ( n  ==   1.0  )    return (  1  );
39       if ( n  >   1.0  )     return ( n  *  Factorial( n  -   1  ) );
40  }