Steel Bank Common Lisp

Hi,

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)

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?

2011/9/2 Peter Keller <psilord@...>

> 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?
>
> 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
> _______________________________________________
> Sbcl-help mailing list
> Sbcl-help@...
> https://lists.sourceforge.net/lists/listinfo/sbcl-help
>



-- 
Poći ću s vama jer volim šalu, hoću da vidim ježa budalu.
Put u Japan - http://ofcan.wordpress.com

Peter Keller <psilord@...> 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.

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

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

Hi,

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

Christophe Rhodes <csr21@...> 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

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@...> 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

Peter Keller <psilord@...> 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@...> 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

>>>>> "PK" == Peter Keller <psilord@...> 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@...>         OpenPGP: 1024D/ED7DAEA6

James Cloos <cloos+sf-sbcl-help@...> writes:

>>>>>> "PK" == Peter Keller <psilord@...> 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

>>>>> "CR" == Christophe Rhodes <csr21@...> writes:

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

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

Yes.  I should've read the whole thread before replying....

Apologies for the noise.

-JimC
-- 
James Cloos <cloos@...>         OpenPGP: 1024D/ED7DAEA6

On Fri, Sep 02, 2011 at 08:39:11PM +0100, Christophe Rhodes wrote:
> Peter Keller <psilord@...> writes:
> 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.

Hello,

So it turns out your fix isn't quite right:

Linux > sbcl
This is SBCL 1.0.51.41-6113d10, 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.

;; This works
* (round (* 10.004939019999991d0 1e8))

1000493902
-8.344650268554688d-7

;; this works too.
* (round 100000000000000.6d0)

100000000000001
-0.40625d0

;; uh oh!
* (round 100000000000000.6)

debugger invoked on a TYPE-ERROR in thread #<THREAD "initial thread" RUNNING
                                              {AAA58B1}>:
  The value NIL is not of type INTEGER.

Type HELP for debugger help, or (SB-EXT:QUIT) to exit from SBCL.

restarts (invokable by number or by possibly-abbreviated name):
  0: [ABORT] Exit debugger, returning to top level.

(SB-KERNEL:%UNARY-ROUND 1.e14)
0] 


When the 100000000000000.6 value is read as a single-float, round messes up
and thinks it isn't an integer. Same thing happens with (round 1.e14) and
(round 1e14).

Thank you.

-pete

Peter Keller <psilord@...> writes:

> * (round 100000000000000.6)
>
> debugger invoked on a TYPE-ERROR in thread #<THREAD "initial thread" RUNNING
>                                               {AAA58B1}>:
>   The value NIL is not of type INTEGER.

Um, how did that happen?  There's a missing alternative branch in the IF
in the binding of ROUNDED in %unary-round?  People shouldn't let me
drive.

It should read
  (if (minusp exp)
      (... complicated fractional-bits stuff ...)
      shifted)
as any fule kno.

Sorry.  Juho, can I slip this in under your freeze radar?

Cheers,

Christophe