剖析Intel IA32架构下C语言及CPU浮点数机制 zt

wanglej

剖析Intel IA32架构下C语言及CPU浮点数机制      作者：谢煜波 文章来源：本站原创 点击数： 更新时间：2004-9-10

（转载请注名原作者及出处）
pdf格式下载：http://www.binghua.com/Soft/Class2/Class5/200409/63.html

剖析Intel IA32架构下C语言及CPU浮点数机制

Version 0.01

哈尔滨工业大学 谢煜波

（email: xieyubo@126.com 网址:http://purec.binghua.com）

（QQ:13916830 哈工大紫丁香BBSID:iamxiaohan）



前言

    这两天翻看一本C语言书的时候，发现上面有一段这样写到

例：将同一实型数分别赋值给单精度实型和双精度实型，然后打印输出。

       #include <stdio.h>

       main()

       {

              float a;

              double b;

              a = 123456.789e4;

              b = 123456.789e4;

              printf(“%f\n%f\n”,a,b);

       }

运行结果如下：

       1234567936.000000

       1234567890.000000

    为什么同一个实型数据赋值给float型变量和double型变量之后，输出的结果会有所不同呢？这是因为将一个实型常量赋值给float型变量与赋值给double型变量，它们所接受的有效数字位是不同的。

    这一段的说法是正确的，但实在是太模糊了！为什么一个输出的结果会比原来的大？为什么不是比原来的小？这之间到底有没有什么内存的根本性原因还是随机发生的？为什么会出现这样的情况？上面都没有对此进行解释。上面的解释是一种最普通的解释，甚至说它只是说出了现象，而并没有很深刻的解释原因，这不免让人读后觉得非常不过瘾！

    书中还有下面一段：

（1）两个整数相除的结果仍为整数，舍去小数部份的值。例如，6/4与6.0/4运算的结果值是不同的，6/4的值为整数1，而6.0/4的值为实型数1.5。这是因为当其中一个操作数为实数时，则整数与实数运算的结果为double型。

    非常遗憾的说，“整数与实数运算的结果为double型”，这样的表述是不精确的，不论从实际程序的反汇编结果，还是从对CPU硬件结构的分析，这样的说法都非常值得推敲。然而在很多C语言的教程上我们却总是常常看见这样的语句：“所有涉及实数的运算都会先转换成double，然后再运算”。然而实际又是否是这样的呢？

    关于浮点数运算这一部份，绝大多数的C教材没有过多的涉及，这也使得我们在使用C语言的时候，会产生很多疑问。

先来看看下面一段程序：

/* -------------a.c------------------ */

#include <stdio.h>

double f(int x)

{

        return 1.0 / x ;

}



void main()

{

        double a , b;

        int i ;

        a = f(10) ;

        b = f(10) ;

        i = a == b ;

        printf( "%d\n" , i ) ;

}

这段程序使用 gcc –O2 a.c 编译后，运行它的输出结果是 0，也就是说a不等于b，为什么？

再看看下面一段，几乎同上面一模一样的程序：

/*---------------- b.c ----------------------*/

#include <stdio.h>

double f(int x)

{

        return 1.0 / x ;

}



void main()

{

        double a , b , c;

        int i ;

        a = f(10) ;

        b = f(10) ;

        c = f(10) ;

        i = a == b ;

        printf( "%d\n" , i ) ;

}

同样使用 gcc –O2 b.c 编译，而这段程序输出的结果却是1，也就是说a等于b，为什么？

