Check if a double is evenly divisible by another double in C?

Question

How can I check if a double x is evenly divisible by another double y in C? With integers I would just use modulo, but what would be the correct/best way to do it with doubles?

I know floating point numbers carry with them imprecision, but I'm getting the double from standard input. Maybe I should not scan it as a double straight away but as two integers instead, but where would I go from then?

Answer 1

The standard header math.h defines the following functions:

• double fmod(double x, double y);
• float fmodf(float x, float y);
• long double fmodl(long double x, long double y);

These functions return the result of the remainder of x divided by y. The result has the same sign as that of x. You can use r = fmod(x, y); for double numbers x and y, and check if r == 0. If you want to not test for exact divisibility but add some tolerance, then you can check if r is "close enough" to 0 or y (thanks caf).

fmodf() and fmodl() are new in C99.

Edit: C99 also defines a separate remainder(double x, double y) function, that returns the remainder of x/y. From http://docs.sun.com/source/806-3568/ncg_lib.html:

The remainder(x,y) is the operation specified in IEEE Standard 754-1985. The difference between remainder(x,y) and fmod(x,y) is that the sign of the result returned by remainder(x,y) might not agree with the sign of either x or y, whereas fmod(x,y) always returns a result whose sign agrees with x. Both functions return exact results and do not generate inexact exceptions.

...

When y ≠ 0, the remainder r = x REM y is defined regardless of the rounding mode by the mathematical relation r = x - ny, where n is the integer nearest the exact value of x/y; whenever | n - x/y | = 1/2, then n is even. Thus, the remainder is always exact. If r = 0, its sign shall be that of x. This definition is applicable for all implementations.

(Either fmod() or remainder() should work for you.)

Answer 2

The fmod() family of functions give terrible results. Suppose you want to determine if 42 is evenly divisible by 0.4. It is, 105 times. However, fmod does the division and gets a result like 104.99999 which it then rounds down to 104 resulting in a remainder of 0.399999 which gives a false negative result. remainderl(), however, seems to work. Even 0.4 itself is represented inexactly in floating point.

For the folks who don't grok the concept of "evenly divisible", it has nothing to do with the result being an even number - you probably have your etymology backwards. Even numbers are those numbers which are evenly divisible by 2. And the concept of divisibility is entirely valid for non-integers. Evenly divisible means the result of the division is an integer regardless of whether the dividend or divisor are. An example application is if you have a metal lathe with a 3mm pitch leadscrew and are cutting a 0.4mm pitch bolt. 14 threads at 3mm line up with 105 threads at 0.4mm. The divisibility calculation is used to tell where the various moving parts of the lathe sync up again so you can reengage for the next cutting pass. Another example is imperial measurements which have been converted to metric. 50.8mm (2") is evenly divisible by 25.4mm (1"). Even without metric conversions, dimensions are often non-integers yet divisibility is often an issue: 0.5" is evenly divisible by 0.1", 0.125". and 0.250". Converting a floating point number (such as 0.375") to a fractional representation (3/8") is one more application of divisibility to non-integer numbers.

The two alternative calculations in this sample function give the same results for hundreds of different number pairs. However, replacing remainderl() with fmodl() or roundl() with floorl() gives lots of invalid results. I originally used a fuzz of 0.001. Actual calculation error seems to usually be of order 1E-15 so a smaller fuzz can be used. However, comparing the result to 0.0 will give false negative results. You might want to express your fuzz in terms of your denominator in case you are working with very small numbers. divisible(42, 0.4) and divisible(41,0.4) should give the same results as divisible(0.000000042, 0.0000000004) and divisible(0.000000041, 0.0000000004). I.E. are 42nm and 41nm divisible by 0.4nm? With the version of the function given here, they do. With a fixed fuzz, they do not necessarily. However, divisible(42, 0.0000000004) still gives a false negative (error is 1.53003e-15 which is larger than the fuzz of 4E-19) so comparing numbers that differ by 9 orders of magnitude is not reliable. IEEE floating point has its limitations. Notice I used long double calculations to minimize calculation and representation errors. This function was not tested with negative numbers.

int divisible(long double a, long double b) 
{
  int result;
#if 1
   if(fabsl(((roundl(a/b)*b)- a)) <= (1E-9*b) ) {
    result=TRUE;
  } else {
    result=FALSE;
  }
#else
  if( fabsl(remainderl(a,b)) <= (1E-9*b ) ){
    result=TRUE;
  } else {
    result=FALSE;
  }
#endif
  // printf("divisible(%Lg, %Lg): %Lg, %Lg,%d\n", a, b, roundl(a/b), fabsl(((roundl(a/b)*b)-a)), result);
  return(result);
}