x86 double-float issue

[Raymond Toy]

Consider this sample code

(defun rbug (z tt betain beta)
 (declare (double-float z tt betain beta)
      (optimize (speed 3) (safety 0)))
 (= (* (* (* z tt) betain) beta) z))

Compile it and run:

(rbug 5.562684646268004d-309 (1+ double-float-epsilon) .5d0 2d0)

This is t on sparc, nil on x86.

This is caused by x86 issues with double-float (53-bit) precision and the extra range of the exponent in long-double (64-bit) format. Denormals in this format are not handled the same as denormals on sparc or ppc.

There is a fix for this. When computing x*y, we scale x by an appropriate value, multiply by y, and scale back. This will produce the correctly rounded denormal. The only issue would be the exponent range. That is fixed by storing the number to memory and reloaded. (This solution taken from a proposed solution for Java numerics).

This could be easily implemented, but potentially slows down double-float arithmetic by a factor of 2-4 times. I don't think we want to do that.

[Raymond Toy]

This should be fixed in the 2008-12 snapshot when sse2 is used. This issue will remain for x87; we have no intention of fixing this.