2.2　浮点数据类型

浮点类型是实数在计算机上的近似实现。当需要对数据的小数部分进行计算时，或者数据的绝对值非常大，以至超过整型数据的表示范围时，就需要使用浮点数。

2.2.1　浮点数据的表示方法

浮点类型数据的表示方法不同于整型。在浮点数的表示中，二进制位被分成不同的域，以分别表示浮点数的符号，数值的有效位，以及小数点的位置。目前所有常见的计算机系统都采用IEEE关于浮点数的标准IEEE-754来表示浮点数和进行浮点运算。

根据IEEE-754，一个浮点数被表示成V=(-1) ^s ×M×2 ^E 的形式，其中s表示符号，占一位二进制位。当s为0时数值是正数，s为1时数值是负数。M表示尾数，也就是浮点数的有效数字，是一个介于0和1之间的二进制小数。E表示指数，可以是正数，也可以是负数，用以确定尾数的权值，也就是规定小数点所在的位置。对于单精度浮点类型（float），尾数占23位，指数占8位。对于双精度浮点类型（double），尾数占52位，指数占11位。单精度浮点数和双精度浮点数在计算机中的表示如图2-3所示。

浮点数可以分为规格化的浮点数和非规格化的浮点数。当一个浮点数指数部分的各位既不是全0也不是全1时，这个浮点数就是一个规格化的浮点数。这时，浮点数指数部分的值等于指数部分按照无符号数解释，然后再减去偏移量。单精度浮点数的偏移量是127，双精度浮点数的偏移量是1023，因此float的指数取值范围是-126～127，而double的指数取值范围是-1022～1023。当一个浮点数指数部分的各位是全0时，这个数就是一个非规格化的浮点数。在非规格化的情况下，浮点数指数的值等于1减去偏移量。此外，在规格化的表示方式下，在尾数M的左端隐含有一个整数1，因此浮点数尾数部分的值实际上等于1.M。而在非规格化的表示方式下，在尾数M的左端没有隐含的整数1，因此浮点数尾数部分的值等于0.M。这样，规格化的浮点数尾数的精度就比非规格化的浮点数多了一位二进制位。非规格化数的用处有二：一是表示0.0；二是表示非常接近于0.0的数。因为在规格化的浮点数的尾数的左端隐含有一个整数1，所以无法表示尾数全为0的数。而非规格化浮点数因为没有这个隐含的整数1，所以可以表示尾数全为0的数。浮点数中的0.0有两种：一种是+0.0，另一种是-0.0。+0.0的所有二进制位均为0，-0.0的符号位为1，其余的所有二进制位均为0。当一个浮点数的指数部分是全1并且小数部分是全0时，这个浮点数表示无穷大。无穷大的符号由符号位决定。当一个浮点数的指数部分的二进制位是全1并且小数部分不为全0时，这个浮点数表示一个被称为NaN的特殊值。NaN表示“不是一个数”（Not a Number）。当一些计算的结果既不是实数也不是正负无穷时，就会产生这个值。

采用IEEE-754标准表示浮点数有一个重要的优点，即对浮点数大小的比较可以直接采用整数比较的方式进行。首先，符号位处于最高位，与整数中符号所在的位置一样。其次，指数部分安排在尾数部分的左侧，具有更高的权重，符合浮点数中指数对绝对值的大小影响更大的事实。最后，指数部分使用偏移量，使得这部分被表示为永远大于等于0的无符号数，因此在按整数二进制位解释的情况下也能够正确地区分不同数值的大小。这样，除了当参与比较的浮点数中有NaN的情形外，采用整数比较的方式对浮点数进行比较的结果也是正确的。规格化和非规格化浮点数的表示范围如表2-1所示：

