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

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，而是一种浮点类型的任意精度（8bit，16bit，32bit，甚至可能是24bit），由它的位计数和指数长度指定（并遵循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。
到目前为止我的想法是：
我们知道在给定的浮点类型中正正规值（不包括零）的计数，并且我们知道负正规值空间与正正规值空间是对称的。
我们还知道，离散值在值范围（远离零）内的分布遵循指数函数，其相对epsilon遵循相关阶跃函数：
应该可以计算给定浮点类型的正常值范围内的给定实值会落在哪个离散正常值上（通过某种对数投影，或其他方式？），不是吗？考虑到这一点，我们应该能够从阶跃函数中计算出相应值的epsilon，并将其与指定的最大误差进行比较，不是吗？
我觉得这实际上可能足以计算（或至少准确估计）预期的铸造误差。我只是不知道怎么把这些东西放在一起。
你会怎么处理？（实际代码的加分：P）
注：为了提供更多的上下文：我正在研究一个 nth实现，为了找出给定值的最小有损（或给定epsilon内的有损）可转换表示，我目前正在使用上述简单的往返逻辑执行二进制搜索，以找到正确的大小。它是有效的，但缺乏效率和冷静的部门。尽管这绝不是性能瓶颈（yada yada premature optimization yada yada），但我很好奇是否能找到一个更符合数学基础和更优雅的解决方案。；）

Answer 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，因为我希望nans限定为“bad”。 releps是该变量的坏名称，因为阈值实际上比使用“舍入到最近”时指定的数字大半个ulp。