现代电子计算机通常使用浮点数表示小数,最常用的数据类型是 double。 计算机中浮点数表示形式如下(^ 表示乘方):

m * 2^e

浮点数类型使用注意

早几年的时候,如果我们在路边摊买一本编程的书,书上会温馨提醒我们不要直接使用 == 比较两个浮点数(如 double)是否相等。

  • 同一个数值可能存在多种表示形式,而底层语言可能是按 bit 进行 == 相等比较,同一数值的多种表示形式可能被误认为不相等(?)。

  • 实数有无穷多个,而计算机浮点数只能使用有限个 bit 表示(如 double 通常使用 64 bit),只能表示有限个数值, 因此必然有一些数值是无法使用计算机浮点数表示的。

    上式中 2^e 影响浮点数的精度,最小精度之间的数值均不能表示,只能保存为计算机能表示的与之最接近的数值。 多次计算时,每次计算结果都可能丢失一定精度,因此最终结果数值可能与理论值存在一定误差。

因此,比较两个浮点数是否相等的正确姿势是计算他们的差值,如果差值小于某个可以容忍的误差范围(如 0.00001),则认为其相等。

C/C++ 代码示例如下:

double d1;
double d2;
double precision = 0.00001;
// d1 == d2 ?
if (fabs(d1 - d2) <= precision) {
}

误差精度范围如何确定呢?上例中的 0.00001 显然给的很随意。 遗憾的是精度范围跟具体计算数值范围有关,甚至可能跟具体应用场景有关,没有一个普适的值。

epsilon 值

  • C 语言中给出了一个参考精度值,用宏 DBL_EPSILON 表示(参考 float.h)。
  • C++ 进行了类型安全的封装,用 std::numeric_limits<double>::epsilon() 表示(参考 numeric_limits)。

尝试打印 epsilon 值,测试代码如下:

#include <stdlib.h>
#include <limits.h>
#include <float.h>
#include <math.h>
#include <stdio.h>

#include <limits>
#include <string>
#include <sstream>

int main() {
	double d1 = DBL_EPSILON;
	double d2 = std::numeric_limits<double>::epsilon();
	double precision = pow(2, -100);
	bool equal = false;
	// d1 == d2 ?
	if (fabs(d1 - d2) <= precision) {
		equal = true;
	}
	printf("d1=%lf d2=%lf precision=%lf equal=%d\n", d1, d2, precision, (int) equal);

	std::ostringstream ss;
	ss.setf(std::ios_base::boolalpha);
	ss << "d1=" << d1 << " d2=" << d2 << " precision=" << precision << " equal=" << equal;
	std::string s = ss.str();
	printf("%s\n", s.c_str());

	return 0;
}

输出结果如下:

d1=0.000000 d2=0.000000 precision=0.000000 equal=1
d1=2.22045e-16 d2=2.22045e-16 precision=7.88861e-31 equal=true
  • C 语言 printf "%lf" 默认限制了输出精度(?),未能有效输出 epsilon 数值。
  • C++ 默认自动使用科学计数法输出,设置 boolalpha 还可以用文本方式输出 bool 值。

稍微看一下文档就会知道,C/C++ 中 epsilon 值是浮点数表示数值 1 时,所能表示的最小精度值。

计算机浮点数表示形式

wiki 可以知道, 现代计算机的 64 bit double 数据类型通常表示如下(IEEE 754):

IEEE_754_Double_Floating_Point_Format.svg

  • 最高位 1 bit 符号位。符号位单独表示,数值 0 存在两种浮点数表示形式:+0-0
  • 高位 11 bit 指数位。数值 e0 范围 [0, 2047],实际指数 e 使用 e0 - 1023(范围 [-1023, 1024]),最小值和最大值留作特殊用途,有效指数 e 范围 [-1022, 1023]
  • 低位 52 bit 小数位。表示 m 的二进制小数部分,m 取值为 1.x,整数部分约定固定为 1。

前面提到的浮点数表示可进一步细化为:

m * 2^e