尽管浮点数所表示的数值范围很大，但是由于尾数的位数有限，其所能表示的数值精度是有限的。对于float类型，其23位二进制尾数大约相当于7位十进制的有效数字；double类型的52位二进制尾数大约相当于16位十进制的有效数字。因此对于绝大多数实数来说，C语言的浮点表示只是一种近似的表示。对于不带小数部分的有理数，只要其有效数字的位数不超过尾数的位数，就都是可以准确表达的。对于带有小数部分的有理数，C语言浮点数所能精确表示的只是那些能够在浮点数各字段有效数字范围内表示成整数与2的整数次幂乘积的有理数，即能够表示成 X * 2 ^y 形式的数。其中 X 是能够在尾数有效数字位数的范围内准确表达的整数， y 是在指数范围内所能准确表示的整数。例如，1.5、0.625等都是可以准确表示的浮点数，因为1.5=3 * 2 ^-1 ，0.625=5 * 2 ^-3 。0.1是不能被准确表示的浮点数，这是由于0.1=1/10，而10不能表示成为2的整数次幂，因此没有任何正整数 y 能使得0.1=1 * 2 ^{- y} 成立。

2.2.2　有效数字和最低位当量

单精度浮点数所能表示的最小绝对值约为10 ^-44 ，双精度浮点数所能表示的最小绝对值约为10 ^-323 。但这并不是说一个单精度或双精度浮点数在任何时候都可以准确表示出位于小数点后44位或323位的有效数字。单精度浮点数和双精度浮点数分别只有大约7位和16位十进制有效数字，其最低位有效数字所表示数值的大小取决于浮点数本身的绝对值。例如，对于一个大于10 ¹⁶ 的浮点数，其小数部分所有的数字都是不准确的。

C语言浮点数的表示精度有限，对于在其尾数的有效位数内所不能精确表示的数据会产生一定的舍入误差。当然，这些初始的舍入误差只发生在浮点数二进制尾数部分的最低位。在一次转换或计算中所产生的误差最大不超过最低有效位单位当量的1/2。对于float型的数，最低有效位单位当量只相当于尾数所能表示的最大值的2 ^-23 ，即约相对于数据大小的800万分之一。对于double型的数据，初始误差只相当于尾数所能表示的最大值的2 ^-52 ，这就更微不足道，完全可以满足一般数值计算精度的要求。

尽管对于一般的数值计算来说，C语言浮点数的精度已经足够了，但是对于复杂的计算，如果没有正确处理计算的精度，使得误差累积或者扩散，也会产生难以预期的结果。在浮点计算中经常会遇到的一个问题就是，计算所产生的结果有可能与预期不同：我们认为等价的计算会得出不同的结果，我们认为相同的数值却不相等，而这也常常是由于计算精度有限所造成的误差引起的。因此在描述浮点数计算时，必须深入理解浮点计算的特殊性，并根据程序的需要进行必要的处理。忽视了程序中浮点计算的这些特殊性，对于一般的程序运行看起来影响似乎不大，但是对于关键性的应用来说，则有可能产生不可预期的，甚至是灾难性的后果。一个被经常引用的例子是，在1991年海湾战争期间，美军的一枚爱国者导弹由于控制计算机中一个数据计算误差的累积，导致发射失败，击中了美军自己的兵营，击毙了近30名美军士兵。

从理论上讲，C语言浮点数对实数的近似表示使得C语言中的浮点计算不具有实数计算的一些基本性质。在运算结果不出现无穷大的情况下，浮点数的加法和乘法运算满足单调性的要求，即当a>b且c>=0时，浮点运算可以保证a+c>b+c和a * c>b * c。浮点数的加法和乘法运算也满足交换性，即a+b等于b+a，以及a * b等于b * a。但是，C语言的浮点加法运算不具有可结合性，也就是说，在C语言的浮点计算中，不能保证(a+b)+c等于a+(b+c)。例如，(1.23+1e16)-1e16就不等于1.23+(1e16-1e16)。前者的结果为2.0，而后者的结果为1.23。类似地，(1.23+1e17)-1e17的结果为0.0，而1.23+(1e17-1e17)的结果仍然为1.23。这是因为表达式的不同的结合方式改变了运算的顺序，因而也有可能改变舍入误差，并最终影响表达式的求值。当首先将1.23与1e16或1e17相加时，由于1.23接近或小于1e16或1e17中最低有效位的数值当量，因此计算中产生了较大的误差，或者干脆就导致了1.23被忽略。而当首先计算1e16或1e17相减时，因为减数与被减数相等，结果等于0，所以1.23再与这个结果相加时，可以得到完整的表达，不会在计算过程中损失精度，因此计算的结果是准确的。