国内几乎没有一本C语言书（至少我还没看见），解释了这个问题，在C语言对浮点数的处理方面，国内的C语言书几乎都是浅尝即止，蜻蜓点水，而国外的有些书对此就有很详尽的描述，上面的例子就是来源于国外的一本书《Computer Systems A Programmer’s Perspective》（本文参考文献2，以下简称《CSAPP》），这本书对C语言及CPU处理浮点数描写的非常细致深入，国内很多书籍明显不足的地方，就在于对于某些细节我们是乎并没有某种深入的精神，没有一定要弄个水落石出的气度，这也注定了我们很少出版一些Bible级的著作。一本书如果值得长期保留，能成为Bible，那么我认为它必须把某一细节描述的非常清楚，以至于在读了此书之后，再也不需要阅读其它的书籍，就能对此细节了如指掌。

    《CSAPP》这本书的确非常经典，遗憾的是此书好像目前还没有电子版，因此我打算以此书为基础（一些例子及描述就来自此书），再加上自己看过的一些其它资料，以及自己对此问题的理解与分析，详细谈一下C语言及Intel CPU对浮点数的处理，以期望在此方面，能对不清楚这部分内容的学弟学妹们有些许帮助。

    要无障碍的阅读此文，你需要对C语言及汇编有所了解，本文的所有实验，均基于Linux完成，硬件基于Intel IA32 CPU，因此，如果你想从此文中了解更多，你最好能熟练使Linux下的gcc及objdump命令行工具（非常遗憾的是，现在少有C语言教材会对此进行讲述），另外，你还需要对堆栈操作有所了解，这在任何一部讲解数据结构的书上都会提到。

    由于自身知识及能力有限，如果书中有描述不当的地方或错误，请你与我联系，我也会在哈工大纯C论坛上（http://purec.binghua.com）对� ... �跟踪及反馈。



一、Intel CPU浮点运算单元的逻辑结构

    在很久以前，由于CPU工艺的限制，无法在一个单一芯片内集成一个高性能的浮点运算器，因此，Intel还专门开发了所谓的协处理器配合主处理器完成高性能的浮点运算，比如80386的协处理器就是80387，后来由于集成电路工艺的发展，人们已经能够将在一个芯片内集成更多的逻辑功能单元，因此，在80486DX的时候，Intel就在80486DX这个芯片内集成了很强大的浮点处理单元。下面，我们就来看看，被集成到主处理器内部之后，这个浮点处理单元的逻辑结构，这是理解Intel CPU对浮点数处理机制的前提条件。


（图1 Intel CPU 浮点处理单元逻辑结构图）

    上图就是Intel IA32 架构CPU浮点处理单元的逻辑结构图，从图中我们可以看出它总共有8个数据寄存器，每个80位（10B）；一个控制寄存器（Control Register），一个状态寄存器（Status Register），一个标志寄存器（Tag Register），每个16位（2B）；还有一个最近一次指令指针（Last Instruction Pointer），及一个最近一次操作数指针（Last Operand Pointer），每个48位（6B）；以及一个操作码寄存器（Opcode Register）。

    状态寄存器用处与常见的主CPU的程序状态字差不多，用来标记运算是否溢出，是否产生错误等，最主要的一点是它还记录了8个数据寄存器的栈顶位置（这点在下面将会有详细描述）。

    控制寄存器中最重要的就是它指定了这个浮点处理单元的舍入的方式（后面将会对此详细描述）及精度（24位，53位，64位）。Intel CPU浮点处理器的默认精度是64位，也称为Double Extended Precision（中文也许会译为：双扩展精度，但这种专有名词，不译更好，译了反而感觉更别扭）。而24位，与53位的精度，是为了支持IEEE所定义的浮点标准（IEEE 754标准），也就是C语言中的float与double。

    标志寄存器指出了8个寄存器中每个寄存器的状态，比如它们是否为空，是否可用，是否为零，是否是特殊值（比如NaN：Not a Number）等。

    最后一次指令指针寄存器与最后一次数据指针寄存器用来存放最后一条浮点指令（非控制用指令）及所用到的操作数在内存中的位置。由于包括16位的段选择符及32位的内存偏移地址，因此，这两个寄存器都是48位（这涉及到Intel IA32架构下的内存地址访问方法，如果对此不清楚的，可以不用太在意，只需知道它们指明了一个内存地址就行，如果很想弄清楚可以参考看本文的参考文献1）。

    操作码寄存器记录了最后一条浮点指令（非控制用指令）的操作码，这很简单，没什么可多说的。

    下面我们将详细描述一下，浮点处理单元中的那8个数据寄存器，它们与我们通常用的主cpu中的通用寄存器，比如eax，ebx，ecx，edx等相比有很大的不同，它们对于我们理解Intel CPU浮点处理机制非常关键！



