如何获得浮点除法的正确下限？答案

【问题标题】：How to get the correct floor of a floating-point division?如何获得浮点除法的正确下限？
【发布时间】：2022-08-22 18:11:38
【问题描述】：

我想获得两个正浮点数相除的浮点下限。特别是 I\'m 之后的最大浮点数不大于该除法下限的精确值。股息可以很大，除数可以很小，但在我的应用程序中，除法中没有溢出或下溢的风险。

如果我这样做：

quotient = floor(dividend / divisor);

我的问题是，当商大于尾数的精度时，除法的结果总是一个整数，所以 FPU 舍入它而不是取整它，因为它是舍入到最近的或- 偶数模式； floor() 也什么都不做，因为它已经输入了一个整数。由于它是四舍五入的，有时结果会大于确切的下限，这不是我所追求的。

在除法期间更改 FPU 的舍入模式将是一个解决方案，但这不是一个选项，所以除此之外，我怎样才能获得正确的地板？

您是否已经尝试过更正步骤，例如 e=fma(-75.0, b, a); if (e < 0.0) b = nextafter (b, 0.0); ？我并不是说这个特定的更正步骤总是有效的，只是为了澄清你所尝试的。
除了改变舍入模式之外没有其他解决办法。您正在处理的数字只是近似值。在 Python 中，您可以选择将计算作为整数 (237261451793987450000000000000) 进行，这将产生一个准确的答案，但这可能是对您没有的精度的断言。
要找到余数，在 Python 中使用 %: 2.3726145179398745e+29 % 75。在 C 中，使用 fmod: fmod(2.3726145179398745e+29, 75)。假设使用 IEEE-754 binary64 浮点格式，这两个都产生 58，这是 237261451793987452973306871808 模 75 的正确余数，而 237261451793987452973306871808 是“2.37261451793948”转换为“2.37261451793948”格式的结果。对于正操作数，正确实现的余数没有舍入误差。对于任何操作数，正确实现的 C 语言 fmod 没有舍入误差。
（当操作数具有不同的符号时，Python % 可能会出现舍入错误，因为它可能需要返回大于第一个操作数的结果，并将其置于浮点格式的不同指数区间中。）
你真的想要地板还是只想要剩下的？当地板不可代表时你想做什么？ floor(237261451793987452973306871808 / 75) 是 3163486023919832706310758290，但这不能用二进制 64 表示。最接近的可表示值为 3163486023919832955533393920，最接近的可表示值为 3163486023919832405777580032。所以，如果你真的想要floor，不使用扩展精度算术是不可能的。

标签： floating-point ieee-754

【解决方案1】：

我最终使用整数进行除法。以下函数仅适用于 IEC-559 浮点数或双精度数：

#include <stdint.h>
#include <math.h>

#ifdef __GNUC__
#define int_fast128 __int128
// other compilers pending
#endif

double truncdiv(double a, double b)
{
  int ae, be, re, sh, sh2;
  int_fast64_t am, bm;
  int_fast64_t rm;
  am = 9007199254740992. * frexp(a, &ae);
  bm = 9007199254740992. * frexp(b, &be);
  sh = 52 + (am < bm);  // add 1 if quotient is 1 bit short
  re = ae - be - sh;
  // Truncate the mantissa when the exponent is in range -52..0
  sh2 = re >= 0 ? 0 : -re;
  rm = re < -52 ? 0 : (((int_fast128)am << sh) / bm) >> sh2 << sh2;
  return ldexp(rm, re);
}

请注意，此函数不是为处理有符号零、NaN、无穷大、溢出或被零除而编写的。它也是截断除法而不是下除法，即它向零舍入，而不是向负无穷大。它需要 128 位整数类型，可能并非在所有平台上都可用。对于单精度，它只需要一个 64 位整数类型，它得到更广泛的支持：

#include <stdint.h>
#include <math.h>

float truncdivf(float a, float b)
{
  int ae, be, re, sh, sh2;
  int_fast32_t am, bm;
  int_fast32_t rm;
  am = 16777216.f * frexpf(a, &ae);
  bm = 16777216.f * frexpf(b, &be);
  sh = 23 + (am < bm);  // add 1 if quotient is 1 bit short
  re = ae - be - sh;
  // Truncate the mantissa when the exponent is in range -23..0
  sh2 = re >= 0 ? 0 : -re;
  rm = re < -23 ? 0 : (((int_fast64_t)am << sh) / bm) >> sh2 << sh2;
  return ldexpf(rm, re);
}

【讨论】：