除了计算顺序的不同可能会产生不同的计算结果外，在算术上等价，但表示方法不同的计算表达式也有可能产生不同的结果。下面就是一个这样的例子。

【例2-3】表示方法不同的等价浮点表达式产生不相等的计算结果

t = 1.234;
x = t * 2.01;
y = t + t + t / 100.0;
if (x == y)
    printf("X (%.10f) == Y (%.10f)\n", x, y);
else
    printf("X (%.10f) != Y (%.10f)\n", x, y);

这段代码的执行结果如下：

X (2.4803400000) != Y (2.480340000)

从表面上看，这一结果不仅出乎意料，而且难以理解：两个等价的运算过程应该产生两个相等的结果，从显示出来的数值看，结果也明明是一样的，但比较结果却说它们不相等！这一示例典型地说明了C语言中浮点运算的特点：等价的表达式所得出的结果不一定相等。之所以会产生这样的结果，是由于在不同的运算过程中，不同的计算步骤所带来的舍入过程是不一样的，舍入误差的累积也可能是不同的，因而导致了最终结果的不相等。把这两个数放大一些，就可以更清楚地看到这一点。为此可以改变t的初始赋值：

t = 1.234e16;

变量t的值被放大了10 ¹⁶ 倍，因此x和y的值也相应地放大了10 ¹⁶ 倍。这时程序的执行结果如下：

X（24803399999999996.0000000000) != Y (24803400000000000.0000000000)

可以清楚地看出，x和y之间的误差为4，也就是2 ² ，因此这两个数不相等就是很显然的了。因为准确表达x和y这两个整数需要55位二进制位，超过了C语言中double类型的53位尾数长度，所以x和y的最低有效位的单位当量是2 ² 。这说明这两个等价但运算方式不同的表达式之间的差异产生在最低有效位上：一个结果的最低有效位是0，而另一个是1。因为最低有效位的单位当量是整数4，所以在这一位上的差异在打印输出时可以显示出来。而当数值的绝对值较小时，其最低有效位的单位当量也很小，在精度有限的打印转换过程中，这一差异无法表现出来，因此就产生了前面看似难以理解的结果。下面我们再看一个由有效数字引起程序错误的例子。

【例2-4】以float类型变量作为循环控制变量引起无限循环

float a, n;
......
for (a = 0.0; a < (float) n; a += 1.0) {
    ......
}

上面的代码是一个简单的for循环语句，如果在循环体中没有break和goto之类的跳转语句，我们预期程序在执行完n次循环之后会结束对for语句的执行。但是实际上，这一预期只有在n是一个比较小的数值时才是正确的。当n大于16 777 216，也就是2 ²⁴ 时，这一语句会无限循环地执行下去。检查一下变量a的值就会发现，尽管float类型的数据可以表示的最大值大于10 ³⁸ ，但是变量a的值在增加到16 777 216.0之后就不再增加了，因此循环的条件永远满足，循环也就无限地进行下去。变量a的值之所以停止在16 777 216，是因为a的尾数最低位当量随a中值的增加而增加。当a的值增加到了2 ²⁴ 之后，对于尾数只有23位的float类型来说，在表示这一数值时其尾数的最低位当量等于2，这时a+=1.0对a的增加值不大于其尾数最低位当量的1/2，因此不会改变数据原有的值，a的值也就不再增加了。

