如何最小化double类型x*y/z运算的精度损失？

现有三个double类型值x、y、z，需要计算x*y/z的结果，目前有两种常见实现方式：

double mul_div_1(double x, double y, double z) {
    return x * y / z;
}

double mul_div_2(double x, double y, double z) {
    return x * (y / z);
}

但这两种实现无法始终给出最精确的结果，比如：

• 当x = 0.1，y = 0.1，z = 0.1时：
  • mul_div_1的结果为0.10000000000000002
  • mul_div_2的结果为0.1
• 当x = 23.0，y = 13.0，z = 23.0时：
  • mul_div_1的结果为13.0
  • mul_div_2的结果为12.999999999999998

已知isfinite(x) && isfinite(y) && isfinite(z) && z != 0.0，我们可以从x*y/z、(x/z)*y、x*(y/z)三种运算顺序中选出精度损失最小的结果，具体方法如下：

核心思路 #

真实值true_val = x*y/z满足true_val * z = x*y，因此可以通过残差校验判断每个计算结果的精度：计算每个结果res对应的|res * z - x*y|，残差越小，说明res越接近真实值，精度损失越小。

实现代码 #

#include <math.h>

double mul_div_opt(double x, double y, double z) {
    // 计算三种运算顺序的结果
    double res1 = x * y / z;
    double res2 = x * (y / z);
    double res3 = (x / z) * y;
    
    // 计算每个结果的残差
    double err1 = fabs(res1 * z - x * y);
    double err2 = fabs(res2 * z - x * y);
    double err3 = fabs(res3 * z - x * y);
    
    // 选择残差最小的结果
    double min_err = err1;
    double best_res = res1;
    if (err2 < min_err) {
        min_err = err2;
        best_res = res2;
    }
    if (err3 < min_err) {
        min_err = err3;
        best_res = res3;
    }
    return best_res;
}

验证示例 #

• 对于x=0.1,y=0.1,z=0.1：
  res2的残差为0，是最优结果，符合预期；
• 对于x=23.0,y=13.0,z=23.0：
  res1的残差为0，是最优结果，符合预期。

这种方法通过简单的残差比较，就能在三种运算方式中选出精度最高的结果，有效最小化浮点计算的精度损失。

内容的提问来源于stack exchange，提问作者EFanZh