What shall I do, if C compares two doubles wrong?

Question

I want to write a fmod() function

double fmod(double x, double y) {
    double mod = x;

    while(mod >= y)
    {
        mod -= y;
    }
    
    return mod;
}

But fmod(1.2, 0.05) returns 0.05

Answer 1

Although the title asks about incorrect comparison, the comparison in the program shown is correct. It is the only correct floating-point operation in the program; all the others have errors (compared to real-number arithmetic).

In fmod(1.2, 0.05), neither 1.2 nor 0.5 are representable in the double format used in your C implementation. These numerals in source code are rounded to the nearest representable values, 1.1999999999999999555910790149937383830547332763671875 and 0.05000000000000000277555756156289135105907917022705078125.

Then, in the subtraction in mod -= y;, the exact real arithmetic result, 1.14999999999999995281552145343084703199565410614013671875 is not representable, so it is rounded to 1.149999999999999911182158029987476766109466552734375.

Similar errors continue during the calculations, until eventually 0.0499999999999994615418330567990778945386409759521484375 is produced. At each point, the comparison mod >= y correctly evaluates whether mod is greater than or equal to y. When mod is less than y, the loop stops.

However, due to intervening errors, the result produced, 0.0499999999999994615418330567990778945386409759521484375, is not equal to the residue of 1.1999999999999999555910790149937383830547332763671875 divided by 0.05000000000000000277555756156289135105907917022705078125. The correct result can be calculated with the standard fmod function, which returns 0.04999999999999989175325509904723730869591236114501953125.

Note that, when you define a function named fmod, the C standard does not define the behavior because this conflicts with the standard library function of that name. You ought to give it a different name, such as fmodAlternate.

Inside the fmod routine, errors can be avoided. It is possible to implement fmod so that it returns an exact result for the arguments it is given. (This is possible because the result is always in a region of the floating-point range that is fine enough [has a low enough exponent] to represent the real arithmetic result exactly.) However, the errors in providing the arguments cannot be corrected: It is not possible to represent 1.2 or 0.05 in the double format your C implementation uses. The source code fmod(1.2, .05) will always calculate fmod(1.1999999999999999555910790149937383830547332763671875, 0.05000000000000000277555756156289135105907917022705078125), which is 0.04999999999999989175325509904723730869591236114501953125.

An alternative is to represent the numbers differently. For example, you could scale these numbers by a factor of 100, and fmod(120, 5) will return 0. What solution is appropriate depends on the circumstances of the problem you are trying to solve.