浮点数据尾数有效数字较少还有可能带来计算误差和其他不易察觉的问题。例如，级数 是一个发散的级数，当n趋于∞时，它的值也趋于∞。但是当编写程序来计算这个无穷级数时，计算结果却停留在一个固定的极限上，而且这个极限取决于计算时所使用的浮点数据类型。当使用double类型进行计算时，这个极限大约不到35，而当使用float类型时，这个极限约等于15.403 683。导致计算出现这种结果的原因也很简单：因为C语言浮点数的精度有限，所以其所能表示的最小值也是有限的，无穷级数中小于这一限度的项就一律被表示成0.0，这样，对无穷级数的求和就变成了对有穷级数的计算。因为float类型的精度低于double类型，其所能表示的最小值大于double类型，所以在使用float类型进行计算时，级数中的有效项要远远少于使用double类型时，所以求和的过程更早地收敛，结果也自然要小于使用double类型了。不仅如此，计算的顺序也同样会影响计算的结果。例如，使用float类型时的极限15.403 683是按照对i的升序计算1/i得到的，而当对i按降序计算时，如果从i=2 ²⁷ 开始计算，则这个极限约等于18.807 918 5。这是因为在按照升序计算1/i时，首先被累加的是数值较大的项。这些项的累加结果迅速增加，使得数据的最低位当量也随之增大，因此其后面数值较小的项就会较早地被抛弃。而当采用逆向顺序计算时，累加从数值较小的项开始。因为这时数据的累加结果较小，float类型有较高的精度来表示这些数值较小的项的值，所以这些数值较小的项的值，或者至少它们的一部分会被计入最终的累加结果中。

浮点计算的结果除了受到数值表示精度和计算方式的影响外，也有可能受到硬件及编译系统的影响。例如，下面的代码在IA32平台上，由不同的编译系统产生的可执行文件的执行结果是不同的：

double t = 1.0, x = 10.0;
printf("%d\n", 1.0 / 10.0 == t / x);

当使用Linux/GCC 3.2.2时，程序打印出0；而当使用Windows/VC++ 6.0时，程序打印出1。

从上面这些例子中可以看出，因为C语言中浮点数是一种对实数的近似，所以其计算过程具有一些与普通算术运算不同的性质，会产生各种与普通算术运算不同的结果。因此，在处理浮点数时应当充分考虑到计算精度对计算结果的影响。当这种影响大到足以改变结果的正确性时，就需要采取必要的措施。可能采取的措施包括采用更高精度的数据类型，调整计算的次序，修改算法等，以保证计算结果的误差在可以接受的范围之内。

2.2.3　浮点数的比较

从上面的例子中可以看出，在C语言中不同的计算顺序和不同的等价表达式可能产生不同的结果。尽管这些结果之间的误差可能很小，但是有时仍有可能对后续的计算产生较大的影响，特别是对于数值之间是否相等的逻辑判断。因此，对于浮点数的比较既不同于一般算术中实数之间的比较，也不同于程序中整数之间的比较。这是因为，无论是算术中实数之间的比较，还是程序中整数之间的比较，都是可以精确表达的数值之间的比较，而浮点数之间的比较是不一定能够精确表达的数值之间的比较。在进行浮点数的比较时，必须要考虑到数值表示的精确程度，以及在计算过程中可能引入的舍入误差及其积累，并对这些误差进行妥善的处理。这就要求我们根据计算任务的精度要求，在容许的误差范围之内，适当地选择比较的方法，控制比较的过程。下面我们再看一个浮点数比较的例子。

【例2-5】循环语句中使用浮点数比较 按照0.1的步长打印出从0.0到0.9的数值。

下面这段代码使用double类型的i作为循环控制变量，将其初始化为0.0，利用打印语句输出i的值，在每次输出后将i的值增加0.1。只要i不等于1.0，循环就一直进行下去。10个0.1相加等于1.0，因此我们似乎有理由期望程序在输出10个数之后结束循环。

double i;
for (i = 0.0; i != 1.0; i += 0.1)
    printf("%.2f\n", i);

