双精度和四精度浮点数笔记(C++)

                     

贡献者: addis

  • 本文处于草稿阶段。

   另见 “双精度和变精度浮点数测试(Matlab)”。

double

   C++ 浮点相关函数:

ULP

   unit in the last place (ULP) 就是 mantissa 增加 1 bit 时表示的数增加多少。这等于 $2^{-52}\times 2^{e}$。$e$ 是双精度的指数(以 2 为底)。其中 $2^{-52}$ 一般记为 eps

   所以双精度的 mantissa 表示 1 时,相对精度就是 eps,表示 1-eps 时,相对精度就非常接近 eps/2最小相对误差)。所以一般来说相对精度在二者之间。但一个例外就是当绝对值非常小时(小于 $2^{-1022}$)也就是 $2.225 \times 10^{-308} $),也就是 subnormal,mantissa 的第一 bit 可能并不是默认的 1,此时相对误差可能会变得很大。

需要多少位十进制数

   需要多少位十进制数可以通过四舍五入精确表示任意的双精度数呢?$N$ 位有效数字的十进制数的相对误差介于 $10^{-N}$(前几位是 999... 时)和 $10^{-N+1}$(前几位是 1000... 时)之间。所以为了满足 eps/2($1.11 \times 10^{-16} $)的相对误差,最多需要 17 位有效数字

   那么一个 double 能不能精确地用十进制表示呢?可以的,因为 eps 可以,但可能需要非常多的位数尤其是 2 的指数非常大或者非常小的时候可能会需要几百位,没有太大必要。

long double

   在 x86 机器上,long double 一般是 1 bit 符号 + 15 bit 指数 + 64 bit 小数 = 80 bit,十进制大约是 18-19 位有效数字2。虽然 sizeof 是 16 字节,但实际上只使用了 10 字节3long 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,小数项应该是 $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$。当指数为 $-1023$ 时,仍然相当于指数为 $-1022$ 但小数部分前面不加 1(如果小数也都是零,那就表示 0,所以 0 除了第一 bit 其他都必须是 0)。当指数为 $1024$ 时,表示 inf(要求小数部分都是零,正负号仍然有效)或者 nan0b1... 表示安静的 qNaN0b0..1.. 表示会发错误信号的 sNaN)。

   双精度类型可以表示 -2e532e53 的每一个整数,更小或更大的整数有的也可以表示(事实上此时只可以表示整数),只是不连续。

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. 四精度

  4使用四精度比使用任意精度类型(如 MPFR 或 Flint 库)要快得多,因为任意精度浮点数是动态分配内存的。

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

图
图 1:见 Wikipedia

   加上头文件 #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> 比较好。

   另见 “C++ 高精度计算资源列表”。


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


致读者: 小时百科一直以来坚持所有内容免费无广告,这导致我们处于严重的亏损状态。 长此以往很可能会最终导致我们不得不选择大量广告以及内容付费等。 因此,我们请求广大读者热心打赏 ,使网站得以健康发展。 如果看到这条信息的每位读者能慷慨打赏 20 元,我们一周就能脱离亏损, 并在接下来的一年里向所有读者继续免费提供优质内容。 但遗憾的是只有不到 1% 的读者愿意捐款, 他们的付出帮助了 99% 的读者免费获取知识, 我们在此表示感谢。

                     

友情链接: 超理论坛 | ©小时科技 保留一切权利