准确预测任意浮点格式之间转换的舍入误差

Posted 2023-02-16

【中文标题】准确预测任意浮点格式之间转换的舍入误差

【英文标题】：Accurately predicting rounding error of cast between arbitrary floating-point formats

【发布时间】：2015-04-15 11:09:07

【问题描述】：

假设您有一个具有任意值的 float64_t 数字，并且您想知道该数字是否可以安全地向下转换为 float32_t，但限制是生成的舍入误差不得超过给定的 epsilon .

可能的实现如下所示：

float64_t before = 1.234567890123456789;
float64_t epsilon = 0.000000001;
float32_t mid = (float32_t)before;  // 1.2345678806304931640625
double after = (float64_t)mid; // 1.2345678806304931640625
double error = fabs(before - after); // 0.000000009492963526369635474111
bool success = error <= epsilon; // false

---

为了让事情更有趣，我们假设您不应该对这两种类型之间的手头值执行任何实际类型转换（如上所示）。

只是为了提高一个档次：假设您不是向下转换为float32_t，而是任意的浮点类型精度（8 位、16 位、32 位，甚至可能是 24 位）由其位数和指数长度指定（并遵循 IEEE 754 的约定，例如舍入到偶数）。

所以我正在寻找的是更类似于此的通用算法：

float64_t value = 1.234567890123456789;
float64_t epsilon = 0.000000001;
int bits = 16;
int exponent = 5;
bool success = here_be_dragons(value, epsilon, bits, exponent); // false

举个例子，将 64 位数字 1.234567890123456789 向下转换为较低的精度会导致以下舍入误差：

 8bit: 0.015432109876543309567864525889
16bit: 0.000192890123456690432135474111
24bit: 0.000005474134355809567864525889
32bit: 0.000000009492963526369635474111
40bit: 0.000000000179737780214850317861
48bit: 0.000000000001476818667356383230
56bit: 0.000000000000001110223024625157

---

已知情况：

所讨论的两种精度类型（一种精度低于另一种）的规范： 总长度（以位为单位）（例如，浮点数为 32） 指数长度（以位为单位）（例如，浮点数为 8） 每种类型的 min 和 max 值（因为这些值可以从上面派生）。 正正常值（不包括零）的数量 (((2^exponent) - 2) * (2^mantissa)) 指数的bias ((2^(exponent - 1)) - 1) 实际的value（在给定的更高精度类型中提供）。 错误阈值epsilon 允许向下转换以被视为成功（也在给定的更高精度类型中提供）。

（预期误差的近似值可能就足够了，这取决于它的准确性和偏差因素。但显然更喜欢精确的计算。）

不需要涵盖的案例（因为它们很容易单独解决）：

如果输入值是任何

非正态值

（次正态、无穷大、nan、零……），则答案应在此定义为true。 如果输入值落在

已知边界之外

（+- 给定的 epsilon）的给定类型的较低精度，那么答案应在此定义为false。

---

到目前为止我的想法：

我们知道给定浮点类型中正正常值（不包括零）的计数，并且我们知道 负 值空间与 对称 >正面一个。

我们还知道，离散值在取值范围内（远离零）的分布遵循指数函数，其相对ε是相关阶跃函数：： p>

应该可以计算出nth 离散正常值一个给定的真实值在给定浮点类型的正常值范围内将落在上（通过某种对数投影或其他方式？），不是吗？鉴于此n，然后应该能够从其step function计算相应值的epsilon并将其与指定的最大误差进行比较，不是吗？

我觉得这实际上可能足以计算（或至少准确估计）预期的铸造误差。我只是不知道如何将这些东西放在一起。

你会如何处理这个问题？ （实际代码加分：P）

