3

关于浮点数的运算,彻底弄懂为何 101.4 - 80.0 != 21.4

 3 years ago
source link: http://yangxikun.com/2020/08/18/float.html
Go to the source link to view the article. You can view the picture content, updated content and better typesetting reading experience. If the link is broken, please click the button below to view the snapshot at that time.
neoserver,ios ssh client

关于浮点数的运算,彻底弄懂为何 101.4 - 80.0 != 21.4

18 August 2020

关于浮点数运算精度问题网上已有很多资料解释,但笔者还是想再造下“轮子”,加深理解。

先看一段简单的代码:

package main

func main() {
	var a float32 = 101.4
	var b float32 = 80.0
	var c float32 = 21.4
	var d float32 = 21.4
	if a-b == c {
		println("a-b = c")
	}
	if c == d {
		println("c = d")
	}
}

上述代码的运行结果为:

c = d

为何 a-b != c ? 如果我们直接在纸上用笔算,或者直接脑算都能得出 101.4 - 80.0 = 21.4。那么为何计算机算错了呢?

一句话说明原因:计算机在存储浮点数时存在精度误差。

之所以我们能算出 101.4 - 80.0 = 21.4,是因为在纸上或者脑子里可以精确的表示 101.4 和 80.0。那么为何计算机表示浮点数时会存在精度误差,并且这个精度误差是如何影响到计算结果的?

十进制整数可以直观地转为二进制数,在计算机中采用补码的形式存储,方便做加减乘除。而浮点数因为有小数点,计算机只认识0、1,可不认识小数点,一开始各家计算机公司的各个型号的计算机,有着千差万别的浮点数表示。为了解决数据交换、计算机协同工作问题,IEEE 设计了一套业界通用的标准,即 IEEE 二进制浮点数算术标准(IEEE 754)。

IEEE 754 规定二进制浮点数 V 可以表示为如下等式:

V = (-1)^S * M * 2^E

  • (-1)^S 为符号位,当 S = 0,V 为正数,当 S = 1,V 为负数;
  • M 表示有效数的小数部分
  • 2^E 为指数部分

IEEE 754 规定的浮点数在计算机中存储时的内存布局如下:

exponent = E + 127,即实际存储时,需要把 E 加上 127 当表示 < 1 的十进制小数时,E < 0,E + 127 后的取值范围为 [0, 255],之所以要这么转换,应该是为了方便比较浮点数大小

对于 32bit 浮点数,IEEE 754 规定:sign(符号位)占 1 位,exponent(指数位)占 8 位,fraction(有效数的小数部分)占 23 位。

一个十进制小数是如何转换为 IEEE 754 规定的浮点数表示存储在计算机中的?举个例子,十进制小数 100.25 的转换过程如下:

首先,对于整数部分,即 100,采用“除2取余,逆序排列”法可得:1100100

100 / 2 = 50 余 0
50 / 2 = 25 余 0
25 / 2 = 12 余 1
12 / 2 = 6 余 0
6 / 2 = 3 余 0
3 / 2 = 1 余 1
1 / 2 = 0 余 1 // 商为 0 停止

对于小数部分,即 0.25,采用“乘2取整,顺序排列”法可得:01

0.25 * 2 = 0.5 取出整数部分 0
0.5 * 2 = 1 取出整数部分 1 // 无小数部分停止

所以十进制小数 100.25 转换为二进制小数 1100100.01 = 1.10010001 * 2^6。那么 S=0,M=10010001,E=6,存储为二进制格式 01000010110010001000000000000000。

可以看到 M 只取了小数部分,整数 1 被省略掉了。这是因为 IEEE 754 规定,exponent 在 [1, 254] 范围的时候,那么这个浮点数将被称为规约形式的浮点数,在这种情况下,整数部分始终为 1,因此在实际存储的时候可以省略。 当 exponent = 0 或 255 的时候被用来表示一些特殊的数,详情可查看wikipedia IEEE 754

回到文章开头的代码,如果将 101.4 转换为二进制小数,会发现小数部分 0.4 的二进制小数是无限循环的:0.0110,0110循环。

0.4 * 2 = 0.8
0.8 * 2 = 1.6
0.6 * 2 = 1.2
0.2 * 2 = 0.4
......