但实际情况却不是这样。代码的执行结果是，大量的数据在屏幕上不断地滚过，程序进入无限循环而无法停止。如果有办法让程序输出的速度慢一些的话，我们可以看到，程序的输出从0.00开始，然后是0.10、0.20、……，一路增长上去，在输出了0.90之后，程序并没有停住，而是继续输出了1.00。此后，程序继续不断地按每次0.1的增幅输出新的数值。由此可知，循环的持续条件是一直被满足的。也就是说，即使程序在屏幕上输出了1.00，i的值依然不等于双长浮点数1.0。造成这一结果的原因是，由于误差的累积，10个0.1相加并不严格等于1.0。这两个值之间有着微小的差异，因此条件i!= 1.0依然成立。在此之后，i的值继续累加，就更不可能等于1.0了，因此循环的条件永远满足，程序就进入了无限循环。

既然==和!=不能满足对浮点数的比较的要求，可能的解决方法是使用<和>来进行浮点数的比较。在上面的例子中，我们可以把for语句中的循环条件改为i<1.0：

for (i = 0.0; i < 1.0; i += 0.1)
    printf("%.2f\n", i);

对于浮点数的比较来说，使用<和>确实要比使用==和!=好得多。但是这仍然没有从根本上解决浮点数中微小误差对比较结果产生影响的问题。例如，上面这段代码的运行结果是：

0.00
0.10
......
0.90
1.00

这说明变量i每次递增0.1所产生的误差是舍弃误差，因此10个0.1相加的结果依然小于1.0，只有在加上第11个0.1后，循环条件才不再满足。但是，如果在这个例子中使用其他一些数，例如0.5、0.25这样能够使用浮点数精确表达的数值，或者使用能够产生正向累计误差的数值，这段代码的运行结果就是正确的了。这样，这类代码的执行结果是否正确，在很大程度上要取决于参与计算的具体数据，因此使得程序具有很大的不确定性，而这种情况在程序设计中是绝对应该避免的。除了浮点数之间的直接比较外，涉及返回浮点数值的函数的比较也存在相同的问题。

在进行浮点数的比较时一个最重要的原则是，应尽量避免两个浮点数之间的直接比较，特别是在判断两个浮点数是否相等或不相等时。也就是说，在判断两个浮点数是否相等时，不应直接使用运算符==；在判断两个浮点数是否不相等时，也不应直接使用运算符!=。当判断两个浮点数a和b是否相等时，正确的方法应当是使用语句 if(a-b<ESP && a-b>-ESP) 或者 if(fabs(a-b)<ESP) 。类似地，对于判断一个浮点数是否等于0，也不应当使用 if(x==0.0) ，而应当使用 if(x<ESP && x>-ESP) ，或者 if(fabs(x)<ESP) 。在上面各个语句中，ESP是一个在程序中自行定义的正浮点数常量，用作比较误差的控制值，函数fabs()是一个标准库函数，其功能是返回浮点数的绝对值。ESP的值取决于参与比较的数值的绝对值的大小和对计算精度的要求。ESP应该大于该数在计算中正常的舍入误差，因此任何由于浮点数的基本计算误差引起的数值变化都不会影响比较的结果；同时，ESP也应该足够小，以便任何由浮点数的基本计算误差之外的其他因素所产生的数值变化都不会被掩盖。需要注意的是，在一个程序里，不同的浮点数计算过程中参与计算的数值范围和对计算精度的要求不一定相同，因此可能需要不同的比较误差控制常数。

为了便于计算ESP的值，C编译器在头文件<float.h>中定义了两个常数宏，一个是FLT_EPSILON，另一个是DBL_EPSILON，分别表示对于float类型和double类型的浮点数，1.0+x不等于1.0的最小值。在IA32结构下的gcc中，这两个常数定义如下：

#define FLT_EPSILON 1.19209290e-07F
#define DBL_EPSILON 2.2204460492503131e-16

它们分别是当数值为1.0时相应类型的浮点数尾数的最低有效位所代表的数值。参照这两个数值，结合程序中进行比较的数值范围和对计算精度的要求，就可以较好地选择比较误差控制常数ESP的值了。类似地，在需要准确地进行大于和小于的比较时，也应留出足够的误差余量，以保证比较结果的正确性。