这里有个很巧妙的设计:

  • m 取值范围为 [1, 2)(最大值 2 - e-52`),给定指数 e,可表示浮点数范围为 [2e, 2e+1)
  • 数值继续增大时,指数 e 加 1,m 回到最小值 1。如此周而复始。

如此可知,给定数值 d,即可确定指数 e2e <= d < 2e+1),进一步确定 m 的值,从而得出其浮点数表示。 如数值 1 浮点数(double)表示如下:

0 01111111111 0000000000000000000000000000000000000000000000000000 ≙ +1 * 21023 - 1023 = +1 * 20 = 1

此时 2e == 20 == 1m 最低小数位为 2-52, 这个值即是前面提到的 epsilon 值(即浮点数表示数值 1 时的最小精度),约 2 * 10-16

  • C 语言中同样定义了计算机浮点数使用的进制 FLT_RADIX(通常为常量 2),有效指数范围 DBL_MIN_EXPDBL_MAX_EXP,常规数值范围 DBL_MINDBL_MAX
  • C++ 进行了类型安全的封装了,int 常量 std::numeric_limits<double>::radix, min_exponent, max_exponent,double 函数 min()max()
  • Java 中定义了相关常量 Double.MIN_EXPONENT, MAX_EXPONENT, MIN_NORMAL, MAX_VALUE, MIN_VALUE, POSITIVE_INFINITY, NEGATIVE_INFINITY, NaN 等。

我们知道 e 最小值为 -1022 (MIN_EXPONENT),m 最小值为 1, 因此常规最小值 MIN_NORMAL 为:

0 00000000001 0000000000000000000000000000000000000000000000000000 ≙ +1 * 2−1022 ≈ 2.2250738585072014 * 10−308

此时 m 最低位为 2-52,此时浮点数的最小精度为 2-1022-52 = 2-1074。 Java 中定义了另一个最小值常量 MIN_VALUE 表示这个值。 计算机浮点数如何表示这个值呢? 前面提到原始指数最小值(指数 bit 全 0)留作特殊用途,即在此使用,此时有效指数 e 仍为 -1022,但 m 整数部分变成 0, 这样计算机还能表示 2^52 个小于 MIN_NORMAL 的浮点数,最小的一个正数即为 MIN_VALUE。 这也是计算机浮点数(double)所能表示的最小精度。

0 00000000000 0000000000000000000000000000000000000000000000000001 ≙ +2-52 * 2−1022 = 2−1074 ≈ 4.9 * 10−324

C++ 打印相关常量:

void printDoubleConst() {
	int const radix = std::numeric_limits<double>::radix;
	int const minExp = std::numeric_limits<double>::min_exponent;
	int const maxExp = std::numeric_limits<double>::max_exponent;
	double const min = std::numeric_limits<double>::min();
	double const max = std::numeric_limits<double>::max();

	std::ostringstream ss;
	ss.setf(std::ios_base::boolalpha);
	ss << "radix=" << radix << " minExp=" << minExp << " maxExp=" << maxExp
			<< " min=" << min << " max=" << max;
	std::string s = ss.str();
	printf("%s\n", s.c_str());
}

结果如下:

radix=2 minExp=-1021 maxExp=1024 min=2.22507e-308 max=1.79769e+308

Java 打印相关常量:

public static void main(String[] args) throws Exception {
	String s = String.format("minExp=%d maxExp=%d minNormal=%s max=%s min=%s", 
			Double.MIN_EXPONENT, Double.MAX_EXPONENT, Double.MIN_NORMAL, Double.MAX_VALUE, Double.MIN_VALUE);
	System.out.println(s);
}
  • Java 中使用 "%f" 与 C 中 "%lf" 类似,默认会限制精度。但使用 "%s" 转字符串会自动使用科学计数法,与 C++ 类似。

结果如下:

minExp=-1022 maxExp=1023 minNormal=2.2250738585072014E-308 max=1.7976931348623157E308 min=4.9E-324
  • 可看到 C/C++ 和 Java 对最小指数的常量定义存在轻微差异。

给定数值 d,计算浮点数指数 e

  • C 语言中可使用 double logb(double x),返回值为取整后的 double,这个函数相当于 floor(log2(x)),但速度更快。
  • Java 中可使用 int Math.getExponent(double d),返回值为整数。

前面提到:表示不同范围的数值时,浮点数所能表示的最小精度值是不同的。 很显然最小精度即 2e-52。 当最小精度为 1 时,e52,由此可知使用 double 表示小于 2^53 的整数时都不会丢失精度(精度为 1 时可表示 2^52 个整数,精度更小时还可表示另外 2^52 个整数), 但表示更大的整数时就会丢失精度,且精度为 2 的整数次方。 因此 C/C++/Java 中将 int 转为 double 是安全的,将 long 转为 double 则是不安全的,可能丢失精度,需要特别注意。

控制和访问浮点数 bit 位

C/C++ 可巧妙利用 union 类型呈现 double 类型的内部表示。 代码示例如下:

union Double {
	double val;
	uint64_t longBit;
	struct Bit {
		static_assert(sizeof(double) == 8, "double must be 64bit (8 byte)");
		static int constexpr Sign = 1;
		static int constexpr Exponent = 11;
		static int constexpr Fraction = 52;
	};
	struct {
		uint64_t fraction:Bit::Fraction;
		uint16_t exponent:Bit::Exponent;
		uint8_t sign:Bit::Sign;
	};
	// ... ...
}
static_assert(sizeof(Double) == sizeof(double), "size of Double must same as double");
  • static_assert 是 C++11 引入的编译时断言,检查确保 Double 是对 double 的单纯包装,没有增加任何额外的 bit 位。
  • 经典 C++ (C++98) 自动把某些 const 优化(可选?)为编译时常量,C++11 对 const 进一步细化,引入 constexpr 明确表示编译时常量(必须)。
  • union 类型中不允许定义 static 静态字段,在嵌套 struct Bit 中定义。
  • union 类型中的匿名 struct 自动成为 union 成员,其字段暴露到 union 中直接访问。
  • 使用位域将 double 内部表示拆分为 sign, exponent, fraction 3 段,全部使用无符号整数表示。 因为内存地址是从低位到高位,所以这里 struct 定义是从 fraction 到 sign。

添加构造函数,支持自定义构造 double 值:

explicit Double(double const& val) {
	this->val = val;
}

explicit Double(uint64_t longBit) {
	this->longBit = longBit;
}

Double(uint8_t sign, uint16_t exponent, uint64_t fraction) {
	this->sign = sign;
	this->exponent = exponent;
	this->fraction = fraction;
}
  • 构造函数添加 explicit 关键字,避免参数类型自动转换(隐式)转换为 Double 类型。

支持查看和打印 double 内部结构:

int exp() {
	return exponent ? ((int) exponent) - 1023 : -1022;
}

std::ostream& printBit(std::ostream& out) const {
	out << std::bitset<Bit::Sign>(sign) << " "
			<< std::bitset<Bit::Exponent>(exponent) << " "
			<< std::bitset<Bit::Fraction>(fraction) << " "
			;
	return out;
}
  • 如前文所说,有效指数 e 为 原始 exponent - 1023,但原始 exponent 为 0 时特殊表示最小指数 -1022
  • 使用 std::bitset<N> 打印 bit 位二进制表示。其构造函数入参 bit 位多于 N 时,使用最低 N 位。

尝试打印一组数据:

void printDoubleBit() {
	std::ostringstream ss;
	double dList[] = {
			+0.0, -0.0, Double((uint64_t) 1L).val, Double(0, 0, ((uint64_t) 1L) << 51).val, std::numeric_limits<double>::min(),
	};
	for (double d : dList) {
		Double& wrap = Double::wrap(d);
		ss << std::setw(12) << d << " -> ";
		wrap.printBit(ss)
				<< " e=" << std::setw(5) << wrap.exp()
				<< std::endl;
	}
	std::string s = ss.str();
	printf("%s", s.c_str());
}

结果如下:

           0 -> 0 00000000000 0000000000000000000000000000000000000000000000000000  e=-1022
          -0 -> 1 00000000000 0000000000000000000000000000000000000000000000000000  e=-1022
4.94066e-324 -> 0 00000000000 0000000000000000000000000000000000000000000000000001  e=-1022
1.11254e-308 -> 0 00000000000 1000000000000000000000000000000000000000000000000000  e=-1022
2.22507e-308 -> 0 00000000001 0000000000000000000000000000000000000000000000000000  e=-1022
  • C/C++ 字面值区分 +0.0-0.0,符号为不一样,其他位全为 0。
  • fraction 取最小正整数 1,其他位全 0,得到最小正浮点数。这个值即为 Java 定义的常量 Double.MIN_VALUE
  • fraction 最高位为 1,其他位全 0,这个值为常规最小值的一半(此时特殊表示有效指数 e 为最小值 -1022,m 值为二进制 0.1,十进制 0.5)。
  • 原始 exponent 为常规最小值 1,其他位全 0,即常规最小值,C++ 中为 min(),Java 中为 Double.MIN_NORMAL

再测试打印一组值:

double constexpr epsilon = std::numeric_limits<double>::epsilon();
double constexpr max = std::numeric_limits<double>::max();
Double maxLess = Double(max);
--maxLess.fraction;
double dList[] = {
		epsilon, 1.0, 1.0 + epsilon, max, maxLess.val, max - maxLess.val, Double(0, 52 + 1023, ~((uint64_t)0L)).val,
};

结果如下:

 2.22045e-16 -> 0 01111001011 0000000000000000000000000000000000000000000000000000  e=  -52
           1 -> 0 01111111111 0000000000000000000000000000000000000000000000000000  e=    0
           1 -> 0 01111111111 0000000000000000000000000000000000000000000000000001  e=    0
1.79769e+308 -> 0 11111111110 1111111111111111111111111111111111111111111111111111  e= 1023
1.79769e+308 -> 0 11111111110 1111111111111111111111111111111111111111111111111110  e= 1023
1.99584e+292 -> 0 11111001010 0000000000000000000000000000000000000000000000000000  e=  971
  9.0072e+15 -> 0 10000110011 1111111111111111111111111111111111111111111111111111  e=   52
  • epsilon 值为 2^-52
  • 1 + epsilon,fraction 最低位加 1(最小精度),其他 bit 位不变。此时打印浮点数时,最低有效位被忽略,打印为 "1"。
  • 原始 exponent 为最大常规值 2046(有效指数 e 为 1023),fraction 为最大值时,为最大浮点数。
  • 在最大值范围区间,最低精度为 2^(1023-52) = 2^971,fraction 最低位减少 1,表示数值减少 2^971
  • 有效指数 e 为 52,fraction 为最大值时,为 double 能表示的不丢失精度的(最小精度 <= 1)最大整数,值为 2^53 - 1
  • 无论表示多大/多小的数,m 有效小数位始终为二进制 52 位 fraction,十进制约 16 位。