rounding issue

Hello,

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))

1000493901
0.9999991655349731d0

I don't expect that output....

On CLISP 2.48 I get this:

[1]> (round (* 10.004939019999991d0 1d8))
1000493902 ;                                                                    
-8.344650268554688d-7

On ecl 11.1.1 I get this:

> (round (* 10.004939019999991d0 1e8))

1000493902
-8.344650268554688d-7

Needless to say, I spent quite a few days tracking this down. Is there
a known bug inside of SBCL concerning double-float representation that
I'm hitting? Or am I so close to the available manitissa bits that I'm
hitting undefined behavior?

Re: rounding issue

------------------------------------------------------------------------------
Special Offer -- Download ArcSight Logger for FREE!
Finally, a world-class log management solution at an even better 
price-free! And you'll get a free "Love Thy Logs" t-shirt when you
download Logger. Secure your free ArcSight Logger TODAY!
http://p.sf.net/sfu/arcsisghtdev2dev

_______________________________________________
Sbcl-help mailing list
Sbcl-help <at> lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/sbcl-help

Re: rounding issue

Peter Keller <psilord <at> cs.wisc.edu> 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.

------------------------------------------------------------------------------
Special Offer -- Download ArcSight Logger for FREE!
Finally, a world-class log management solution at an even better 
price-free! And you'll get a free "Love Thy Logs" t-shirt when you
download Logger. Secure your free ArcSight Logger TODAY!
http://p.sf.net/sfu/arcsisghtdev2dev

Re: rounding issue

On Fri, Sep 02, 2011 at 09:59:31AM +0200, Slobodan Milnovi?? wrote:
> I'm just an sbcl user like you, and I have tried the same with the latest
> stable 1.0.51 (compiled from source on x86-64 debian)

I forgot to mention in the first email (and I apologize) that
I'm on an x86 ubuntu 10.04.3 LTS box (and those other tests were on a
32-bit box). 

> CL-USER> (round (* 10.004939019999991d0 1e8))
> 1000493902
> -8.344650268554688d-7
> 
> So, my guess is that whatever issue has been with the 1.0.46/47, sbcl
> developers seem to have solved it with 1.0.51. Perhaps you could test your
> code with that one and see if it works?

Sure, I'll give 1.0.51 a try.

Hrm. Nope. It still doesn't work.

Linux black > sbcl
This is SBCL 1.0.51, an implementation of ANSI Common Lisp.
More information about SBCL is available at <http://www.sbcl.org/>.

SBCL is free software, provided as is, with absolutely no warranty.
It is mostly in the public domain; some portions are provided under
BSD-style licenses.  See the CREDITS and COPYING files in the
distribution for more information.
* (round (* 10.004939019999991d0 1e8))

1000493901
0.9999991655349731d0

If I start it with --no-userinit, the round call still doesn't work.

Since it worked for you, this has the signature of a 32/64 bit problem.

My specific processor is this:

Linux black > cat /proc/cpuinfo 
processor	: 0
vendor_id	: GenuineIntel
cpu family	: 6
model		: 13
model name	: Intel(R) Celeron(R) M processor         1.60GHz
stepping	: 8
cpu MHz		: 1596.238
cache size	: 1024 KB
fdiv_bug	: no
hlt_bug		: no
f00f_bug	: no
coma_bug	: no
fpu		: yes
fpu_exception	: yes
cpuid level	: 2
wp		: yes
flags		: fpu vme de pse tsc msr pae mce cx8 apic mtrr pge mca cmov clflush dts acpi mmx fxsr sse sse2 ss tm pbe nx
up bts
bogomips	: 3192.47
clflush size	: 64
cache_alignment	: 64
address sizes	: 32 bits physical, 32 bits virtual
power management:

Thank you for your help!

-pete

------------------------------------------------------------------------------
Special Offer -- Download ArcSight Logger for FREE!
Finally, a world-class log management solution at an even better 
price-free! And you'll get a free "Love Thy Logs" t-shirt when you
download Logger. Secure your free ArcSight Logger TODAY!
http://p.sf.net/sfu/arcsisghtdev2dev

Re: rounding issue

On Fri, Sep 02, 2011 at 09:17:45AM +0100, Christophe Rhodes wrote:
> 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?)

I'll try out the suggested patch as soon as I can. I need sleep badly and
tomorrow is a busy day. So it might be a day before I get a chance to try
it.

> Thank you for the report!

Thank you for the amazingly fast debugging of it! I very much appreciate it!

-pete

------------------------------------------------------------------------------
Special Offer -- Download ArcSight Logger for FREE!
Finally, a world-class log management solution at an even better 
price-free! And you'll get a free "Love Thy Logs" t-shirt when you
download Logger. Secure your free ArcSight Logger TODAY!
http://p.sf.net/sfu/arcsisghtdev2dev

Re: rounding issue

Hi,

I'm glad I put the PS in and sent this to Raymond Toy, because he
pointed out that...