一般地说，在循环控制语句中，特别是当对循环次数进行控制时，应当尽量避免使用浮点数作为循环控制变量和条件，而应该使用可以精确表达的int类型数据。例如，在【例2-5】中，我们可以把代码改成下面的形式：

int i;
for (i = 0; i < 10; i++)
       printf("%.2f\n", i/10.0);

这样就把不精确的浮点数比较转换为精确的整型数比较，从根本上避免了浮点数比较可能引起的错误。

2.2.4　浮点数值计算中的上溢和下溢

在数值计算中，数据的类型和计算方法往往取决于所要计算的数值的范围以及对计算精度的要求。因为任何数据类型的精度都是有限的，所以对于有效数字较多的数据，就可能会损失精度。这种情况往往出现在绝对值很大的数据和绝对值很小的数据同时参加运算的时候。虽然对于大量的常规计算，double类型的浮点数的数值表示区间和精度一般都可以满足任务对计算精度的要求，其所产生的误差都在可以允许的范围之内。但是仍然有不少实际的计算任务，在其计算过程中，运算的中间结果会超越double类型的数值表示区间或精度的允许范围，造成上溢或下溢，导致最终计算结果的错误。下面我们看一个例子。

【例2-6】泊松（Poisson）分布的概率计算 泊松分布是一种常用的离散型概率分布，数学期望为 m 的泊松分布的分布函数定义如下：

P ( m , k )= m ^k ×e ^{- m} / k ! ( k =0, 1, 2, 3, …)

设计一个函数，对于给定的 m 和 k (0< m <2000, 0≤ k ≤2500)，计算其概率。

因为 m 是实数， k 是整数，概率是一个介于0和1之间的实数，所以函数的原型可以定义如下：

double poisson(double m, int k);

根据概率的分布函数，可以直接写出函数的定义如下：

#include <math.h>
double factorial(int n)
{
    int i;
    double f = 1.0;

    for (i = 1; i <= n; i++)
            f *= (double) i;
    return f;
}

double p_1(double m, int k)
{
    //  P(m, k) = m ^ k / k! * e ^ -m
    return pow(m, k) * exp(-m) / factorial(k);
}

在这段代码中，pow()和exp()是在<math.h>中说明的标准数值计算函数。pow(m, k)返回 m ^k ，exp(-m)返回e ^{- m} 。这段代码是对服从泊松分布随机概率的直接描述，在表达形式上是正确的，对于较小的参数也可以进行计算，但是却不能满足对计算范围的要求。这是因为尽管函数的计算结果是一个介于0和1之间的实数，位于double类型所能够正确表示的数值区间，但是计算过程的中间结果却可能产生溢出。无论是 m ^k 还是 k !，远在 m 和 k 达到2000之前就会由于数值过大而产生向上的溢出。例如，尽管在函数factorial()的计算过程中采用了double类型，函数也只能计算到170!。因此等不到 k 等于2500时，函数factorial()就只能将结果表示成无穷大（INF），使得后续的计算无法进行。同时，因为double类型的指数部分只有11位，只能表示到2 ^±1024 ，所以e ^{- m} 远在 m 达到2000之前就会由于数值过小而变成0.0，使得所有后续计算的结果都等于0。即使 m 和 k 的值小于上述两个极端的限制， m ^k 仍有可能在计算中产生溢出。例如，当 m 等于100， k 等于160时， m ^k =100 ¹⁶⁰ =10 ³²⁰ ，已经超过了double所能表示的10 ³⁰⁸ 。为了能够正确地在给定范围内计算泊松分布的概率，需要对计算过程的描述进行修改，在符合计算公式要求的条件下调整计算的顺序，以避免在计算过程中产生中间结果的上溢和下溢。

