Re: [Sbcl-help] rounding issue

[Christophe R. <cs...@ca...>]

Peter Keller <ps...@cs...> writes:

> I'm writing some math software using SBCL 1.0.47 (and 1.0.46 gives the
> same response) and I'm curious why I get this output:
>
> * (round (* 10.004939019999991d0 1e8))

(It's pretty clearly a bug, though finding out where is a bit more
miserable; I'm always nervous about things on x86 because of the
internal 80-bit precision mode.) This happens to me on x86 but not
amd64; the answer from sb-kernel:%unary-round is different.

amd64:
* (round (* 10.004939019999991d0 1d8))
  0: (SB-KERNEL:%UNARY-ROUND 1.0004939019999992d9)
  0: SB-KERNEL:%UNARY-ROUND returned 1000493902

x86:
* (round (* 10.004939019999991d0 1d8))
  0: (SB-KERNEL:%UNARY-ROUND 1.0004939019999992d9)
  0: SB-KERNEL:%UNARY-ROUND returned 1000493901

The answer I think is that the not-fixnum branch in %unary-round is
completely wrong.  The logic looks like a weird mix of truncate and
round-to-even, getting I think neither of them right.  I think if we
replace the ROUNDED binding with
  (if (minusp exp)
      (let ((fractional-bits (logand bits (lognot (ash -1 (- exp)))))
            (0.5bits (ash 1 (- -1 exp))))
        (cond
          ((> fractional-bits 0.5bits) (1+ shifted))
          ((< fractional-bits 0.5bits) (1+ shifted))
          (t (if (oddp shifted) (1+ shifted) shifted)))))
then it'll start working.  (It passes my minimal tests, but maybe you
could check that the rest of your application, not just this one round
statement, works?)

Thank you for the report!

Best,

Christophe

PS: if there are any cmucl maintainers reading this, the version in
cmucl looks differently wrong but still wrong -- it will only ever round
up if the integer part of the float is odd.