但 IEEE 754 中定义的 M 只有 23 bit,所以在存储的时候会发生截断。测试了下 Golang 中的截断规则:当截断后第一位为 1 时,进位,否则直接舍弃。

func main() {
	// 逗号表示之后的 bit 会被截断掉
	var f float32 = 16777218.0 // 二进制表示为:100000000000000000000001,0
	printFloat32Binary(f)      // IEEE 754 格式存储的二进制内容:01001011100000000000000000000001 = 16777218.0
	f = 16777219.0             // 100000000000000000000001,1
	printFloat32Binary(f)      // 01001011100000000000000000000010 = 16777220.0
	f = 33554437.0             // 100000000000000000000001,01
	printFloat32Binary(f)      // 01001100000000000000000000000001 = 33554436.0
	f = 67108875.0             // 100000000000000000000001,011
	printFloat32Binary(f)      // 01001100100000000000000000000001 = 67108872.0
	f = 67108879.0             // 100000000000000000000001,111
	printFloat32Binary(f)      // 01001100100000000000000000000010 = 67108880.0
	f = 67108871.0             // 100000000000000000000000,111
	printFloat32Binary(f)      // 01001100100000000000000000000001 = 67108872.0
}

func printFloat32Binary(f float32) {
	fmt.Printf("%f\n", f)
	byteSliceRev := *(*[4]byte)(unsafe.Pointer(&f))
	for i := 0; i < 4; i++ {
		b := byteSliceRev[3-i]
		for j := 0; j < 8; j++ {
			if b&0b10000000 == 0 {
				print(0)
			} else {
				print(1)
			}
			b = b << 1
		}
	}
	print("\n")
}

那么来看下 101.4,80.0,21.4 的二进制存储(逗号分隔了 sign,exponent,fraction):

101.4: 0,10000101,10010101100110011001101
 80.0: 0,10000101,01000000000000000000000
 21.4: 0,10000011,01010110011001100110011

现在模拟 101.4 - 80.0 的过程,因为它们的阶码是相同的,直接把 M 相减可得:0.01010101100110011001101,为了确保最高有效位为 1,需要把阶码减 2,得到最终结果:0,10000011,01010110011001100110100,显然 fraction 与 21.4 不同。

上述计算过程是笔者直观地算出来的,但跟实际结果是一致的。

package main

import (
	"fmt"
	"unsafe"
)

func main() {
	var a float32 = 101.4
	var b float32 = 80.0
	var c float32 = 21.4
	fmt.Printf("%.32f\n", a-b) // 21.40000152587890625000000000000000
	printFloat32Binary(a - b)  // 01000001101010110011001100110100
	fmt.Printf("%.32f\n", c)   // 21.39999961853027343750000000000000
	printFloat32Binary(c)      // 01000001101010110011001100110011
}

func printFloat32Binary(f float32) {
	byteSliceRev := *(*[4]byte)(unsafe.Pointer(&f))
	for i := 0; i < 4; i++ {
		b := byteSliceRev[3-i]
		for j := 0; j < 8; j++ {
			if b&0b10000000 == 0 {
				print(0)
			} else {
				print(1)
			}
			b = b << 1
		}
	}
	print("\n")
}

所以 101.4 - 80.0 != 21.4 的原因就在与 101.4 和 21.4 在计算机浮点数中无法被精确表示。那么对于能够精确表示的浮点数,实际上就不会存在问题:

func main() {
	var a float32 = 1.5
	var b float32 = 0.5
	var c float32 = 1.0
	if a-b == c {
		println("a - b = c")
	}
	fmt.Printf("%.32f\n", a-b) // 1.00000000000000000000000000000000
	printFloat32Binary(a - b)  // 00111111100000000000000000000000
	fmt.Printf("%.32f\n", c)   // 1.00000000000000000000000000000000
	printFloat32Binary(c)      // 00111111100000000000000000000000
}
a - b = c
1.00000000000000000000000000000000
00111111100000000000000000000000
1.00000000000000000000000000000000
00111111100000000000000000000000

About Joyk


Aggregate valuable and interesting links.
Joyk means Joy of geeK