_{Ps：为了提供更多背景信息：我正在研究 var_float 实现，并且为了找出给定值的最小无损（或在给定 epsilon 内有损）可转换表示，我' m 目前使用上述简单的往返逻辑执行二进制搜索以找到正确的大小。它有效，但在效率和冷静部门缺乏。尽管它绝不是性能瓶颈（yada yada 过早优化 yada yada），但我很好奇是否可以找到更基于数学和优雅的解决方案。 ;)}

【问题讨论】：

你不能只看一下 8 字节双精度的最后一位是否适当数量的零吗？还要检查指数。较小格式的范围更小。

那行不通——我认为——因为“0.328125”可以无损地转换为float8_t（0 001 0101），不留尾随零。我的“不需要涵盖”规则涵盖了范围问题。 ;)

在较长的格式中，21/64=(1.0101)_2*2^(-2) 的尾数以零继续。这正是将其压缩成较短格式所需的条件。

啊，我看错了你的评论。以为您是在专门谈论 8 位表示中的尾随位。现在说得通了（有点像 DrKoch 的回答，不是吗？）。

是的，这就是为什么我不添加自己的答案。但是，这种方法只适用于真正的无损压缩。如果允许一些有限的错误，则需要更多的努力。

【参考方案1】：

类似以下的方法可能会起作用：

double isbad(double x, double releps) 
  double y = x * (1 + 0x1.0p29);
  double z = y-x-y+x;
  return !(fabs(z/x) < releps);

这使用了一个技巧（我相信由于 Dekker）将浮点数拆分为“大半”和“小半”，它们的总和与原始数字完全一致。我希望“大半”有 23 位，而“小半”有其余的，所以我使用常数 1 + 2^(52-23) 进行分割。

注意事项：您需要通过检查上限和下限来处理更有限的指数范围。次正规（特别是在小类型中结果为次正规但不是大类型的情况）需要不同的特殊处理。我写了!(fabs(z/x) < releps) 而不是fabs(z/x <= releps，因为我希望NaN 被认为是“坏的”。 releps 是该变量的错误名称，因为阈值实际上比您在使用四舍五入时指定的数字大半个 ulp。

【讨论】：

哇，就是这样！将函数的最后一行更改为 return fabs(z); 使其返回确切的错误。碉堡了。这就是我所希望的那种魔法。

此页面有两个参考文献，前一个是 Veltkamp (1968)：proval.lri.fr/gallery/Dekker.en.html

@PascalCuoq：你有 Veltkamp 论文的电子版吗，或者知道在哪里可以找到？

@tmyklebu 我看到她会问网页的作者。

【参考方案2】：

向下转换相当于将尾数的最低有效位设置为零。

因此，对于给定的浮点数，只需提取尾数的最低有效位（宽度取决于向下转换类型）并使用当前指数进行缩放。这应该（非常准确地）是向下转换时会发生的“舍入误差”。

---

编辑

如 cmets 中所述，上述情况仅适用于所有情况的 50%。 （当向下转换导致向下舍入）。在向下转换导致四舍五入的情况下，稍微修改的方法将有所帮助：

（极端/极端情况：示例：向下转换类型的五位尾数）

Rounding down: 0x1.00007fff -> 0x1.0000 
               -> Err == 0x0.00007fff

Rounding up:   0x1.00008000 -> 0x1.0001 -> Err == 0x1.00010000 - 0x1.00008000
               -> Err == 0x0.00008000

【讨论】：

“向下转换相当于将尾数的最低有效位设置为零”并不完全正确，因为从 double 到 float 的转换 rounds 到最接近的 float编译平台。您所描述的是如果转换截断会发生什么。

例如，将双精度 0x1.000001fffffffp0 舍入到 float 时的错误比您的方法所相信的要少数百万倍，因为前面提到的 double 只是简单地转换为 0x1.000002p0f。

好点，我的var_float 类型确实做到了平局。我猜应该特别提到这一点。将更新我的问题。

但是如果最后一位真的为零，那么向下转换只是截断。如果它们不为零，那么无论舍入规则如何，向下转换都会给出不同的数字。