Cracking C++(6): 准确打印浮点数

1. 目的

给 float 或 double 类型的变量赋值后， 打印出来的值和赋值时传入的值可能不一样， 也就是说有有精度损失。常见的几个疑问是：

• 为啥有精度损失？
• 为啥 float 类型精度损失这么大， 我们老师以前说是1e-6的？
• 为啥明明有好几位小数， printf 和 cout 只打印出6位？

答案：

• float 和 double 类型是 IEEE 754 标准规定的， 首先要转为二进制表示， 再按格式算出能表示的值， 而转为二进制表示的过程中受限为位数， 存在精度损失的问题；
• 在得到有精度损失的数值后， printf 和 cout 默认打印的精度不是很友好， 可以用 fmt 库执行打印

本文不涉及浮点数转二进制时的 rounding 细节， 考虑的是得到 rounding 后的二进制后， 逐步算出被表示的浮点数数值的过程， 以及各个部分的二进制表示。支持 float, double, fp16 三种类型.

2. 准确打印浮点数： 使用 fmt 库

CMakeLists.txt

cmake_minimum_required(VERSION 3.25)
project(x)
set(CMAKE_CXX_STANDARD 20)
add_executable(testbedcoutRealNumber.cpp
)
add_subdirectory("/Users/chris/work/github/fmt" fmt.out)
target_link_libraries(testbed PUBLIC fmt::fmt)

C++ 关键代码：

    float pi_f32 = 3.1415926;std::cout << fmt::format("{}", pi_f32);

运行结果

printf pi=3.141593
std::cout pi=3.14159
fmt::format pi=3.1415925

3. 准确算出被表示的值

3.1 直观感受IEEE-754: float-toy

如果你关心”为什么赋值和打印结果不一样“， 那就需要按 IEEE 754 标准， 按步骤算出取值。

不妨先直观感受下 π 的 fp16, float32, float64 类型的二进制表示， 以及计算出的结果, 用到的在线工具是 float-toy:

3.2 获取浮点数二进制表示

使用 std::bitset<N> 来表示浮点数的二进制表示。其实整数也可以用它来获取二进制表示。

以 float32 类型为例

class Float32
{
public:explicit Float32(float _f) :f(_f){memcpy(&u, &f, sizeof(float));b = std::bitset<32>(u);}int getSignBit() const{return b[31];}std::bitset<8> getExponent() const{std::bitset<8> exponent;for (int i = 31, j = 7; i >= 24; i--, j--){exponent[j] = b[i-1];}return exponent;}std::bitset<23> getSignificand() const{std::bitset<23> significand;for (int i = 23, j = 22; i >= 1; i--, j--){significand[j] = b[i-1];}return significand;}std::bitset<32> getBinary() const{return b;}private:std::bitset<32> b;unsigned int u;float f;
};int main()
{float pi_f32 = 3.141592653589793;{std::cout << "IEEE 754 single precision example" << std::endl;Float32 r(pi_f32);std::cout << "sign: " << r.getSignBit() << "\n";std::cout << "exponent: " << r.getExponent().to_string() << "\n";std::cout << "significand: " << r.getSignificand().to_string() << "\n";}
}

运行结果如下：

IEEE 754 single precision example
sign: 0
exponent: 10000000
significand: 10010010000111111011011

3.3 float 类型

核心公式是:

V = SP * FP * EP= (-1)^s * M * 2^E

其中 SP 意思是 sign part, 符号部分；
FP 意思是 fraction part, 小数部分；
EP 意思是 exponent part, 指数部分。

M， E 的具体计算可以翻《CSAPP》这本书。这里只考虑常规的浮点数， 也就是说像 NAN， INF 这样的这里没考虑。

对应的代码实现，在 Float32 类类型中增加成员函数

