问题
对于实现精确IEEE 754算法的C99编译器，类型f，divisor的值float是否存在f / divisor != (float)(f * (1.0 / divisor))？
编辑：通过“实现精确的ieee 754算法”，我指的是一个正确定义flt_eval_方法为0的编译器。
上下文
提供符合ieee 754标准的浮点的c编译器，如果所述逆本身可以表示为afloat，则只能用逆的单精度乘法代替常数的单精度除法。
实际上，这只发生在二次方的情况下。因此，程序员alex可能有信心将f / 2.0f编译为f * 0.5f，但如果alex可以用0.10f而不是除以10来进行乘法，alex应该通过在程序中编写乘法或使用编译器选项（如gcc的-ffast-math）来表示它。
这个问题是关于将单精度除法转换为双精度乘法。它总是产生正确的四舍五入结果吗？它是否有可能更便宜，从而成为编译器可能进行的优化（即使没有-ffast-math）？
我比较了1和2之间(float)(f * 0.10)的所有单个精度值f / 10.0f和f，但没有找到任何反例。这应该包括产生正常结果的所有正常floats的划分。
然后我用下面的程序把这个测试推广到所有除数：
#include <float.h>
#include <math.h>
#include <stdio.h>
int main(void){
for (float divisor = 1.0; divisor != 2.0; divisor = nextafterf(divisor, 2.0))
{
double factor = 1.0 / divisor; // double-precision inverse
for (float f = 1.0; f != 2.0; f = nextafterf(f, 2.0))
{
float cr = f / divisor;
float opt = f * factor; // double-precision multiplication
if (cr != opt)
printf("For divisor=%a, f=%a, f/divisor=%a but (float)(f*factor)=%a\n",
divisor, f, cr, opt);
}
}
}
搜索空间足够大，可以让这个有趣（246）。程序当前正在运行。有人能告诉我，在它完成之前，它是否会打印一些东西，也许会解释为什么或为什么不打印？
最佳答案
你的程序不会打印任何东西，假设圆形连接到偶数圆形模式。这一论点的实质如下：
我们假设f和divisor都在1.0和2.0之间。因此f = a / 2^23和divisor = b / 2^23对于a范围内的一些整数b和[2^23, 2^24)。案例divisor = 1.0并不有趣，因此我们可以进一步假设b > 2^23。(float)(f * (1.0 / divisor))可能给出错误结果的唯一方法是，精确值f / divisor非常接近中间值（即，两个单精度浮点之间正好有一个数字），以致表达式中的累积错误将我们推到中间值的另一边。
但这不可能发生。为了简单起见，我们首先假设f * (1.0 / divisor)，因此精确的商在f >= divisor中。现在，对于区间[1.0, 2.0)中的单个精度，任何中间条件的形式都是一些带[1.0, 2.0)的奇数整数c / 2^24。c的确切值是2^24 < c < 2^25，因此差异的绝对值f / divisor限定在a / b以下，因此至少是f / divisor - c / 2^24（因为1 / (2^24 b)）。因此，我们比任何中间情况下的双精度ulps都要远得多，而且应该很容易证明，双精度计算的误差永远不会超过16ulps。（我还没有做过算术，但我想很容易在错误上显示3个ulps的上限。）
因此1 / 2^48不可能足够接近产生问题的中间案例。请注意，b < 2^24也不可能是一个确切的中间情况：因为16是奇数，f / divisor和f / divisor是相对素数，所以我们得到c的唯一方法是如果c是2^24的倍数。但是c / 2^24 = a / b在b范围内，所以这是不可能的。2^24的情况是相似的：中间的情况有b的形式，类似的论点表明(2^23, 2^24)大于f < divisor，这再次给了我们一个c / 2^25双精度ulps的余地。