浮点数及其乘法运算(基于IEEE754标准)

浮点数及其乘法运算(基于IEEE754标准)

---

背景：
最近心血来潮想学习下浮点数乘法，基于IEEE754标准完成一个浮点数乘的代码实现；且笔者在查阅相关资料过程中，发现很多关于相关知识点的介绍零零散散，并没有太多较为完整或详细的笔记。于是，决定对相关知识点进行整合，进行一个完整的浮点数乘法实现。

---

说明1：如果文章有误，欢迎大家指出、讨论，笔者也会积极改正，希望大家一起进步！
说明2：为了书写方便，本文中的一些符号相关的写法，可能是基于verilog语言的写法。另有一些流程，为了更直观的的表达，用的是verilog伪代码的方式进行展示。若给大家的阅读带来不便，还请多多包涵。

---

---

前言

浮点数，Floating Point ，相比于日常所接触到的定点数，浮点数的小数点是可以改变的，类似于科学计数法。常见为浮点数类型包括单精度浮点数(Single Floating Point， binary 32 format)、双精度浮点数(Double Floating Point， binary 64 format)，以及其它类型的表示方式。

本文主要围绕32位二进制的单精度浮点数展开研究与实现。

---

一、单精度浮点数(binary 32)

1.binary 32 format

一个binary 32浮点数A[31:0]分为三部分，从左到右依次：符号位(sign，1 bit，对应的是A[31] )、阶码位(exponent，8 bit，对应的是A[30:23] )、尾数位(fraction，23 bit，对应的是A[22:0])。
由于此时exponent的位宽为8 bit ，则 偏移bias = 2^(8-1)-1 =127 。

wire [31:0] A;
assign A[31:0] = {sign,exponent,fraction};

//   sign: 1 bit ，0表示负数、1表示正数;
//	 exponent: 8 bit ，则exponent-bias的取值范围是[-127, 128];
//	 fraction: 23 bit，为浮点数小数部分的二进制表示；默认整数位为1，但不直接写到binary32中;

则，该binary 32浮点数的十进制表示为：

注：
双精度浮点数，binary 64，包括符号位(sign，1 bit)、阶码位(exponent，11 bit)、尾数位(fraction，52 bit)，且偏移bias = 2^(11-1)-1 =1023。
单精度与双精度，主要是位数不一样，其它基本一致，或者说是可以进行类推。

2.特殊情况

由上可知，可以调整sign、exponent、fraction的二进制数值，即可得到对应不同的十进制数值。但是，在exponent、fraction赋值过程中，也会随着赋值的不同，产生一些特殊值。
在进行浮点数计算时，也会产生一些特殊情况，需要进行相应的判断与处理。

2.1.特殊值

2.1.1.特殊值分类

阶码位的取值变化，需要留意。

// 注：Verilog伪代码，不是规范的代码书写，只是为了方便 阅读 与 区分。

input [31:0] A;

if   A[30:23]==8’b00000000   &&   A[22:0]==23'd0   begin
	即exponent、fraction的值都为0;
	此时，浮点数A为0，记作zero；
	
end else if   A[30:23]==8’b00000000   &&   A[22:0]!=23'd0   begin
	即exponent的值为0，fraction的值不为0;
	此时，浮点数A为非规格数，记作denormal；
	引入非规格数，可以获得更小数；
	此时，尾数位fraction是0.开头，而非1.开头。
	（0.开头即非规格化浮点数，1.开头即规格化浮点数，注意区分，详见2.1.2节）。
	
end else if    A[30:23]==8’b11111111   &&   A[22:0]==23'd0   begin
	即exponent的值为127，fraction的值为0;
	此时，浮点数A为无穷大，记作infinity；
	注：无穷大可分为 正无穷大 和 负无穷大 ，根据符号位的值而定。
	
end else if   a[30:23]==8’b11111111   &&    a[22:0]!=23'd0   begin
	即exponent的值为127，fraction不为0;
	此时，浮点数A为非数，记作nan，Not A Number；
	
end

2.1.2.非规格数与规格数的取值范围

