From: Peter K. <ps...@cs...> - 2011-09-02 07:24:25
|
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 |
From: Slobodan M. <slo...@gm...> - 2011-09-02 07:59:57
|
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 <ps...@cs...> > 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 > Sbc...@li... > 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 |
From: Peter K. <ps...@cs...> - 2011-09-02 08:30:50
|
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 |
From: Christophe R. <cs...@ca...> - 2011-09-02 08:17:55
|
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. |
From: Peter K. <ps...@cs...> - 2011-09-02 08:44:27
|
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 |
From: Christophe R. <cs...@ca...> - 2011-09-02 18:46:29
|
Hi, I'm glad I put the PS in and sent this to Raymond Toy, because he pointed out that... Christophe Rhodes <cs...@ca...> 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 |
From: Peter K. <ps...@cs...> - 2011-09-02 19:24:33
|
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 <cs...@ca...> 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 |
From: Christophe R. <cs...@ca...> - 2011-09-02 19:39:20
|
Peter Keller <ps...@cs...> 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 <cs...@ca...> 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 |
From: Peter K. <ps...@cs...> - 2011-10-07 03:42:27
|
On Fri, Sep 02, 2011 at 08:39:11PM +0100, Christophe Rhodes wrote: > Peter Keller <ps...@cs...> 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 |
From: Christophe R. <cs...@ca...> - 2011-10-07 08:12:25
|
Peter Keller <ps...@cs...> 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 |
From: James C. <clo...@jh...> - 2011-09-04 14:30:37
|
>>>>> "PK" == Peter Keller <ps...@cs...> 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 <cl...@jh...> OpenPGP: 1024D/ED7DAEA6 |
From: Christophe R. <cs...@ca...> - 2011-09-04 15:00:32
|
James Cloos <clo...@jh...> writes: >>>>>> "PK" == Peter Keller <ps...@cs...> 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 |
From: James C. <clo...@jh...> - 2011-09-06 17:32:11
|
>>>>> "CR" == Christophe Rhodes <cs...@ca...> 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 <cl...@jh...> OpenPGP: 1024D/ED7DAEA6 |