写在前面

在PA_2019fall中有一项任务是完成CPU中的浮点数运算，这也是我第一次认真的思考了一下真实的计算机中CPU是如何进行的浮点数运算

在写PA的过程中一头雾水，从迷茫，到困惑，到弄懂，到完成，中间经历了各种坎坷，又无奈于手上的资料仅仅只有一个Guide和i386，剩下不会的地方全靠百度

于是就诞生了这一篇博客，把其中的过程给大家讲的明明白白的！

预备知识

什么是浮点数？浮点数表示的是一个数字，其小数点所在的位置是不确定的，也就是浮动的，因此称之为浮点数

在IEEE 754标准中，单精度的浮点数由3部分组成

符号位，在首位，用1表示负数，0表示正数

阶码

部分，一共八位，采用移码形式，偏置常数为127

尾数部分，一共23位。

规格化于非规格化

规格化的数表示阶码部分不是全都为0，其23位的尾数实际上表示了24位，最高缺省位为1。而非规格化的数没有缺省位。

尾数加上最高位可能存在的缺省的1后构成的有效数字称为

实现之前的思考

排除特殊情况

Simple

根据表格我们可以先把边界条件排除掉，比如

加减

如何让两个浮点数相加？根据上面的知识我们知道两个浮点数的形式基本可以这样给出

，我们很容易就能想到两个浮点数是不能直接相加的，应该让他们同阶后再将尾数相加，这个操作我们叫对阶

对阶

对阶有两种方式，小阶往大阶对齐和大阶往小阶对齐，该如何选择呢？

不妨考虑一下这样的事情，在IEEE754表示的浮点数中，对于相同的阶数来说，尾数中越靠后面的1对最后结果的影响越小。

在对阶中，如果大阶向小阶看齐，那么实际上是将尾数左移，小阶向大阶看齐的话就是尾数右移。

如果我们选择大阶向小阶看齐，那么我们很容易就将尾数中靠前的1给左移没了，因为尾数一共只有23位，这样对计算造成的误差相当大，所以我们选择小阶向大阶看齐

注：考虑如何减小误差是浮点数运算中非常重要的一环

对阶中的特殊情况

在处理非规格化浮点数的时候，它们的阶码全都是0，但是尾数并非是0，其表示的真实值为

，不能将阶码理解成-127，在后面会讲为什么会有这样的结果。

shift的计算

在此，我们已经可以给出右移的位数

了，在这里我们假设

shift = (fb.exponent==0 ? fb.exponent+1:fb.exponent) - (fa.exponent==0? fa.exponent+1:fa.exponent)

最后我们得到的临时的阶码就是

流程图

先让我们把目前的东西总结一下，后面结合实际栗子来讲也会更加清晰

流程图

在这个图中，我们可以看到在对阶之前先判断了

的值是否为0，如果不是那就开始对阶

在对阶的过程中，小阶慢慢变大，这会导致尾数右移(前面已经提到过)，在尾数右移的过程中，我们有可能把尾数移为23个0，如果是这样，我们简单的判断为是阶数小的数实在是太小了，就像10亿+1，那个1加不加我们可能都没有直观的体验

再接下来，我们对尾数进行相加，这个过程中要考虑符号位，如果是做减法，那就对减数的尾数位采用广义上的 取反加一

相加完后，如果尾数变成0，那结果就是0，如果不是0，那就考虑一下尾数有没有溢出(尾数溢出指相加后最高位缺省位达到了2)。如果溢出了，那我们就右移一次，然后判断一下阶数的情况即可

要是没有溢出，那当然万事大吉，但是我们接下来还需要将我们的浮点数规格化一下，意思是规格化的浮点数为

，但是我们计算出来的浮点数可能是

，那么我们需要小数点移位，降低阶数，将最高缺省设置成1

当然，规格化的过程中阶数也有可能等于0，这里需要一次特判。

小小的总结

到这里浮点数的加减基本上已经介绍完了，对于乘除法就很简单了，只需要将指数相加，尾数相乘就行了，当然，其中涉及到的一些溢出情况都需要单独判断，在这里不多做讨论

附加位与舍位

一个贴近生活的栗子

加入你是亿万富翁，你的银行卡中有1000亿元，你每天稳定收入1000元，但是由于银行的计算机用的是上面的实现方法，这1000块钱在放进你的卡中的过程中需要进行对阶，对着对着尾数变成了0，相当于你存进去0块

那么该怎么办呢？当然你可以设置一些if语句，让1000加多点，加到足够影响1000亿的时候再一起加上去，这是目前计算机的一个研究方向之一。

