给定三个“双精度”值“x”、“y”和“z”，如何最小化计算“x * y / z”的精度损失？

精度损失不必绝对最小

对于您给出的示例，精度的最大损失为 1 位。

---

x=0.1 y=0.1 z=0.1

无论你采取哪种方法，你最多只能偏离一点。

---

x=23 y=13 z=23

无论你采取哪种方法，你最多只能偏离一点。

---

对于加法和减法，您希望使用类似量级的数字。但是乘法和除法没有这样的限制。您只需要注意溢流和下溢。

如果您只关心精度损失（不关心上溢/下溢），则可以使用融合乘法 （） 以更高的精度计算乘积，并手动对其进行长除法：fma

#include <math.h>
double mul_div(double x, double y, double z)
{
    double xy = x*y;    // higher bits of product
    double xy_ = fma(x, y, -xy); // remaining bits of product

    double q = xy / z;  // higher bits of quotient

    double r= xy_ + fma(-q, z, xy); // remainder of division, added to lower part of product
    q += r / z;    // add quotient of remainder
    return q;
}

如果您还必须处理上溢/下溢，您可以首先使用 ，对分数进行处理并加/减指数，然后使用 将数字重新组合为指数。frexpmul_divldexp

有一些特殊情况：

• fabs(x) == fabs(z)，确定符号并返回完成或 。y-y

• fabs(y) == fabs(z)，确定符号并返回完成或 。x-x

• 有效位和的已用位之和（或总位减去尾随的零位）小于或等于有效位的位数;因此，乘法不会导致精度损失，应该先做（因为先做除法可能会导致乘法中可避免的精度损失）。例如，1.5 是 2 位（二进制为 1.1），1.0625 是 4 位（二进制为 1.001），2+4 = 6 位，这将适合 53 位有效，没有任何精度损失。xy

• x或是 的倍数 ;首先进行相应的除法可以减少乘法中的精度损失。yz

除了这些情况（可能代替这些情况）之外，最简单的选择是使用更大的东西进行计算 - 要么解压缩成“符号、有效、指数”整数，然后自己用整数进行数学运算，或者 80 位扩展精度（也许如果你幸运的话）或 ，或成对的 .long double__float128double

如何将计算的精度损失降到最低，其中，并且都是类型？a * b / cabcdouble

无需求助于更高精度的数学运算

乘积和商，当与无限数学不完全相同时，每一步都会更高或更低 - 取决于舍入模式。*/

将普通回合与最接近，与偶数模式并列，结果将是最佳答案≤ 0.5 ULP。

给定两个操作，则最大预期误差在 1.5 ULP 以内。@njuffa。

当（中间）结果接近 2 的幂时，可能会出现一个极端情况，即 2 个步骤（然后或然后）将导致 2.0 ULP。*//*

使用 either 或 ，略微偏向于第一个。当 in 中的有效数字很少时，乘法是精确的，这与许多除法不同，因此代码只有 1 舍入的风险。mul_div_1()mul_div_2()a, b

具有更高精度的数学运算

当精度高于：long doubledouble

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

测试 2 条路径中的一条是否没有错误

虽然不是 100% 保证测试来检测除法是否准确，但它很接近。fmod(y,z) == 0

#include <math.h>
double mul_div_3(double x, double y, double z) {
  if (fmod(y,z) == 0) {
    puts("x * (y / z)");
    return x * (y / z);
  }
  puts("x * y / z");
  return x * y / z;
}

#include <float.h>
#include <stdio.h>
int main() {
  printf("%.*g\n", DBL_DECIMAL_DIG, mul_div_3(0.1, 0.1, 0.1));
  printf("%.*g\n", DBL_DECIMAL_DIG, mul_div_3(23, 13, 23));
}

输出

x * (y / z)
0.10000000000000001
x * y / z
13