最大规约数： (-1)^s * (2^127) * 1.111 1111 1111 1111 1111 1111，其阶码位为8’b11111110。
最小规约数： (-1)^s * 2^(-126) * 1.000 0000 0000 0000 0000 0000，其阶码位为8’b00000001。
最大非规约数：(-1)^s * 2^(-126) * 0.111 1111 1111 1111 1111 1111，其阶码位为8’b00000000。
最小非规约数：(-1)^s * 2^(-126) * 0.000 0000 0000 0000 0000 0000，其阶码位为8’b00000000。

“如果浮点数的指数部分的编码值是0，分数部分非零，那么这个浮点数将被称为非规约形式的浮点数。一般是某个数字相当接近零时才会使用非规约型式来表示。 IEEE 754标准规定：非规约形式的浮点数的指数偏移值比规约形式的浮点数的指数偏移值小1。例如，最小的规约形式的单精度浮点数的指数部分编码值为1，指数的实际值为-126；而非规约的单精度浮点数的指数域编码值为0，对应的指数实际值也是-126而不是-127。实际上非规约形式的浮点数仍然是有效可以使用的，只是它们的绝对值已经小于所有的规约浮点数的绝对值；即所有的非规约浮点数比规约浮点数更接近0。规约浮点数的尾数大于等于1且小于2，而非规约浮点数的尾数小于1且大于0。”(引用来源：https://baike.baidu.com/item/IEEE%20754/3869922?fr=aladdin，转侵删)

关于非规格化数取值范围的详细分析，可参考https://zhuanlan.zhihu.com/p/564458140。

2.1.3.特殊值的出现条件

a. nan(非数)
1 任何数值与nan作运算，其运算结果均为nan
2 对负数开平方之类的运算
3 infinity-infinity、infinity/infinity
4 infinity * 0、infinity/0、0/0
注：nan != nan、nan^0=1

b. infinity(无限大)
1 任意正数N除以零，即N/0。
2 任意正数N乘以无穷大，即infinity * N
3 infinity+infinity、infinity * infinity

2.2.溢出( underflow和overflow )

“当指数和有效数都是0是，就是数学上的0，然而如果指数是0，有效数不是0，则称为underflow，意思非常接近于0，无法用正规的浮点来表示，这种就称为非正规化浮点。”(引用来源：https://zhuanlan.zhihu.com/p/558540055，转侵删)

反之，overflow则意味着数值过大，超出正常的浮点数范围，亦无法用正规的浮点数表示，这是被称为无穷大。

3.IEEE754舍入方式

为了对浮点数数进行舍入计算，且尽可能地减小计算误差，IEEE754标准给出了以下4 种舍入方式。

同时，假设需要对浮点数计算结果C的后n位（即，C[n-1:0]）进行舍入。

3.1.就近舍入

就近舍入，也叫 向偶数舍入，Round to nearest (even)，对应于十进制下的“四舍五入”。
赛灵思在对浮点数进行乘法运算后，所采用的便是这种舍入方式，具体流程如下。

1   如果C[n-1:0] > {1'b1,{(n-1)'b0}}
		即，C[n-1:0] 大于 2^n，则舍去低n位时，产生进位，需要C[n]+1'b1；
2   如果C[n-1:0] < {1'b1,{(n-1)'b0}}
		即，C[n-1:0] 小于 2^n，则舍去低n位时，不产生进位，C[n]不变；
3   如果C[n-1:0] = {1'b1,{(n-1)'b0}}
		即，C[n-1:0] 等于 2^n，此时：
			如果C[n]=1’b1，则产生进位，需要C[n]+1'b1；
			如果C[n]=1’b0，则不产生进位，C[n]不变；

3.2.朝0舍入

朝0舍入，Round toward zero，即直接截尾，将C[n-1:0]舍去。
此时，无论C为正数还是负数，舍入后的c值相比于舍入前，在数值上更靠近0点了些。

3.3.朝正无穷舍入

朝正无穷舍入，Round up toward +∞ 。

1   如果C是负数
		即，C小于0，则直接舍去c的低n位C[n-1:0]，不需要C[n]+1'b1；
2   如果C是正数
		即，C大于0，则
			当C的低n位不全为零（即，C[n-1:0]>0）时，产生进位，C[n]+1；否则，直接舍去C[n-1:0]。

