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

Question

假设您有一个具有任意值的 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.

到目前为止我想到了什么：

我们知道给定浮点类型中正正常值（不包括零）的计数并且我们知道负值space是对称到正一个。

我们也知道离散值在取值范围内（远离零）的分布遵循一个指数函数及其相对epsilon一个相关的阶跃函数:

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

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

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

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

Answer 1

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

因此，对于给定的浮点数，只需提取尾数的最低有效位（宽度取决于向下转换类型）并按当前指数缩放。这应该（非常准确地）是 "rounding error" 将在向下转换中发生。

编辑

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

(extreme/corner个案例：例子：五位尾数为向下转换型)

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

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

Answer 2

类似以下的方法可能有效：

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

这使用了一个技巧（我相信是 Dekker 的功劳）将一个浮点数拆分为一个 "big half" 和一个 "little half"，它们的总和正好等于原始数字。我希望 "big half" 有 23 位而 "little half" 有其余部分，所以我使用常量 1 + 2^(52-23).

进行拆分

注意事项：您需要通过检查上限和下限来处理更有限的指数范围。次正规（尤其是结果在 little 类型中而不是 big 类型中次正规的情况）需要不同的特殊处理。我写 !(fabs(z/x) < releps) 而不是 fabs(z/x <= releps 是因为我希望 NaN 符合条件，因为 "bad." releps 是该变量的错误名称，因为阈值实际上比使用舍入到最近时指定的数字。

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

Accurately predicting rounding error of cast between arbitrary floating-point formats

c

algorithm

math

floating-point

floating-point-conversion