附加位

根据IEEE754的规定，所有浮点数运算的中间结果右边都必须至少保留两位的附加位，依次是保护位(guard,G)和舍入位(round,R)，为了进一步提高精度，舍入位的右侧还有一个粘位(sticky,S)，只要舍入位右边有任何非0的数字，粘位就是1，否则为0。

我们简称这三位为

，下面用一个实例来展示一下它的作用

保护位

x = 1.000..00 * 2^1 y = 1.11...11 * 2^0

// without guard bits

x = 1.000...00 * 2^1

- y = 0.111...11 * 2^1

z = 0.000...01 * 2^1

= 1.000...00 * 2^-22

// with guard bits

x = 1.000...00 0000 * 2^1

- y = 0.111...11 1000 * 2^1

z = 1.000...00 0000 * 2^-23

可以很明显的观察到，在多了保护位后，计算结果明显更精确了

但是解决了一个问题后又出现了第二个问题，那就是我们的尾数只有23位，就算我们使用了保护位，也只能保护计算过程中的误差，那么对于结果来说，依然有可能不精确

所以接下来，我们引入舍位

舍位

为了保证运算过程中的最大精确度，我们会把23位+3位附加位的浮点数扩展到64位，相加完后再进行规格化，规格化完后我们要对这三位进行舍入操作。舍入的方法采用就近舍入，中间值4按照奇偶来舍入，确定舍入后，完成尾数的后三位右移操作。

缺陷

我们不妨考虑这个栗子

\\省略对阶

1 010001100 1 010001100 000000

-0 000000111 0 000000111 011110

1 010000101 1 010000100 100010

645.0 644.0

就算有了舍位，在单精度下我们依然很难给出两个浮点数加减的准确答案。

关于浮点数的表示

如果都是规格化的浮点数，那么我们很容易知道我们是没法表示0的，因为最小就是

，所以为了表示0，我们将

和

的和用来表示0，那

该怎么办呢？

为了填上0到

这段上的实数，我们设置了非规格化的浮点数，并且将

到

的数扩展到了整个0到

上，这之间的每个数间隔都是一样的

最大的数为

其中每个数的间隔都是

和最大数之间的距离是

Summary

Floating-point representation

Operations

Addition

Precision Consideration

按：计算机组成课程第四周作业

算法证明

图表 1       浮点数的表示

浮点数的表示如上图所示，我们要做的是将按如上方式存储的两个浮点数相乘，将其结果用如上的方式表示。

符号位只是两者的异或，指数位基本上是两指数的相加，而尾数位就需要在正常的乘法之余考虑到移位(和随之而来的指数为的溢出)和进位的情况。

下面就来讨论一下尾数的运算：

在尾数前补上1，进行无符号数乘法。小数点仅作为算法分析时的记号，实际上不参加运算。用64位的long long类型的数来存储运算结果。

图表 2       尾数相乘(用“画图”程序画的)

如上图所示，后46个数字是结果的小数点后的数据，小数点的的问号处可能只有1，可能是1X(10或11，在算法中没有太大区别)。若问号处只有1，说明已经规格化完成；若是1X，需要将整个数右移一位，质数加一，从而使问号处只有1。

规格化后要进行round操作，如下图所示。如果第23位是0，这一位之前的23位就是所需要的尾数。如果第23位是1，截取这位之前的数，加一(类似于四舍五入)。判断结果是否仍规格化。若没有规格化，从规格化开始继续做，直到找到规格化后的23位尾数。

图表 3       有效数处理(用“画图”程序画的)

程序框图

主要流程如下图所示。在代码注释中的Step1到step5即分别对应流程图中的1至5部分。有所不同的是，因为后面只会增加exponent而不会减少，所以我的程序在step1就判断是否overflow，在step3也判断是否overflow就行了，规格化完成后再判断underflow。

图表 4       浮点乘法程序框图

使用方法

输入两个数，程序用自带的scanf将其作为单精度浮点数读入。通过一系列的无符号数的计算，得到用单精度浮点数表示的结果，将其输出(为了更详细地看到整个数，输出的数小数点后有20位)。如果overflow或underflow，输出的是相应的提示。为了验证计算是否正确，在输出的数据后面有个括号，里面放了是直接让机器进行的单精度浮点数的乘法的结果。

特殊处理

特殊输入数据

有些输入的数据是在太大或太小，以至于数据本身就已经overflow或underflow，从而根本无法存贮到单精度浮点数中，不是合法的单精度浮点数。对于这样的数，我们不应该对其提供任何服务，不然只会闹出笑话。

