[Openmcl-devel] Fwd: Double precision Floor bug?

[Sudhir Shenoy]

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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.clozure.com/pipermail/openmcl-devel/attachments/20120604/904775c1/attachment.htm>