Christophe Rhodes <csr21 <at> cantab.net> writes:

> 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))

... this line is wrong.  "Why increment if the fractional part is less
that 0.5?"  (Why indeed?)

>           (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?)

Cheers,

Christophe

------------------------------------------------------------------------------
Special Offer -- Download ArcSight Logger for FREE!
Finally, a world-class log management solution at an even better 
price-free! And you'll get a free "Love Thy Logs" t-shirt when you
download Logger. Secure your free ArcSight Logger TODAY!
http://p.sf.net/sfu/arcsisghtdev2dev

Re: rounding issue

On Fri, Sep 02, 2011 at 07:46:20PM +0100, Christophe Rhodes wrote:
> I'm glad I put the PS in and sent this to Raymond Toy, because he
> pointed out that...
> 
> Christophe Rhodes <csr21 <at> cantab.net> writes:
> 
> > 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))
> 
> ... this line is wrong.  "Why increment if the fractional part is less
> that 0.5?"  (Why indeed?)

So I should just remove the < condition code path from the COND clause?

> >           (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!

-pete

------------------------------------------------------------------------------
Special Offer -- Download ArcSight Logger for FREE!
Finally, a world-class log management solution at an even better 
price-free! And you'll get a free "Love Thy Logs" t-shirt when you
download Logger. Secure your free ArcSight Logger TODAY!
http://p.sf.net/sfu/arcsisghtdev2dev

Re: rounding issue

Peter Keller <psilord <at> cs.wisc.edu> writes:

> On Fri, Sep 02, 2011 at 07:46:20PM +0100, Christophe Rhodes wrote:
>> I'm glad I put the PS in and sent this to Raymond Toy, because he
>> pointed out that...
>> 
>> Christophe Rhodes <csr21 <at> cantab.net> writes:
>> 
>> > 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))
>> 
>> ... this line is wrong.  "Why increment if the fractional part is less
>> that 0.5?"  (Why indeed?)
>
> So I should just remove the < condition code path from the COND clause?

no, sorry, that line should still be there but return "shifted" instead
of "(1+ shifted)".  So: if the fractional bits are bigger than 0.5; then
use (1+ shifted); if they're smaller than 0.5, use shifted; if they are
exactly 0.5 then if the shifted value is odd, use (1+ shifted) otherwise
use shifted.  Proper round-to-even logic, not at all confusing.

Cheers,

Christophe

------------------------------------------------------------------------------
Special Offer -- Download ArcSight Logger for FREE!
Finally, a world-class log management solution at an even better 
price-free! And you'll get a free "Love Thy Logs" t-shirt when you
download Logger. Secure your free ArcSight Logger TODAY!
http://p.sf.net/sfu/arcsisghtdev2dev

Re: rounding issue

>>>>> "PK" == Peter Keller <psilord <at> cs.wisc.edu> writes:

PK> * (round (* 10.004939019999991d0 1e8))

PK> 1000493901
PK> 0.9999991655349731d0

I get that on my ia32 box, but my amd64 box gives:

1000493902
-8.344650268554688d-7

It looks to be an fpmath=387 vs fpmath=sse issue.

The ia32 ABI passes all floats, doubles and long doubles in the x387
registers, which are 80 bit registers.

The amd64 ABI passes floats and doubles via the xmm registers, so floats
stay 32 bits and doubles stay 64 bits.

In short, this is a rounding issue which is dependent on how floating
point values are stored on each CPU.  

-JimC
--

-- 
James Cloos <cloos <at> jhcloos.com>         OpenPGP: 1024D/ED7DAEA6

------------------------------------------------------------------------------
Special Offer -- Download ArcSight Logger for FREE!
Finally, a world-class log management solution at an even better 
price-free! And you'll get a free "Love Thy Logs" t-shirt when you
download Logger. Secure your free ArcSight Logger TODAY!
http://p.sf.net/sfu/arcsisghtdev2dev

Re: rounding issue

James Cloos <cloos+sf-sbcl-help <at> jhcloos.com> writes:

>>>>>> "PK" == Peter Keller <psilord <at> cs.wisc.edu> writes:
>
> PK> * (round (* 10.004939019999991d0 1e8))
>
> PK> 1000493901
> PK> 0.9999991655349731d0
>
> I get that on my ia32 box, but my amd64 box gives:
>
> 1000493902
> -8.344650268554688d-7
>
> It looks to be an fpmath=387 vs fpmath=sse issue.

It is not.  It is a size-of-most-positive-fixnum issue, coupled with a
bug.

Christophe

------------------------------------------------------------------------------
Special Offer -- Download ArcSight Logger for FREE!
Finally, a world-class log management solution at an even better 
price-free! And you'll get a free "Love Thy Logs" t-shirt when you
download Logger. Secure your free ArcSight Logger TODAY!
http://p.sf.net/sfu/arcsisghtdev2dev