(define (fmod numerator denominator)
(- numerator
(* (floor (/ numerator denominator))
denominator)))
(define (test-fmod numerator denominator)
(let ((answer (fmod numerator denominator)))
(if (> answer denominator)
(begin (display "Whoops: ")
(display (list numerator denominator answer))
(newline)))))
1 ]=> (do ((i 0 (+ i 1))) ((>= i 1000)) (test-fmod (random 1.0) (random 1.0)))
; Value: #t
1 ]=> (test-fmod .59 .01)
Whoops: (.59 .01 1.0000000000000009e-2)
;Unspecified return value
Is this a question, a quiz, or a comment? The answer is that floating point is imprecise, and powers of 10 are one of the places it's especially imprecise, due to truncation of repeating decimals.
ReplyDeleteTo do this right, you have to specify an error bound, not use absolute equality.
It is an observation. We have a constraint that the remainder should be less than the divisor. For randomly chosen arguments, this is almost always true (you can pick thousands of random examples and not find a violation). It turns out that only a very small handful of argument pairs give a remainder larger than the divisor. This pair is the simplest one I've found with an exhaustive search.
ReplyDeleteI don't think that `floating point is imprecise' is a good argument here. I'm not measuring the accuracy of the answer. I don't expect it to be exactly right. But I do expect the remainder to be smaller than the divisor even if it is not exactly correct. In fact, since I am taking the floor, I am explicitly truncating the dividend to less accuracy than the implicit rounding. If rounding is a problem, then explicitly lopping off even more digits should produce even more problems, right? Yet the problematic argument pairs are quite rare.
The issue is more complicated than just `round off error'.
Curiously, trying (fmod ...) in Chicken (csi):
ReplyDelete(fmod 0.59 0.01) => 0.01
Didn't test exhaustively, but as in the above, the result was never greater than the denominator.
So maybe it is an implementation-specific issue. You didn't say which implementation you were running. If I had the time (and space on this laptop's smallish drives) it would be interesting to test several Schemes with (fmod). Which Scheme are you using?
Not sure what else I could add.
JRA.
Again: binary floating point doesn't represent 0.01 (or any power of ten) accurately. When you mix addition, multiplication, division, and integer operations the errors propagate.
ReplyDeleteHere's a visual binary analyzer:
http://babbage.cs.qc.edu/IEEE-754/
Enter your numbers, and you'll see that both 0.59 and 0.01 are repeating decimals in binary.
I don't understand why you think the errors won't propagate. I understand that you have that constraint, and you should either code so that the constraint is true or loosen it (I admit I was wrong to say that loosening the constraint was the only solution). Do you have some reason why your algorithm is immune to the usual floating point errors?
It's been too long since I took my applied linear algebra, so I'm not up to analyzing exactly why your error is negative rather than positive; but I don't suffer from any surprise at all in seeing that it is. It's the sort of thing I've seen many, many times before.
-Wm
In guile 2.0:
ReplyDelete(floor (/ .59 .01)) ; -> 59.0
But:
(define n .59)
(define d .01)
(floor (/ n d)) ; -> 58.0 (whoopsy!)
Why the problem happens is easy to understand - if division produces a quotient that is very very close to a multiple of 1.0, the floating-point error (lack of precision) might send it just below that multiple, and then floor imposes a huge error on the answer.
ReplyDeleteI reproduced the problem in my own Scheme, and this change "fixed" the problem by paying attention to what ROUND says:
(define (fmod numerator denominator)
(let* ((q (/ numerator denominator))
(md (- numerator (* denominator (round q)))))
(if (zero? md)
0.0
(- numerator
(* (floor q) denominator)))))
I was initially surprised that "zero?" was a sufficient check, but it caught every case that cropped up with the prior version. It no longer works once the quotient exceeds 1e16 or so (yes, double-precision floating point arithmetic).
Nevermind, my previous comment doesn't solve another problem of returning *negative* results even when n and d are both positive. This version is simpler, and behaves much better over a wide dynamic range.
ReplyDelete(define (fmod n d)
(let ((q (/ n d)))
(* d (- q (floor q)))))
I'm late to the party here, but I don't think the explanation has been sufficiently clearly pointed out as yet. So:
ReplyDeleteThe underlying cause is that the call to floor happens too late -- a rounding has already occurred in the division step, and that rounding is in the wrong direction. It so happens with these particular values and precisions that the representations have the following properties:
* The value a ~= 0.59 is very slightly greater than 59/100
* The value b ~= 0.01 is very slightly greater than 1/100
* The quotient a/b is very slightly less than 59
The result a/b is, of course, not exactly representable. Consequently the result of the division is rounded to the nearest representable value; this happens to be a rounding up to 59. Thereafter the application of floor applies to a value which is larger than the intended result it should apply to, and the above unexpected behaviour occurs.
(Note that this is the opposite of what Jim Rees suggested; the error is due to rounding up, not down. In fact, I'd be surprised if it were possible for a rounding down to cross a positive integer barrier.)
In order to fix this code you need control over the rounding behaviour of the division operator and set it explicitly to round to -infinity (or perhaps round to zero, but I'm avoiding thinking about negative values for the moment). In general, trying to apply specific rounding to only parts of a computation -- as has implicitly been done here -- is not going to work. One needs to ensure that the desired rounding properties are held at every step.
I am a bit surprised that you could find this behaviour with such small values; thanks for sharing that, it's a nice example that I hope to remember the next time I am explaining these kinds of issues to someone.