float getValue() const{//return value;//return f;// V = SP * FP * EP//   = (-1)^s * M * 2^E// SP: OKint s = getSignBit();int SP = (s == 0) ? 1 : -1;// FP: OKunsigned int significand = getSignificand().to_ulong();float f = significand * (1.0 / (1 << 10));float FP = 1.0f + f;printf("significand: %d\n", significand);// EP: OKunsigned int e = getExponent().to_ulong();unsigned int Bias = 15; // 2^(k-1) - 1, k = 5unsigned E = e - Bias;float EP = (1 << E);printf("SP: %d\n", SP);printf("FP: %lf\n", FP);printf("EP: %f\n", EP);// TODO: 这里打印的结果， 和 float-toy 对不上// 考虑使用 https://github.com/Maratyszcza/FP16/blob/master/include/fp16/fp16.h 作为验证float V = SP * FP * EP;return V;}

3.4 double 类型

和 float 类型的 getValue() 函数大同小异。

这里的插曲是，原版 float-toy 有 bug，至少对于页面默认显示的 π 的 fp16 类型来说， 结果是错的。具体讨论见 https://github.com/evanw/float-toy/issues/9。

class Float64
{
public:explicit Float64(double _lf) :lf(_lf){memcpy(&u, &lf, sizeof(double));b = std::bitset<64>(u);}int getSignBit() const{return b[63];}std::bitset<11> getExponent() const{std::bitset<11> exponent;for (int i = 63, j = 10; i >= 53; i--, j--){exponent[j] = b[i-1];}return exponent;}std::bitset<52> getSignificand() const{std::bitset<52> significand;for (int i = 52, j = 51; i >= 1; i--, j--){significand[j] = b[i-1];}return significand;}std::bitset<64> getBinary() const{return b;}double getValue() const{// return lf;//return f;// V = SP * FP * EP//   = (-1)^s * M * 2^E// SP: OKint s = getSignBit();int SP = (s == 0) ? 1 : -1;// FP: OKunsigned long long significand = getSignificand().to_ullong();double f = significand * (1.0 / (1ULL << 52));double FP = 1.0 + f;// printf("significand: %d\n", significand);// EP: OKunsigned long long int e = getExponent().to_ullong();unsigned long long int Bias = 1023;unsigned long long E = e - Bias;double EP = (1ULL << E);// printf("SP: %d\n", SP);// printf("FP: %lf\n", FP);// printf("EP: %f\n", EP);double V = SP * FP * EP;return V;}private:std::bitset<64> b;uint64_t u;double lf;
};

3.5 fp16 类型

class Float16
{
public:explicit Float16(float f){memcpy(&u, &f, sizeof(float));std::bitset<32> b32(u);b[15] = b32[31];for (int i = 0; i < 5; i++){b[15 - 1 - i] = b32[31 - 1 - i];}for (int i = 0; i < 10; i++){b[10 - 1 - i] = b32[23 - 1 - i];}}int getSignBit() const{return b[15];}std::bitset<5> getExponent() const{std::bitset<5> exponent;for (int i = 15, j = 4; i >= 11; i--, j--){exponent[j] = b[i-1];}return exponent;}std::bitset<10> getSignificand() const{std::bitset<10> significand;for (int i = 10, j = 9; i >= 1; i--, j--){significand[j] = b[i-1];}return significand;}std::bitset<16> getBinary() const{return b;}float getValue() const{//return value;//return f;// V = SP * FP * EP//   = (-1)^s * M * 2^E// SP: OKint s = getSignBit();int SP = (s == 0) ? 1 : -1;// FP: OKunsigned int significand = getSignificand().to_ulong();float f = significand * (1.0 / (1 << 10));float FP = 1.0f + f;printf("significand: %d\n", significand);// EP: OKunsigned int e = getExponent().to_ulong();unsigned int Bias = 15; // 2^(k-1) - 1, k = 5unsigned E = e - Bias;float EP = (1 << E);printf("SP: %d\n", SP);printf("FP: %lf\n", FP);printf("EP: %f\n", EP);float V = SP * FP * EP;return V;}private:std::bitset<16> b;unsigned int u;float value;
};

3.6 验证

int main()
{float pi_f32 = 3.141592653589793;double pi_f64 = 3.141592653589793;{std::cout << "IEEE 754 single precision example" << std::endl;Float32 r(pi_f32);std::cout << fmt::format("{}", r.getValue()) << "(0x" << std::hex << r.getBinary().to_ulong() << ")\n";std::cout << "sign: " << r.getSignBit() << "\n";std::cout << "exponent: " << r.getExponent().to_string() << "\n";std::cout << "significand: " << r.getSignificand().to_string() << "\n";}if (1){std::cout << "\n";std::cout << "IEEE 754 double precision example" << std::endl;Float64 r(pi_f64);std::cout << fmt::format("{}", r.getValue()) << "(0x" << std::hex << r.getBinary().to_ulong() << ")\n";std::cout << "sign: " << r.getSignBit() << "\n";std::cout << "exponent: " << r.getExponent().to_string() << "\n";std::cout << "significand: " << r.getSignificand().to_string() << "\n";}if (1){std::cout << "\n";std::cout << "IEEE 754 half precision example" << std::endl;Float16 r(pi_f32);std::cout << fmt::format("{}", r.getValue()) << "(0x" << std::hex << r.getBinary().to_ulong() << ")\n";std::cout << "sign: " << r.getSignBit() << "\n";std::cout << "exponent: " << r.getExponent().to_string() << "\n";std::cout << "significand: " << r.getSignificand().to_string() << "\n";}return 0;
}

结果:

IEEE 754 single precision example
3.1415927(0x40490fdb)
sign: 0
exponent: 10000000
significand: 10010010000111111011011IEEE 754 double precision example
3.141592653589793(0x400921fb54442d18)
sign: 0
exponent: 10000000000
significand: 1001001000011111101101010100010001000010110100011000IEEE 754 half precision example
significand: 584
3.140625(0x4248)
sign: 0
exponent: 10000
significand: 1001001000

4. 结论和讨论

• 使用 format 库获得准确的 float/double 类型的打印
• 先用 bitset 获得浮点数的二进制表示， 然后根据 IEEE754 标准里的步骤， 算出精确的取值
• 获取二进制表示的时候，是偷懒做法， 是已经包含了 rounding 处理的过程； 如果打算从头算出二进制表示， 需要对整数和小数部分分别处理， 并手动 rounding。

5. References

• https://www.cnblogs.com/zjutzz/p/10140559.html