二、Intel CPU浮点运算单元浮点数据寄存器的组织

    Intel CPU浮点运算单元中浮点数据寄存器总共有8个，它们都是80位，即10字节的寄存器，对于每个字节所带表的含义我将在后面描述浮点数格式的时候详细介绍，这里详细介绍的将是这8个寄存器是怎么组织以及怎么使用的。

    Intel CPU把这8个浮点寄存器组织成一个堆栈，并使用了状态寄存器中的一些标志位标志了这个栈的栈顶的位置，我们把这个栈顶记为st(0)，紧接着栈顶的下一个元素是st(1)，再下一个是st(2)，以此类推。由于栈的大小是8，因此，当栈被装满的时候，就可以访问的元素为st(0)~st(7)，如下图所示：


（图2 装入不同数据时浮点寄存器中栈顶的位置）

    由上图可以很明显的看出浮点寄存器是怎样被组织及使用的。需要注意的是，我们并不能通过指令直接使用R0~R7，而只能使用st(0)~st(7)，这在下边描述浮点运算指令的时候，会有详细描述。

    也许会有朋友对上图产生疑问，当已经放入8个数后，也即当st(0）处于R0的时候，再向里面放入一个数会产生什么情况呢？

当已经有8个数存入浮点寄存器中后，再向里面放入数据，这会根据控制寄存器中的相应的屏蔽位是否设置进行处理。如果没有设置相应的屏蔽位，就会产生异常，就像产生一个中断似的，通过操作系统进行处理，如果设置了相应的屏蔽位，则CPU会简单的用一个不确定的值替换原来的数值。如下图所示：


（图3 装入数据大于八个数时，浮点寄存器状态）

    可见其实浮点寄存器相当于是被组织成了一个环形栈，当st(0)在R7位置的时候，如果还有数据装入，则st(0)会回到R0位置，但这个时候装入st(0)的却是一个不确定的值，这是因为CPU将这种超界看做是一种错误。

    那么上面的说法倒底对不对呢？别急，在下面描述了浮点运算之后，我将会用一段实验代码验证上面所述。



三、Intel CPU浮点运算指令对浮点寄存器的使用

在第二节中，我们指出Intel CPU将它8个浮点寄存器组织成为一个环形堆栈结构，并用st(0)指代栈顶，相应的，Intel CPU的相当一部份浮点运算指令也只对栈首的数据进行操作，并且大多数指令都存在两个版本，一个会弹栈，一个不会弹栈。比如下面的一条取数指令：

fsts 0x12345678

    这就是一个不会弹栈的指令，它只是将栈顶，即st(0)的数据存到内存地址为0x12345678的内存空间中，其中fsts最后的字母s表明这是对单精度数进行操作，也就是说它只会把st(0)中的四个字节存入以0x12345678开始的内存空间中。具体是那四个字节，这就涉及到从80位的st(0)到单精度（float）的一个转换，这将在下面介绍浮点数格式的小节中详细描述。

       上面的指令执行后，不会进行弹栈操作，即st(0)的值不会丢失，而下面就是同种指令的弹栈版本：

fstps 0x12345678

    这条指令的功能几乎于上面一条指令完全相同，唯一不同的地方就在于这是一个会引起弹栈操作的指令，其中fstps中的字母p指明了这一点。此条指令执行后，原来st(0)中的内容就丢失了，而原st(1)中的内容成为st(0)中的内容，这种堆栈的弹压栈操作我想对大家是再熟悉不过了，因此，这里将不再对其进行描述，不清楚的可以参看任一本讲数据结构的书。

本文主旨在于描述一下Intel CPU浮点数处理机制的基本原则，而并非浮点指令的资料，因此本文不再对众多的浮点指令进行描述，在下面的描述中，本文仅对所用到的指令进行简单的解释，如果你想完整了解浮点指令，可以参看本文的参考文献1。

下面，我们将用一个例子结束本节的讲述，这个例子将涉及上节及本节所讲述的内容，它验证了上面的描述是否正确。

请在Linux下输入下面的代码：

/* ---------------------------- test.c ------------------------------------ */

void f(int x[])

{

        int f[] = {1,2,3,4,5,6,7,8,9} ;

/*----------------------------- A 部分 ------------------------------*/

        __asm__( "fildl %0"::"m"(f[0]) ) ;

        __asm__( "fildl %0"::"m"(f[1]) ) ;

        __asm__( "fildl %0"::"m"(f[2]) ) ;

        __asm__( "fildl %0"::"m"(f[3]) ) ;

        __asm__( "fildl %0"::"m"(f[4]) ) ;

        __asm__( "fildl %0"::"m"(f[5]) ) ;

        __asm__( "fildl %0"::"m"(f[6]) ) ;

        __asm__( "fildl %0"::"m"(f[7]) ) ;

        // __asm__( "fildl %0"::"m"(f[8]) ) ;                                (*)

        // __asm__( "fst %st(3)" ) ;                                       (**)



/* ------------------------------ B部分 ---------------------------------*/

        __asm__( "fistpl %0"::"m"(x[0]) ) ;

        __asm__( "fistpl %0"::"m"(x[1]) ) ;

        __asm__( "fistpl %0"::"m"(x[2]) ) ;

        __asm__( "fistpl %0"::"m"(x[3]) ) ;

        __asm__( "fistpl %0"::"m"(x[4]) ) ;

        __asm__( "fistpl %0"::"m"(x[5]) ) ;

        __asm__( "fistpl %0"::"m"(x[6]) ) ;

        __asm__( "fistpl %0"::"m"(x[7]) ) ;

}



void main()

{

        int x[8] , j ;

        f(x) ;

        for( j = 0 ; j < 8 ; ++j )

                printf( "%d\n" , x[j] ) ;

}

    上面的代码通过内嵌汇编，在A部分把一个整数数组中的整数压入浮点寄存器中（fildl指令用于把整数压入浮点寄存器），而后又在B部分将浮点寄存器中的数取到另一个数组中（fistpl指令用于把栈顶数据存入指定内存单元中，指令中的字母p表明这是一个弹栈指令，每次都会弹栈）。程序中我们只压入了f[0]~f[7]的8个数据，而压入的数据顺序是1，2，3，4，5，6，7，8，因此，取出的顺序应当是8，7，6，5，4，3，2，1，在Linux下编译并运行我们会得到同样的结果。

    下面，我们将（*）语句处的注释符号“//”去掉，这个时候我们压入了f[0]~f[8]共9个数据，这将会引起超界，按照上面的描述，当发生这种情况的时候，st(0)会从R0的位置变到R7，并在其中存入一个不确定的值，那么实际情况是不是这样呢？同样请在Linux下编译并运行此程序，并将结果与图3进行比较。这里需要注意的时，我们总是按照st(0)，st(1)，……，st(7)的顺序取出数据的。

    最后，我们再将（**）语句处的注释符号去掉，“fst %st(3)”这条指令的作用是把st(0)中的内容存入st(3)，指令中并没有p字母，因此，这并不是一条会引起弹栈的指令。同样请在Linux下编译运行，并对照图3观察它的结果，以验证前文所述内容。



四、浮点数格式

    C语言及CPU所使用的浮点数格式均遵从IEEE 754标准，下面我们就对此详细的讨论一下。

    IEEE标准指出，一个数可以表示为：。其中S在V是正数时取0，在V是负数是取1。对应到程序中，这显然就是一个符号位，所以S占1位，而V及M的位数由数据类型来决定。如果是单精度型（float），那E占8位（e = 8），M占23位（m = 32），如果是双精度型（double），E占11位（e = 11），M占52位（m = 52），如下图所示：


（图4 IEEE 745 浮点数格式）

    这里需要注意的是，我们用的是S’、E’、M’而非S、E、M，这是因为它们之间存在着一种转换关系。请看下面的描述。

    S在任何时候都等于S’，但E与M就不一样了。当E’=0时，E，M = 0.M’，举个例子：0x00 50 00 00 这个数所表示的浮点数是多少呢？我们把这个数展开

0000 0000 0101 0000 0000 0000 0000

    最开头红色的1位是符号位S’，故S=S’=0，中间8位蓝色的是E’位，由于E’=0，所以按上面公式可算得E = -126，最后的是M’，所以M = 0. 101 0000 0000 0000 0000=0.625，所以这个小数应当是V==7.346840e-39，那么实际是不是这样的呢？我们同样通过下面一个程序进行一下验证：

union U{

        float f ;

        struct{

                unsigned char x1 ;

                unsigned char x2 ;

                unsigned char x3 ;

                unsigned char x4 ;

        } ;

} u ;



int main()

{

        u.x1 = 0x00 ;

        u.x2 = 0x00 ;

        u.x3 = 0x50 ;

        u.x4 = 0x00 ;



        printf( "%e\n " , u.f ) ;

        return 0 ;

}

    程序非常简单，你可以在Linux下编译并执行，以检查前文所述。

    上面谈到了当E’ = 0时的情况，那么当E’不为0，且不全为1（即E’）时又是什么情况呢？这个时候，E = E’，M = 1.M’，举个例子，0xbf 20 00 00这个数表示的浮点数是多少呢？同样把这个数展开：

1011 1111 0010 0000 0000 0000 0000 0000

    下面我们来按上述的要求计算：

S = S’ = 1，E = -1，M = 1.M = 1. 010 0000 0000 0000 0000 0000 = 1.25

所以

V = = -0.625

    同样，我们通过一个程序验证一下，不过这次我们把这个程序变一下，直接输入0.625，看它的输出字节是多少：

union U{

        float f ;

        struct{

                unsigned char x1 ;

                unsigned char x2 ;

                unsigned char x3 ;

                unsigned char x4 ;

        } ;

} u ;



int main()

{

        u.f = -0.625 ;



        printf( "%2x %2x %2x %2x\n",u.x4,u.x3,u.x2,u.x1) ;

        return 0 ;

}

    编译并运行之后，我们得到的输出是：bf  20  0  0 这与我们前面的分析完全一致。

这里还想提请大家注意一点的是，通过这个例子，我们很容易可以看见，IEEE格式中，负数并不是用补码来表示的！

    下面只剩下最后一种情况了：当E’全为1即E’时。这个时候，如果M’ = 0，此数表示无穷（inf，当S = 0时是正inf，当S = 1时是负inf），如果M’ 不是0，此数是NaN（Not a Number），比如，你让计算机计算-1的平方根，就会得到NaN，表明它无法用一个数表示。

    上面描述的是IEEE所定义的浮点格式，而在Intel CPU中使用的是一种扩展精度的浮点数形式，它的E有15位（e = 15），M有63（m = 63）位，加上1位的符号位，恰好等于80位，这同Intel CPU中浮点寄存器的长度是一样的。但是E，M位的确定方法还是同IEEE标准兼容，也即本节所描述的方法完全适用于Inte CPU的80位的扩展精度格式。



五、浮点数的舍入及所带来的程序设计上的问题

    由于存在多种格式的浮点数，因此，从一个高精度的格式转换到一个低精度的格式就会出现舍入，而由于舍入的存在，就会造成程序中出现很多非常有意思的错误，现在我们就来谈谈这个问题。这里将最终解决前言中提到的那个问题。

    首先谈谈从低精度格式转换到高精度格式。这并不会引起精度的丢失，因此，这样的转换是很安全的，随着表示数据的位数的增加，高精度格式可以完全把低精度格式的相应数据位复制过来，并不会丢失任何信息，然而从高精度向低精度进行转换，就会丢失信息，因为低精度格式的数据位数比高精度的要少，容纳不下高精度的所有信息。

    现在我们就来看一下各种精度格式可表示的最大正数或最小正数是多少。

首先来看一下最小正数，最小正数时E’=0，而M’=000…..1，故由前面的描述可知，对于单精度数（float），最小正数为，对于双精度数（double），最小正数数；最大正数时E’=111…10，而M’=111…1，故可知单精度数的最大正数为，双精度数（double）的最大正数为。很明显的可以看到双精度数所能表示的范围，远远大小单精度数所表示的范围。

当一个双精度数所表示的数比单精度所能表示的最大的数还要大时，CPU会让单精度数等于无穷，此时，单精度数的E’全为1，而M’=0（见上节所述）。

当一个双精度数所表示的数比单精度所能表示的最小数还要小时，CPU会让单精度数等于0，此时，单精度数的E’、M’全为0。

注意这样一个事实，M只有两种可表达形式，一种是M=1.M’，一种是M=0.M’，当M=0.M’的时候，指数E只能为，这个特点决定了任何一个非零的数只有一种确定的表示。这种确定性使计算机的处理变得非常方便的。

比如一个double型的浮点数，它在用float型来表示时只能有一种型式，这使我们能很快的确定出float型的数据表示，也一眼就能确定这个数是否造出了float表示数的范围。同样，对于计算机来说，这极大的方便了不同精度之间的转换，使这种转换有唯一而确定的形势。

在CPU的内部，所有的浮点数的在被浮点指令装入浮点寄存器的时候都会发生转换，从单精度、双精度、整数转换为80位的扩展精度，当从浮点寄存器存入内存的时候又会发生转换，从扩展精度转换为相应的精度格式，这种转换是由CPU硬件自动完成的，然而正是由于从扩展精度转换为低精度格式这一行为的存在，会让我们的程序出现某些很奇怪的问题。请看下面一段代码：

int main()

{

        float f = 3.25 + 1e10;

        f = f - 1e10 ;

        printf( "%f\n",f ) ;

        return 0 ;

}

    这段代码就是一个典型的精度丢失的例子，我来现在就来认真的分析一下它：

3.25 + 1e10这个数我们用二进制可以表示为：

1001 0101 0000 0010 1111 1001 0000 0000 11.01

由于M’只能是0.M’的形式或1.M’的形式，如果是0.M’的形式，则上式可以表示为：

0. 1001 0101 0000 0010 1111 1001 0000 0000 1101

    如果是1.M’的形式，则可以表示为：

1.001 0101 0000 0010 1111 1001 0000 0000 1101

    现在我们将它们转换为单精度数，按照IEEE单精度数的格式要求，如果是0.M’的形式，则指数只能是，而上面的E=34，故，上面的数只能采用1.M的形式，这个时候，M’=001 0101 0000 0010 1111 1001 0000 0000 1101，但由于IEEE指出，在单精度格式中，M’只能有23位，因此，最后的0000 0000 1101将会被截断，M’实际为：

M’=001 0101 0000 0010 1111 1001

也就是说，原数中的3.25被丢失了，因此，实际上f=3.25+1e10计算后f=1e10，所以，上面这段代码的输出结果是0。

这里需要提到一点上面所涉及到的截断方法，专业一点的术语称为舍入（Round）。IEEE总共定义了四种舍入规则，分别是“Round to Even”，“Round Toward Zero”，“Round Down”及“Round UP”，到底使用哪一种舍入规则可由程序员在浮点单元的控制寄存器中指定（参见前文所述），默认是使用“Round to Even”，下面我们就来看看这些规则。

先看“Round to Even”，假若有个数x.yyyy0….那么舍入为x.yyyy；如个有个数为x.yyyy1….且1后面的所有位不全为零，那么舍入为x.yyyy+0.0001；如果原数为x.yyyy100….00那么这个时候的舍入就需要分情况讨论了。如果此时最后一个y是1，那么舍入为x.yyyy+0.0001，如果此时最后一个y是0，那么舍入为x.yyyy。比如1.0110 100这个数舍入后为1.0110，而1.0111 100，这个数舍入后就为1.0111+0.0001=1.1000。

在来看看“Round Toward Zero”，这个很简单，它要求舍入后的数的绝对值不大于原数的绝对值。

“Round Down”要求，舍入后的数不大于原数。

“Round UP”要求，舍入后的数不小于原数。

上面关于精度损失的例子是个比较明显而比较简单的例子。但事情不总是这么明显的，下面我们就来解剖一下前言中提到的那个程序：

    先来回顾一下前言中提到的那个问题：

先来看看下面一段程序：

/* -------------a.c------------------ */

#include <stdio.h>

double f(int x)

{

        return 1.0 / x ;

}



void main()

{

        double a , b;

        int i ;

        a = f(10) ;

        b = f(10) ;

        i = a == b ;

        printf( "%d\n" , i ) ;

}

这段程序使用 gcc –O2 a.c 编译后，运行它的输出结果是 0，也就是说a不等于b，为什么？

再看看下面一段，几乎同上面一模一样的程序：

/*---------------- b.c ----------------------*/

#include <stdio.h>

double f(int x)

{

        return 1.0 / x ;

}



void main()

{

        double a , b , c;

        int i ;

        a = f(10) ;

        b = f(10) ;

        c = f(10) ;

        i = a == b ;

        printf( "%d\n" , i ) ;

}

同样使用 gcc –O2 b.c 编译，而这段程序输入的结果却是1，也就是说a等于b，为什么？

    我们现在将第一段代码，也即a.c反汇编，下面是反汇编的结果：

08048328 <f>:

8048328:       55                      push   %ebp

8048329:       89 e5                   mov    %esp,%ebp

804832b:       d9 e8                   fld1  

804832d:       da 75 08                fidivl 0x8(%ebp)//计算，并将结果存入st(0)中

8048330:       c9                      leave

8048331:       c3                      ret   

8048332:       89 f6                   mov    %esi,%esi

上面代码中我们只需要注意一点：fidivl指令表示用st(0)中的数除以某一内存地址中的数，并将结果存在st(0)中。

注意这样一个事实，现在的运算结果是存在st(0)中的，而且这是一个80位的值。下面我们来看看main中的代码：

08048334 <main>:

8048334:       55                      push   %ebp

8048335:       89 e5                   mov    %esp,%ebp

8048337:       83 ec 08                sub    $0x8,%esp

804833a:       83 e4 f0                and    $0xfffffff0,%esp

804833d:       83 ec 0c                sub    $0xc,%esp

8048340:       6a 0a                   push   $0xa

8048342:       e8 e1 ff ff ff          call   8048328 <f>

8048347:       dd 5d f8                fstpl  0xfffffff8(%ebp)

804834a:       c7 04 24 0a 00 00 00    movl   $0xa,(%esp,1)

8048351:       e8 d2 ff ff ff          call   8048328 <f>

8048356:       dd 45 f8                fldl   0xfffffff8(%ebp)

8048359:       58                      pop    %eax

804835a:       da e9                   fucompp

804835c:       df e0                   fnstsw %ax

804835e:       80 e4 45                and    $0x45,%ah

8048361:       80 fc 40                cmp    $0x40,%ah

8048364:       0f 94 c0                sete   %al

8048367:       5a                      pop    %edx

8048368:       0f b6 c0                movzbl %al,%eax

804836b:       50                      push   %eax

804836c:       68 d8 83 04 08          push   $0x80483d8

8048371:       e8 f2 fe ff ff          call   8048268 <_init+0x38>

8048376:       c9                      leave

8048377:       c3                      ret   

代码很长，但我们实际只需关心其中的一小部份，请看：

8048342:       e8 e1 ff ff ff          call   8048328 <f> //计算f(10)，这个时候的

                                                   //计算结果在st(0)中

8048347:       dd 5d f8                fstpl  0xfffffff8(%ebp) //把计算结果存回内存

//a中

804834a:       c7 04 24 0a 00 00 00    movl   $0xa,(%esp,1)

8048351:       e8 d2 ff ff ff          call   8048328 <f> // 计算f(10)，对应于b=f(10)

                                                    // 计算结果在st(0)中

8048356:       dd 45 f8                fldl   0xfffffff8(%ebp) // 直接载入a中的值

                                                     // 这时st(0)= a

                                                     // st（1）方才计算的b的值

8048359:       58                      pop    %eax

804835a:       da e9                   fucompp // 将 st(0)与st(1)进行比较

804835c:       df e0                   fnstsw %ax

    这里我们已经能够看出问题所在了！它先计算了a=f(10)，然后把这个结果存回到内存中了，由于0.1没办法用二进制精确表示，因此，从80位的扩展精度存到内存的64位的double中，产生了精度损失，这之后，计算b=f(10)，而这一值并没有被存回内存中，这个时候，gcc就直接将内存中的a值装入到st(0)中，与方才计算的b值进行比较，由于b值并没有被存回内存中，因此，b值并没有精度损失，而a值是损失了精度的，因此a与b不相等！

    下面我们来把第二段代码反汇编，这里我们只贴出我们最需要关心的那部份代码：

8048342:       e8 e1 ff ff ff          call   8048328 <f> // 计算a

8048347:       dd 5d f8                fstpl  0xfffffff8(%ebp) // 把a存回内存

                                                          // a产生精度损失

804834a:       c7 04 24 0a 00 00 00    movl   $0xa,(%esp,1)

8048351:       e8 d2 ff ff ff          call   8048328 <f> // 计算b

8048356:       dd 5d f0                fstpl  0xfffffff0(%ebp)// 把b存回内存

                                                         // b产生精度损失

8048359:       c7 04 24 0a 00 00 00    movl   $0xa,(%esp,1)

8048360:       e8 c3 ff ff ff          call   8048328 <f> // 计算c

8048365:       dd d8                   fstp   %st(0)

8048367:       dd 45 f8                fldl   0xfffffff8(%ebp) // 从内存中载入a

804836a:       dd 45 f0                fldl   0xfffffff0(%ebp) // 从内存中载入b

804836d:       d9 c9                   fxch   %st(1)

804836f:       58                      pop    %eax

8048370:       da e9                   fucompp // 比较a , b

8048372:       df e0                   fnstsw %ax

    从上面的代码看出，a与b在计算完成之后，都被gcc存回了内存，于是都产生了精度损失，因此，它们的值就是完全一样的了，于是此后再把它们调入浮点寄存器进行比较，得出的结果就是 a = b。这主要是因为程序中多了一个c = f(10) 的计算，它使gcc必须把先前计算的b值存回内存。

    我想，现在你应当对此问题有比较清楚的了解了吧。顺带说一句，gcc中提供了一个long double的关键字，对应于Intel CPU中的80位的扩展精度，它是10字节的。



六、总纲

    据说《九阴真经》在最后有一个总纲，对全书进行总结，这里，我也借鉴一下，对前面的行文进行一下总结。

    这篇文章比较完整的描述了Intel CPU对浮点数的基本的处理机制，希望能对大家理解CPU及在运用C语言进行浮点编程时产生一定有益的影响。在最后，我想对几个比较模糊的说法谈谈自己的看法。

    首先，浮点数不能比较是否相等。其实这个说法是不精确的，从本质上说，浮点数完全可以比较是否相等，CPU也提供了相应的指令，但由于存在舍入问题，因此，浮点数在比较是否相等的时候是不安全的，比如上面分析的那个程序，同样计算的是f(10)，但第一个程序在比较是否相等时是不成立的，而第二个却是成立的，这主要是由于舍入问题的存在，如果我们能够从本质上理解这个问题，那么我们完全可以放心大胆的在程序中比较两个浮点数是否相等。对于上面的这个程序我还想说明一点就是，我们必须在gcc打开优化的前提下才能得出上面的结果，这是因为gcc在非优化的时候每次计算完结果后，都会把结果存回内存。程序优化的最首要的目标就是不能影响程序的执行结果，但非常遗憾的是，gcc并没有做到这一点。

    其次，我们常常在书上看见：“所有的涉及浮点数的运算，都会转成double进行”，从上面的分析中我们也可以看出，这是不精确的。从硬件的角度说，Intel CPU在默认条件下使用的是扩展精度（gcc中的long double类型），因此，它们都是转换成扩展精度进行的。从C语言本身来说，C语言做为一个同硬件非常贴近的语言，他最终的结果就是产生出硬件可以识别的代码，而硬件到底怎么执行是由硬件决定的，C语言无法对此进行控制与干涉，相反C语言必须去适应硬件的规定。因此，不能说C语言本身会将所有浮点运算都转换成某种精度进行计算，到底使用哪种精度这是由cpu决定也是由cpu来完成的。

    第三，由于浮点运算后，结果在存入内存中会产生舍入，这很可能会带来误差，因此，我们应当尽量使用高精度的数据类型，比如用gcc的时候，我们尽量使用long double而非double、float，这会减少相当多的错误机会。不要认为使用它们会带来性能上的下降，它们最主要的弱点在于占用的空间比较大，但现在空间已经不是一个主要考虑的因素。



七、推荐阅读

    C语言是一个与硬件贴得很近的语言，要想真正精通C语言，你应当对硬件结构及指令系统有相当的了解，这你可以参看本文的参考文献1。

    如果你认为参考文献1写得太过繁复，你可以参看一下本文的参考文献2，另外，它还对在编程中可能出现的问题做了很多分析，非常深入细致，在某些细节上描述得非常清楚，具有Bible的品质。

    如果你对C语言不是很了解，你可以参考一下参考文献3，它对C语言有比较简单的描述，特别是书中有一章《C程序设计常见错误及解决方案》对初学者有很大好处。





参考文献：

1. 《IA-32 Intel Architecture Software Develop’s Manual》Vol.1，Vol.2，Intel Corp

2. 《Computer Systems A Programmer’s Perspective》（Randal E.Bryant，David R.O’Hallaron）电子工业出版社（影印版）

3. 《C语言大学实用教程》（苏小红，陈惠鹏，孙志岗）电子工业出版社