10

In using double fmod(double x, double y) and y is an integer, the result appears to be always exact.

(That is y a whole exact number, not meaning int here.)

Maybe C does not require fmod() to provide an exact answers in these select cases, but on compilers I've tried, the result is exact, even when the quotient of x/y is not exactly representable.

  1. Are exact answers expected when y is an integer?
  2. If not, please supply a counter example.

Examples: 

double x = 1e10;
// x = 10000000000
printf("%.50g\n", fmod(x, 100));
// prints 0

x = 1e60;
// x = 999999999999999949387135297074018866963645011013410073083904
printf("%.50g\n", fmod(x, 100));
// prints 4

x = DBL_MAX;
// x = 179769313486231570...6184124858368
printf("%.50g\n", fmod(x, 100));
// prints 68

x = 123400000000.0 / 9999;
// x = 12341234.1234123408794403076171875
printf("%.50g %a\n", fmod(x, 100), fmod(x, 100));
// prints 34.1234123408794403076171875 0x1.10fcbf9cp+5

Notes:
My double appears to the IEEE 754 binary64 compliant.
The limitations of printf() are not at issue here, just fmod().

[Edit]

Note: By "Are exact answers expected", I was asking if the the fmod() result and the mathematical result are exactly the same.

chux - Reinstate Monica
  • 143,097
  • 13
  • 135
  • 256

2 Answers2

13

The result of fmod is always exact; whether y is an integer is irrelevant. Of course, if x and/or y are already approximations of some real numbers a and b, then fmod(x,y) is unlikely to be exactly equal to a mod b.

R.. GitHub STOP HELPING ICE
  • 208,859
  • 35
  • 376
  • 711
11

The IEEE Standard 754 defines the remainder operation x REM y as the mathematical operation x - (round(x/y)*y). The result is exact by definition, even when the intermediate operations x/y, round(x/y), etc. have inexact representations.

As pointed out by aka.nice, the definition above matches the library function remainder in libm. fmod is defined in a different way, requiring that the result has the same sign as x. However, since the difference between fmod and remainder is either 0 or y, I believe that this still explains why the result is exact.

Emilio Silva
  • 1,942
  • 14
  • 17
  • +1 for the citation. I was looking for a quick intuitive explanation of why there's always an exact result for `x - (round(x/y)*y)` within the original precision, but it was getting too complicated so I left it out of my answer. – R.. GitHub STOP HELPING ICE Jan 05 '14 at 01:58
  • 2
    I don't think there is an easy explanation. Naively calculating `x-round(x/y)*y` is not identical to `fmod(x, y)` for some extreme cases that I tried (`x = pow(2, 53)*100` and `y = 100`, for instance). Implementations have to perform tricks like iterative subtraction to keep the result mathematically exact. – Emilio Silva Jan 05 '14 at 02:02
  • The best explanation I know starts by subtracting `r^n * y` where `r` is the radix and `n` is the exponent that yields a result whose exponent is the same as the exponent of `x`; this reduces the number of places of precision by at least one and yields a value congruent to `x` mod `y`. I suspect you apply this argument inductively to get the conclusion, but the details (especially when sign flips) are a pain. – R.. GitHub STOP HELPING ICE Jan 05 '14 at 02:14
  • 3
    @R.. What about the following explanation: let `r` be the mathematical result. Both `x` and `y` are floating-point numbers larger than `r` and therefore multiples of `ulp(r)`. Therefore `r` is a multiple of `ulp(r)`, therefore `r` is exactly representable as a floating-point number. – Pascal Cuoq Jan 05 '14 at 09:39
  • @PascalCuoq: Not all exact multiples of `ulp(r)` are representable. For instance, `LLONG_MAX*ulp(r)` most certainly is not. – R.. GitHub STOP HELPING ICE Jan 05 '14 at 19:25
  • @R.. I should have added “and `r` is smaller than `x`, which we have already assumed to be a finite floating-point number when we started computing a finite result for `fmod(x, y)`”. Also “smaller than” and “larger than” must be understood as “in absolute value”. – Pascal Cuoq Jan 05 '14 at 20:12
  • @PascalCuoq: I still don't get your argument. You're starting with a mathematical result `r` then reasoning with `ulp(r)`, which is not even well-defined for a mathematical result that doesn't already have a floating-point type associated with it. So as far as I can tell your reasoning is circular. – R.. GitHub STOP HELPING ICE Jan 05 '14 at 20:19
  • @R.. It is common to define ULP(x) for x real as “the distance between the two closest straddling floating-point numbers a and b”, to the point that this previous quote is taken from Wikipedia. – Pascal Cuoq Jan 05 '14 at 20:33
  • Why round? Isn't it rather truncate? x-trunc(x/y)*y – aka.nice Jan 07 '14 at 08:51
  • @aka.nice, the standard says precisely "the integer nearest to the exact value x/y", which I believe is `round(x/y)`. I agree that it looks strange. – Emilio Silva Jan 08 '14 at 00:29
  • @EmilioSilva In such case fmod(0.75,1.0) would be -0.25 right? What you describe matches more the function remainder(x,y) of libm. – aka.nice Jan 08 '14 at 00:37
  • @aka.nice, that is correct. The definition of `remainder(3)` matches the IEEE754 operator `REM`. `fmod` is defined slightly different, requiring that the result has the same sign as `x`, which is of course `x-trunc(x/y)*y`. – Emilio Silva Jan 08 '14 at 00:41