[Openmcl-devel] Double precision Floor bug?

[Brian Hayes]

This isn't really a Lisp issue; it's an IEEE-754 floating-point issue. And you needn't invoke anything as complicated as the floor function to see this behavior. Mere addition and subtraction show it clearly enough:

(setf *read-default-float-format* 'double-float)

? (- 0.59 0.58)
0.010000000000000009
? (+ 0.58 0.010000000000000009)
0.59
? (+ 0.58 0.01)
0.59
? (= (+ 0.58 0.010000000000000009) (+ 0.58 0.01))
t
? (= 0.010000000000000009 0.01)
nil

Two excellent papers --- one by Will Clinger, the other by Guy Steele and JonL White --- discuss the decimal-binary conversion problem at the root of this weirdness. The papers are in the 1990 ACM conference on Programming Language Design and Implementation (available through the ACM Digital Library, and probably elsewhere on the net if you go looking).

Another eye-opening paper (not Lisp-oriented): 

 Monniaux, David. 2008. The pitfalls of verifying floating-point computations.
 http://hal.archives-ouvertes.fr/docs/00/28/14/29/PDF/floating-point-article.pdf

The 2008 version of the IEEE-754 standard for floating point includes a decimal floating-point format. This would fix the problem in the specific case (- 0.59 0.58), but of course there are infinitely many real numbers that have no exact decimal representation, so it's not magic.

On the other hand, exact rationals are heavenly:

? (floor 59/100 1/100)
59
0

All those quaint notions like the associative law and the distributive law --- they actually work! (If only we could ignore the diagonals of squares and the circumferences of circles.)

Brian Hayes



On 3 Jun 2012, at 5:58 PM, Sudhir Shenoy wrote:

> Hi all,
> 
> There was a recent discussion on the Lispworks mailing list that I thought may be of interest for Clozure CL users. Essentially, someone was complaining about a possible bug in Lispworks where (floor 0.59d0 0.01d0) returned 0.010000000000000009D0. The reasonable expectation (even given the vagaries of inexact floating point arithmetic) that the remainder is always less than the divisor is violated.
> 
> Joe Marshall pointed out that this is probably due to the way floor is being computed (see below) and CLISP seems to do it the "right" way. CCL, which defines floor in terms of truncate gives
> 
> ? (truncate 0.59d0 0.01d0)
> 58
> 0.010000000000000009D0
> 
> which does seem strange ...
> 
> Cheers
> Sudhir
> 
> Begin forwarded message:
> 
>> From: Joe Marshall <jmarshall at alum.mit.edu>
>> Subject: Re: Double precision Floor bug?
>> Date: June 4, 2012 4:27:15 AM GMT+09:00
>> To: Guy Footring <Guy at footring.net>
>> Cc: lisp-hug at lispworks.com
>> Reply-To: Joe Marshall <jmarshall at alum.mit.edu>
>> 
>> This was bothering me this weekend so I did some more playing around.
>> Computing the quotient part of the answer (the first result) is straightforward:
>> you just divide and drop the fractional part.  It is the remainder that is
>> the problem.  There are several ways to compute the remainder, and
>> it seems that Lispworks does it this way:
>> 
>> (defun floor (n divisor)
>>   (let* ((quotient (/ n divisor))
>> 	 (integer-part (truncate quotient))
>>          (remainder (- n (* integer-part divisor))))
>>         (values integer-part remainder)))
>> 
>> But you could use the fractional part of the quotient instead like this:
>> 
>> (defun floor (n divisor)
>>   (let* ((quotient (/ n divisor))
>> 	 (integer-part (truncate quotient))
>> 	 (remainder (* (- quotient integer-part) divisor)))
>>         (values integer-part remainder)))
>> 
>> If the arithmetic were exact, these would produce the same answers,
>> but in floating point, the first produces
>> 58 0.010000000000000009d0
>> but the second produces
>> 58 0.009999999999999929d0
>> 
>> This second way of computing the remainder satisfies your expectations
>> about the constraints.
>>   
>> Incidentally, clisp appears to do it the second way.
>> [1]> (floor .59d0 .01d0)
>> 58 ;
>> 0.009999999999999929d0
>> 
>> 
>> -- 
>> ~jrm
> 
> _______________________________________________
> Openmcl-devel mailing list
> Openmcl-devel at clozure.com
> http://clozure.com/mailman/listinfo/openmcl-devel

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.clozure.com/pipermail/openmcl-devel/attachments/20120603/002a9d62/attachment.htm>