views:

818

answers:

7

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?

A: 

The concept of "even number" is only defined for integers. You can't apply it to doubles; it does not make mathematical sense. From Wikipedia:

An even number is an integer that is "evenly divisible" by 2, i.e., divisible by 2 without remainder.

I suggest you convert your doubles to ints, by applying whatever method you decide (rounding, truncation) and then use modulo as you suggest.

CesarGon
He wants to know if one number is evenly divisible by a second number, nothing about even numbers. 5.0 is evenly divisible by 2.5
Shane C. Mason
So is what he wants to know that one number is an exact multiple of some other number? Because that makes more sense, although I suspect it'd be lost in machine-representation-of-double land.
Calum
That is incorrect, Shane. "Evenly" is not a concept you can apply to non-integers. Wait until you bump into non-integers that can't be represented in binary and you'll see. :-)
CesarGon
@calum - that is true - but he could determine it within a reasonable tolerance.
Shane C. Mason
@downvoter: Why the downvote? An explanation would help me and others learn what's wrong with my answer. Thank you. :-)
CesarGon
@CesarGon "r1 is evenly divisible by r2" is perfectly well defined for two reals: it is true iff there exists and integer n such that r1 = n*r2. For instance 2π is evenly divisible by π/2.
Pascal Cuoq
@Pascal: *reals*, yes, but not *doubles*. Reals is a mathematical concept. Doubles is a binary-representation concept, and it suffers bit-error representation issues. Not all reals are representable as doubles.
CesarGon
@CesarGon you're grasping at straws now. it makes perfect sense, but the limitations of binary representation won't allow it. What is wrong with saying this as an answer?
San Jacinto
@San: Well, read my answer. I said just that: "I suggest you convert doubles to ints, and then use modulo as you suggest". That is similar to what Laurence suggests in his answer. I still can't see what's wrong with my answer as to deserve two downvotes. Can anyone explain please?
CesarGon
@CesarGon Perhaps the downvotes have to do with going on a tangent with the definition of "even", when the question is about "evenly divisible" and the difficulties of doing it with floating-point are precisely the topic of the question.
Pascal Cuoq
It's perfectly defined for doubles, too. A double `d1` is exactly divisible by a double `d2` iff there is an integer `n` such that `d1 = n*d2` as real numbers.
Stephen Canon
I agree that "evenly divisible" is a rather nebulous concept when applied to floating point numbers. I up voted you. At least the voting works in integer math :)
Mark Wilkins
@Pascal: Thanks for the clarification. I didn't think I was being tangential; I honestly thought I was helping. @Stephen: Can you provide any reliable backing of your claim? And what does "as real numbers" mean? @Mark: Thanks. :-)
CesarGon
Does anyone (Shane, Pascal, Stephen in particular) have a reference that actually shows that the term "evenly divisible" is defined as suggested here? The implication seems to be that the result is an integer, but a bit of googling didn't turn up any such references. Maybe I was using the wrong search terms.
Mark Wilkins
I know I'll get a few more downvotes, but hey. If you check the definitions for Even Number, Divisor and Divisibility Tests on http://mathworld.wolfram.com you'll see that *all* of them are based exclusively on integers. No divisibility tests (even or otherwise) are defined on non-integers.
CesarGon
@Mark: The phrase "evenly divisible" used in this way is admittedly a slight abuse of mathematical terminology, but a very common one, and clearly what the questioner is asking about.
Stephen Canon
@Cesar: Sure I can back up my claim -- I'm one of the authors of the IEEE-754 standard, and I say that it's true =). More on point, what I mean by real numbers is that the numerical value of `d1` (which is a real number) should be exactly equal to the product of the (real number) value of `d2` and the integral real number `n`. i.e. that `n` may not be representable as a double precision number, and that the product of `n*d2` is computed as a real number, without rounding. This corresponds to the notion of remainder used in the IEEE-754 `remainder( )` function.
Stephen Canon
... that is to say that under this notion of "evenly divisible", `d1` is evenly divisible by `d2` if and only if `remainder(d1,d2) == 0.0`.
Stephen Canon
@Stephen: Thanks for the clarification. Still, I find the mixing of doubles and reals to be a problem, for the reasons I explained above. Also, the comparison of double-typed values (as in `remainder(d1,d2) == 0.0`) is strongly discouraged; I am sure you know that. I must agree with your previous comment that the phrase "evenly divisible" used in this way is an abuse of mathematical terminology.
CesarGon
@Cesar: The mixing of doubles and reals is a pain linguistically, but unavoidable when discussing this type of behavior (this is fundamentally a question about the real number values represented by doubles). On the issue of comparing doubles for equality, I'll have to disagree; it should only be discouraged if the programmer doesn't understand what it means. There are too many situations (such as this one) where it is exactly what is needed to proscribe its use entirely.
Stephen Canon
Do you mean that programmers should write code like `remainder(d1,d2) == 0.0``? Or are you using that as an expression to illustrate what you mean only? (just checking)
CesarGon
Yes, programmers who understand what exact comparison of floating-point values actually means should absolutely use it when appropriate. Sometimes exact equality is really what you're testing for, even in floating point, and the usual approximate equality tests will result in errors.
Stephen Canon
@Cesar (much earlier comment): The existence of real numbers that don't have a finite representation as doubles does not imply that the division of two representable doubles is inexact.
Drew Hall
+3  A: 

I am not sure what you are trying to do, but I have used fmod() from math.h in audio synthesis code where I needed my parameters to be floats or doubles and I needed to get a modulo.

Justin Smith
+2  A: 
  1. Scan them in as doubles and call them x1 and x2
  2. Find what x1/x2 is using division and call it x3
  3. Find x1 - (x2*x3) and see if that number is sufficiently close to zero - if it is then x1 is evenly divisible by x2 - (obviously taking into consideration the the possibility of negative values here)

lol - line 3 fixed :)

Shane C. Mason
You may want to fix line 3.
Pascal Cuoq
Sadly, even after fixing line 3, there is still a multiplication which will be rounded to the one of the nearest floating-point numbers (the nearest if you didn't change the processor's rounding mode). Thus the result of the subtraction will appear to be zero even when the difference between x1 and the product could have been represented as a floating-point. Better use `fmod` as suggested elsewhere and hope that it does the right thing.
Pascal Cuoq
+2  A: 

If you want to be absolutely precise, you could used fixed-point math. That is, do everything with ints, but ints that are (in your case) some power of 10 of the value they actually represent.

Say the user enters 123.45 and 6789.1. First, you want to make sure you've got the same number of decimal places, so add trailing zeros to the one with fewer decimal places. That gives us 123.45 and 6789.10 (now both with 2 decimal places). Now just remove the decimal point, to get 12345 and 678910. If one divides into the other evenly, that's your answer.

This works because removing the decimal point multiplies both by the same constant (100 in the example above). (x * 100) / (y * 100) == x / y

A few things to be careful about: if you read the integer part and fractional part as ints, be careful that you don't lose leading zeros on the fractional part. (eg: 0.1 and 0.0001 are not the same number!) Also, if there are enough decimal places you can overflow. You probably want to at least use a long.

You could also do the computation with doubles, but it'll be less precise. To do it that way, do the division, and then compare the difference between the result and the rounded result. If within some small tolerance, then it divides evenly.

Laurence Gonsalves
Fixed-point math is just another distribution of a finite number of values to represent a continuous range of reals, and with floating-point, you can at least hope that `fmod` will be precise at one ulp and that the concentration around zero will be playing in your favor.
Pascal Cuoq
However, in fixed-point math the finite distribution is identical to that implied by the input format itself. If you're taking the user's input in floating-point, then using floating-point to do the calculations is appropriate; but if you're taking the user's input as a limited-length (an hence fixed-point) decimal, then they'll probably expect that it returns "yes" when asked if `0.3` is divisible by `0.1`.
caf
And fmod(0.3/0.1) is 0.0, this stuff about using fixed point is irrelevant here.
Justin Smith
@Pascal and @Justin: In floating point 2.3 / 0.1 = 22.999999999999996. In fixed point you see it's 23 / 1. You can get around this with rounding, as I mentioned, but (as caf said) the input is essentially fixed point, so you can avoid loss of precision by keeping it fixed point.
Laurence Gonsalves
Justin: Try it and see - `fmod(0.3, 0.1)` returns `0.1`.
caf
+10  A: 
Alok
+1 This is exactly why the `fmod` function exists. No tolerance is needed, however -- a correct implementation of `fmod` is always exact.
Stephen Canon
True. What I mean is that the `fmod` function determines exactly if two `double`s divide each other evenly. Depending upon the OP's needs, he may want to add an additional tolerance. I am fixing my wording now.
Alok
@Stephen, thanks for your remark. I have now updated my post.
Alok
@Alok I erased the comment you were replying to in shame (sorry)... Indeed, both `fmod` and `remainder` **can** be exact (if implemented correctly), as the difference of two close enough floating-point numbers can always be exactly represented. My previous comment was therefore of little interest.
Pascal Cuoq
@Pascal. Yeah, I just read F.9.7.1 in n1336, and realized that myself. Thanks for your comments!
Alok
You also need to check if `r` is "close enough" to `y` - for example, try `fmod(0.3, 0.1)`.
caf
@caf: yes, of course.
Alok
+2  A: 

How can I check if a double x is evenly dividable 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?

You would include and link to the math library:

#include <math.h>

Then you would call the floating point modulus function fmod:

if (fmod(5.0, 2.5) == 0.0)
  // evenly divisible
else
  // not evenly divisible

You may want to compare the result of fmod with a small value instead of 0.0 depending on your needs.

+1  A: 

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