对于零要特殊处理，因为0在浮点数中的表示并不是按照一般的规范来的。

溢出

这里的溢出指的是乘数和被乘数本身是合法的单精度浮点数，但相乘后造成了溢出。

因为后面只会增加exponent而不会减少，所以与上面的流程图不同，我的程序在step1就判断是否overflow，在step3也判断是否overflow就行了，规格化完成后再判断underflow。

数位扩展

整个程序都是按照单精度的浮点数来写的，数位扩展可能有点麻烦。如果要扩展为双精度的，需要更改代码中的大量参数。鉴于浮点数一般也就只有单精度和双精度两种，我觉得不是太必要提高做成可扩展的。

实例分析

图表 5       大数的乘法

图表 6       overflow

图表 7       underflow(本程序判为underflow，而C语言自己的浮点数乘法计算结果为0.000000)

图表 8       对0的乘法

图表 9       很大的与很小的正数相乘

图表 10     观察精确度

图表 11     对负数乘法的测试

图表 12     对于无法用浮点数表示的输入数据

图表 13     对于过小的数据

结果分析

我们对过于大的数据和过于小的数据，以及精度太大的数据都做了试验。从上述结果可以看到，这个程序的结果与C程序自带的浮点数乘法结果基本吻合。说明这个程序还是比较可靠的。

通过对精度的测试，我们发现浮点数虽然能表示较大较小的数据，但有效数字还是有限的。如-100与0.6相乘的结果，整数部分60是对的，小数点五位后就出现了偏差。一连串8与1相乘，结果的头部只有7个8。可以大致地认为浮点数的表示与计算的十进制有效数位是7位。

可以更进一步的工作

l  程序中对于无符号整数的加法、乘法，对于exponent的有符号整数运算都是直接使用的，其实可以调用之前写的相关程序，做成更加本质的浮点数乘法。

Source Code

#include

union unionfloat{

float f;

unsigned long l;

};

union unionfloat x,y,product;

/*

disp(union unionfloat x)

{

long i;

unsigned long m;

printf("%f:",x.f);

m=0x80000000;

for(i=0; i<32; i++){

if((i&7)==0)printf(" ");

if(m&x.l)printf("1");

else printf("0");

m>>=1;

}

printf("\n");

}*/

int mul(void){

//Step 1

int exponent;

exponent=(x.l&0x7F800000)>>23;

exponent-=0x7F;

exponent+=(y.l&0x7F800000)>>23;

exponent-=0x7F;

if(exponent>=128) return 1;

//Step 2

unsigned long long raw_product;

unsigned long long raw_x=x.l&0x007FFFFF|0x00800000;

unsigned long long raw_y=y.l&0x007FFFFF|0x00800000;

int carry;

raw_product = raw_x * raw_y;

/*int i=64;

unsigned long long p=0x8000000000000000;

while(i--){

if(i%4==3)printf(" ");

if(raw_product&p)printf("1");

else printf("0");

p=p>>1;

}

printf("\n");*/

//相乘结果若小数点前有两位

do{

//Step 3

if(raw_product>>47){

raw_product = raw_product>>1;

exponent++;

if(exponent>=128) return 1;//printf("Overflow");

}

//Step 4

if(raw_product & 0x0000000000400000){

raw_product = ((raw_product>>23)+1)<<23;

}

}while(raw_product>>47);//still normalized?

//Step 5

product.l = ((x.l>>31)^(y.l>>31))<<31;

if(exponent<=-127)return -1;//printf("Underflow");

//-127

product.l = product.l | ((unsigned long)(exponent+127))<<23;

product.l = product.l|((unsigned long)((raw_product&0x00003FFFFF800000)>>23));

//disp(product);

printf("%.20f(%.20f)\n",product.f,x.f*y.f);

return 0;

}

int main(void){

while(~scanf("%f%f",&x.f,&y.f)){

//printf("Input %.10f %.10f\n",x.f,y.f);

//disp(x);disp(y);

if(x.l==0||y.l==0)

printf("%f(%f)\n",0,x.f*y.f);

else if((((x.l>>23)&0x000000FF)==0) || (((x.l>>23)&0x000000FF)==0x000000FF))

printf("The first number is unavailable\n");

else if((((y.l>>23)&0x000000FF)==0) || (((y.l>>23)&0x000000FF)==0x000000FF))

printf("The second number is unavailable\n");

else{

switch(mul()){

case 1:printf("Overflow(%f)\n",x.f*y.f);break;

case -1:printf("Underflow(%f)\n",x.f*y.f);break;

default:break;

}

}

}

return 0;

}