观察泊松分布的定义可以看出，尽管公式中的每一项在参数值较大时都可能引起中间结果的上溢或下溢，但是它们的增长却是互相抵消的，因此最终的结果是一个介于0和1之间的实数，完全在double类型的表示范围之内。因此我们可以把公式中对各个项的独立计算改为各项之间交叉计算，使得各项的增加互相抵消，把中间结果控制在适当的范围之内。因此可以首先把泊松分布公式改写为 P ( m , k )= m ^k /( k ! * e ^m )，并据此写出如下的代码：

double p_2(double m, int k)
{
    double x = 1.0;
    int i;

    if (m > k) {
        for (i = 0; i < k; i++)
            x *= m / M_E / (k - i);
        for (; i < m; i++)
            x /= M_E;
        return x;
    }
    for (i = 0; i < m; i++)
        x *= m / M_E / (k - i);
    for (; i < k; i++)
        x *= (double) m / (k - i);
    return x;
}

这段代码中的M_E是定义在<math. h>中的符号常量，表示e的近似值。代码首先计算 m ⁿ /( n ! * e ⁿ )，其中 n 是 m 和 k 中的最小值。然后再根据 m 是否大于 k 来决定如何计算公式 m ^k /( k ! * e ^m )与 m ⁿ /( n ! * e ⁿ )的差值部分。采用这种方法，把对 m ^k 、 k !以及e ^m 的计算分散在循环的每一步中，使得分布在分子和分母上的乘积互相抵消，就可以避免直接计算这些函数所可能引起的中间结果溢出问题。但是这样的方式在 m 和 k 都比较大或者两者的差值比较大时，在计算 m ⁿ /( n ! * e ⁿ )的过程中仍然有可能产生上溢或者下溢。例如，当 m 等于2000而 k 等于1000时，函数p_2()的if语句中的第一个for语句执行完毕后，x中的值是实际结果的e ¹⁰⁰⁰ 倍，已经远远超出了double类型的表示范围。即使计算的中间结果没有产生溢出，由于其接近于double类型表示范围的边界，也会引起较大的计算误差，使得计算结果完全错误。例如，当 m 等于2000而 k 等于2100时，正确的结果应为7.442234e-04，而使用上面的代码计算得出的结果是7.322822e+17。为避免计算中的错误，可以进一步调整代码，以便把中间结果控制在double类型数据的表示范围之中。代码调整的方法有多种，下面是一种比较简单的方法：

double poisson(double m, int k)
{
    double x = 1.0;
    int em, mk, nk;

    em = mk = nk = 0;
    while (em < m || mk < k || nk < k) {
        while (mk < k && x < MAX_INTER) {
            x *= m;
            mk++;
        }
        if (x < MIN_INTER)
            return 0.0;
        while (nk < k && x > MIN_INTER) {
            nk++;
            x /= nk;
        }
        while (em < m && x > MIN_INTER) {
            x /= M_E;
            em++;
        }
    }
    return x;
}

在这段代码中，符号常量MAX_INTER和MIN_INTER是中间结果的控制值，可以分别定义为远离double类型表示范围极限的值，如1e280和1e-280。经过这样调整的代码，可以将中间结果的值控制在适当的范围，减少计算过程中的误差。代码中的if语句是为了避免当最终结果小于MIN_INTER时函数进入无限循环。除了上述这种方法之外，更简便的改进方法是使用对数计算来代替对公式的直接计算。使用对数可以把乘除运算转换为加减运算，把乘方运算转换为乘法运算，不但可以确保中间结果不会产生上溢或下溢，而且可以简化函数的代码。我们把这一改进作为练习留给读者。从这个例子可以看出，对数值计算的描述，即使使用浮点数，也不总是一个简单的把计算公式直接转换成C语言的过程。在对计算过程进行描述时，需要对计算过程的中间结果进行分析，并根据分析的结果对计算过程的描述加以调整，以保证在各种条件下计算过程中所产生的中间结果都处在合理的范围之内，避免最终结果的错误。 sg8ESelpCGIy+jbASP2AGisi3xzvDkWVityd0K4vvPGpbdL/ZrJ2LJF1LETcS0y/