双精度算术的 std::numeric_limits

Question

我正在为C++11设计一些双双精度和四双精度算术库。部分代码基于 https://www.davidhbailey.com/dhbsoftware/ 的 QD 和 DDFUN 库。

struct Float64x2 {
    double hi;
    double lo;
};
struct Float64x4 {
    double val[4];
};

Dekker-float 类型具有常规浮点类型所没有的一些属性。例如，要能够准确地表示

FLT_MAX + FLT_DENORM_MIN

，通常需要 ~2048 位的精度来表示 (

2^1023 + 2^-1074

)。由于 15 位指数较大，Float80x2 需要 ~32768 位 （图像）。

我使用 QD 中的一些

std::numeric_limits

值作为起点。尽管我不确定它们是否是最佳使用值。特别是

epsilon()

、

min()

和

max_digits10()

。

不确定是否最好在一篇帖子中询问所有这些问题，或者作为 3 个单独的帖子来询问，所以我现在只关注其中一篇。

std::numeric_limits::epsilon()

epsilon()

是 1.0 和下一个可表示数字之间的差。对于

double

，epsilon 是

0x1.0p-52

，但我不确定对于 Dekker-float 类型它应该是什么。

1. 从字面上理解

   epsilon()

   的定义，

   numeric_limits<double>::min()

   和

   numeric_limits<double>::denorm_min()

   都应该是

   Float64x2

   或

   Float64x4

   。因为它们都可以准确地表示

   1.0 + FLT_MIN

   和

   1.0 + FLT_DENORM_MIN

   。
来自QD的

2. dd_real

   和

   qd_real

   分别将其定义为

   0x1.0p-104

   和

   0x1.0p-208

   ；这可能是

   epsilon()

   更有用的定义。此定义可用于确定何时停止向近似值添加项（即使用泰勒级数来近似

   sin(x)

   ）。
3. DDFUN 对 epsilon 使用了稍微不同的定义（来自

   ddfune.f90 dderfr

   ），即

   ldexp(1.0, -numeric_limits<T>::digits)

   。这为

   0x1.0p-53

   提供了稍小的

   double

   epsilon，并分别为

   0x1.0p-106

   和

   0x1.0p-212

   提供了

   Float64x2

   和

   Float64x4

   值。

我应该通过

epsilon()

为

std::numeric_limits

和

Float64x2

提供

Float64x4

的哪个定义？

---

此外，

std::numeric_limits::<Float64x4>epsilon()

应该返回什么？它应该返回一个

double

（它可以准确地表示 epsilon 的值），还是应该返回相同的类型（

Float64x4

）？

Answer 1

1. 从字面上理解

   epsilon()

   的定义，

   numeric_limits<double>::min()

   和

   numeric_limits<double>::denorm_min()

   都应该是

   Float64x2

   或

   Float64x4

   。因为它们都可以准确地表示

   1.0 + FLT_MIN

   和

   1.0 + FLT_DENORM_MIN

   。

…

我应该通过

epsilon()

为

std::numeric_limits

和

Float64x2

提供

Float64x4

的哪个定义？

从字面上理解

epsilon()

的定义，你根本不应该定义它。 C++ 2020 草案 N4849 17.3.5.1 [numeric.limits.members] 将其指定为：

25 机器 epsilon：1 与可表示的大于 1 的最小值之间的差。

26 对所有浮点类型都有意义。

这告诉我们，

epsilon()

对于

Float64x2

和

Float64x4

没有意义，因为它们不是浮点类型。它们是数字类型，因为它们可以表示数字，但不符合浮点模型。浮点格式表示一个数字 x、一个整数 M（或等效的定点值）和一个指数 e，其中：

        x = M•β^e−p+1

与|M| < β^p 以及整数 β 和 p 的特定值以及 e 上的界限，它们是格式的固定特征。这是来自 Muller 等人的浮点算术手册。 C++ 标准没有明确声明浮点模型，而是从 C 标准继承它，后者声明了等效模型（使用定点 M 而不是整数）。浮点运算的 IEEE 754 标准使用等效模型。

Float64x2

和

Float64x4

不符合此模型。不存在整数 β 和 p 以及 e 上的界限，因此模型中可表示的数字集合就是

Float64x2

中可表示的数字集合，对于

Float64x4

也是如此。一般来说，“浮点”不是实数或分数或非整数数字格式的同义词。它是特定类型数学模型的术语。

通过两个浮点数的未计算和来表示数字是一种不同的模型，并且不是浮点格式，因为它被浮点算术从业者普遍理解，或者如 C 标准中所述并由 C++ 继承。

来自QD的

2. dd_real

   和

   qd_real

   分别将其定义为

   0x1.0p-104

   和

   0x1.0p-208

   ；这可能是

   epsilon()

   更有用的定义。此定义可用于确定何时停止向近似值添加项（即使用泰勒级数来近似

   sin(x)

   ）。

有用并不构成为不同格式的功能赋予相同名称的数字的理由。要告知客户您的格式属性，只需定义新名称即可。