双精度和高精度浮点数测试（C++）

                     

贡献者： addis

　　 double 和 “双精度和变精度浮点数测试（Matlab）” 中的结果一样．

　　 在 x86 机器上，long double 一般是 1 bit 指数 + 15 bit 指数 + 64 bit 小数 = 80 bit，十进制大约是 18-19 位有效数字¹．虽然 sizeof 是 16 字节，但实际上只使用了 10 字节²．long double 的运算是硬件支持的，所以比较快．

sizeof(Doub) = 8
sizeof(Ldoub) = 16

std::numeric_limits<Doub>::max() = 1.797693135e+308
std::numeric_limits<Ldoub>::max() = 1.189731495e+4932

std::numeric_limits<Doub>::min() = 2.225073859e-308
std::numeric_limits<Ldoub>::min() = 3.362103143e-4932

std::numeric_limits<Doub>::denorm_min() = 4.940656458e-324
std::numeric_limits<Ldoub>::denorm_min() = 3.645199532e-4951

std::numeric_limits<Doub>::epsilon() = 2.220446049e-16
std::numeric_limits<Ldoub>::epsilon() = 1.084202172e-19 // 2^(-63)

std::numeric_limits<Doub>::digits10 = 15
std::numeric_limits<Ldoub>::digits10 = 18

　　 然而 sqrt 函数似乎并不会提高精度（仍然是 16 位）：

sqrt(2.)       = 1.4142135623730951454
sqrt(Ldoub(2)) = 1.4142135623730951454
precise val    = 1.4142135623730950488...

　　 另外测试一下双精度的二进制表示（图 3 ），例如双精度的 10 表示为

01000000 00100100 (48 个 0)

第一位 0 说明是正数，1 是复数．接下来 11 bit 先看成 unsigned，然后再减 $2^{10}-1 = 1023$，这里指数 10000000010 表示 $1026-1023 = 3$，于是指数项为 $2^{3} = 8$．小数项应该是 $10/8 = 1+1/4$．第 13 bit 表示 1/2，14 bit 表示 1/4，所以第 14 bit 是 1，后面的小数全是零．另外，由于 intel 的 x86 和 x86-64 都使用 little endian，所以在内存中实际上 byte 的顺序是相反的．

　　 理论上指数部分可以表示 $-1023$ 到 $1024$ 之间的数，但首尾两个数有特殊用途，所以只能表示 $-1022$ 到 $1023$．当指数为 $-1024$ 时，仍然相当于指数为 $-1023$ 但小数部分前面不加 1．当指数为 $1024$ 时，表示 nan 或者 inf．

　　 另外最大和最小的正 double 为

std::numeric_limits<double>::max() = 1.797693135e+308
01111111 11101111 11111111 11111111 11111111 11111111 11111111 11111111
std::numeric_limits<double>::min() = 2.225073859e-308 = 2^-1022
00000000 00010000 00000000 00000000 00000000 00000000 00000000 00000000
std::numeric_limits<double>::epsilon() = 2.22e-16 = 2^-52
00111100 10110000 00000000 00000000 00000000 00000000 00000000 00000000
std::numeric_limits<double>::denorm_min() = 4.940656458e-324 = 2^-1074
00000000 00000000 00000000 00000000 00000000 00000000 00000000 00000001
std::numeric_limits<double>::denorm_min()*11
00000000 00000000 00000000 00000000 00000000 00000000 00000000 00001011
double 1/0.:
01111111 11110000 00000000 00000000 00000000 00000000 00000000 00000000
double max()*2 (正无穷， 同 1/0.):
01111111 11110000 00000000 00000000 00000000 00000000 00000000 00000000
double -max()*2 (负无穷):
11111111 11110000 00000000 00000000 00000000 00000000 00000000 00000000
double sqrt(-1.) (NaN):
11111111 11111000 00000000 00000000 00000000 00000000 00000000 00000000

其中 denorm_min() 比较特殊，当指数 bit 全为 0 的时候，指数为 -1022 而不是 -1023，且小数部分前面不加 1．这叫做 denominize．

1. 四精度

　　³使用四精度比使用任意精度类型（如 MPFR 或 Flint 库）要快得多，因为任意精度浮点数是动态分配内存的．

　　 如果想完全利用 16 字节，见 GCC/g++ 编译器的 __float128 类型（文档）．相当于 FORTRAN 的 QUAD real．有 112 bit 小数，15 bit 指数和 1 bit 正负号（图 1 ）．实测的 bit 表示和上面二进制的同理（也存在两个奇怪的地方）．

　　 加上头文件 #include <quadmath.h>，使用 g++ 编译，编译时加上选项 -lquadmath 和 -fext-numeric-literals 即可．

　　 __complex128 是对应的四精度复数类型，相当于两个 __float128．

FLT128_MAX = 1.1897e+4932
FLT128_MIN = 3.3621e-4932
FLT128_EPSILON = 1.9259e-34 // 2e-112
FLT128_DENORM_MIN = 6.4752e-4966
FLT128_MANT_DIG = 113 // 包括正负位
FLT128_MAX_EXP = 16384 // 2e14
FLT128_DIG = 33 // 33-34 位十进制有效数字
M_PIq = 3.141592653589793238462643383279503 // 精确到最后一位
sqrtq(3) = 1.732050807568877293527446341505872 // 精确到最后一位

　　 给 __complex128 的实部和虚部分别赋值：

__complex128 y;
__real__ y = 实部;
__imag__ y = 虚部;

看来还是用 std::complex<__real128> 比较好．

2. 其他

　　 另见 Boost 的 float128 页面．以及 Intel 编译器的 _Quad 类型（其实也同样定义了 __float128 和 __complex128，其他功能和 g++ 几乎完全兼容）．另外一个比较有名的库是 QD．

3. Lapack

　　 最重要的 BLAS 和 Lapack 的解决方案：

• MPACK（最出名）
• mplapack（C++ 被 Julia 采用）
• Eigen 支持自定义标量类型的矩阵，但需要定义 Eigen::NumTraits<T>，详见这里．
• 可以改一下 Lapack++
• latl（没人气）

---

1. ^ Visual Studio 中据说 long double 和 double 是一样的．
2. ^ 详见这里 的 x86 extended precision format 一节
3. ^ 参考 Wikipedia 相关页面．