3.4.朝负无穷舍入

朝负无穷舍入，Round down toward −∞ 。

1   如果C是负数
		即，C大于0，则
		当C的低n位不全为零（即，C[n-1:0]>0）时，产生进位，C[n]+1；否则，直接舍去C[n-1:0]。
2   如果C是正数
		即，C小于0，则直接舍去c的低n位C[n-1:0]，不需要C[n]+1'b1；

二、浮点数乘法

1.思路

假设需要对两个浮点乘数A[31:0]和B[31:0]进行乘法运算，得到C，则具体流程如下：
C = A * B
= { A_sign, A_exponent, A_fraction} * { B_sign, B_exponent, B_fraction};
其中：

1. C_sign = A_sign ^ B_sign;
2. C_exponent = A_exponent + B_exponent -127;
3. 1.C_fraction = 1.A_fraction * 1.B_fraction;

2.乘法算法

// 注：Verilog伪代码，不是规范的代码书写，只是为了方便 阅读 与 区分。

C_sign = A[31] ^ B[31]; // 这里 ^ 是verilog里面的运算符，表示异或运算，即A[31]与B[31]

if( A或B为 nan ) begin
	则返回nan;
end else if( A或B为 infinity ) begin
	if( 另一个数位 zero或denormal ) begin
		则返回nan;
	end else if( 另一个数为 非zero或非denormal )
		则返回infinity ;
	end
end else if( A或B为 zero或denormal ) begin
	则返回zero;
end else begin
	C_exponent = A[30:23] + B[30:23] - 127;
	1.C_fraction = 1.A[22: 0] * 1.B[22:  0];
	a 规格化处理
	b 使用Round to nearest (even)进行舍入计算
	c 溢出判断：
		若underflow，则返回zero;
		若overflow，则返回nan;
		若无溢出，则返回正常数;
end

3.代码示例

	fp_data  fp_data_inst(	    //   数据读取模块，对存在本地的fp乘法数据进行读取
		.Clk(Clk),
		.rst_N(rst_N),
		.A_S(A_sign),
		.A_E(A_exp[7:0]),   
		.A_M(A_man[23:0]),    //   A_sign、A_exp、A_man，是乘数A的数据，通过matlab生成
		.B_S(B_sign),
		.B_E(B_exp[ 7:0]),
		.B_M(B_man[23:0]),   //   B_sign、B_exp、B_man，是乘数B的数据，通过matlab生成
		.A_B(a_b)                  //   事先通过matlab算出来的A*B结果
	);
	fp_multi fp_multi_inst(	   //   乘数A和乘数B的数据，传递给浮点乘法器
		.A_SIGN(A_sign),
		.A_EXP(A_exp[7:0]),
		.A_MAN({1'b1,A_man[22:0]}), 
		.B_SIGN(B_sign),
		.B_EXP(B_exp[7:0]),
		.B_MAN({1'b1,B_man[22:0]}), 
		.FPM_OUT(fpm_out[31:0])   //   通过浮点乘法器算出来的A*B结果
	);
	
assign PassStatus = (a_b == fpm_out)? 1 : 0 ; //   将matlab算出的与浮点乘法算出的结果进行对比

4. RTL仿真图(vivado)

局部波形图：

整体波形图：

上述图中，
clk、rst_N，分别是时钟信号、复位信号；
addr，是数据的序号，用于表示第x组A与B的数据下标；
a_b，是通过matlab算出来的A*B结果，用于与浮点数乘法计算出来的结果进行比对；
fpm_out，是浮点数乘法计算出来的结果；
PassStatus，是浮点数乘法是否通过计算。如果a_b与fpm_out相等，则PassStatus=1；否则，为0。

从RTL仿真结果来看，PassStatus的一直为1，这意味浮点数乘法器的功能应该是通过了的。

---

最后

在网上看到的浮点数计算的在线工具：
http://weitz.de/ieee/
https://www.toolhelper.cn/Digit/FractionConvert

后续会把本文的verilog代码下载链接挂出来，欢迎大家来下载。同时，本文若有不足之处，欢迎大家在评论区进行交流、指正。