How to get the correct floor of a floating-point division?

user_2428487 0
Viewed 84

I want to obtain the floating-point floor of the division of two positive floating-point numbers. In particular I'm after the greatest floating-point number not greater than the exact value of the floor of the division. The dividend can be big and the divisor small, but in my application there's no risk of overflow or underflow in the division.

If I do this:

quotient = floor(dividend / divisor);

I have the problem that, when the quotient is greater than the precision of the mantissa, the result of the division is always an integer, so the FPU rounds it rather than flooring it because it's in round-to-nearest-or-even mode; also floor() does nothing because it's fed an integer already. Since it's rounded, sometimes the result is greater than the exact floor, and that's not what I'm after.

Changing the FPU's rounding mode during the division would be a solution, but that is not an option, so barring that, how can I obtain the correct floor?

(Related: How to correctly floor the floating point pair sum)

ieee-754 floating-point c
1 Answers

I ended up doing the division using integers. The functions below are only suitable for IEC-559 floats or doubles:

#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);
}

Note that this function is not written to handle signed zero, NaNs, infinities, denormals, overflow or division by zero. It is also truncation division rather than floor division, i.e. it rounds towards zero, not towards minus infinity. It requires a 128-bit integer type, which may not be available in all platforms. For single precision it would require only a 64-bit integer type, which is more widely supported:

#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);
}
user_24284870
Related