From: Peter Keller <psilord@cs...>  20110902 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.344650268554688d7 On ecl 11.1.1 I get this: > (round (* 10.004939019999991d0 1e8)) 1000493902 8.344650268554688d7 Needless to say, I spent quite a few days tracking this down. Is there a known bug inside of SBCL concerning doublefloat 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 Milnović <slobodan.milnovic@gm...>  20110902 07:59:57
Attachments:
Message as HTML

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 x8664 debian) CLUSER> (round (* 10.004939019999991d0 1e8)) 1000493902 8.344650268554688d7 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.344650268554688d7 > > On ecl 11.1.1 I get this: > > > (round (* 10.004939019999991d0 1e8)) > > 1000493902 > 8.344650268554688d7 > > Needless to say, I spent quite a few days tracking this down. Is there > a known bug inside of SBCL concerning doublefloat 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 worldclass log management solution at an even better > pricefree! And you'll get a free "Love Thy Logs" tshirt when you > download Logger. Secure your free ArcSight Logger TODAY! > http://p.sf.net/sfu/arcsisghtdev2dev > _______________________________________________ > Sbclhelp mailing list > Sbclhelp@... > https://lists.sourceforge.net/lists/listinfo/sbclhelp >  Poći ću s vama jer volim šalu, hoću da vidim ježa budalu. Put u Japan  http://ofcan.wordpress.com 
From: Peter Keller <psilord@cs...>  20110902 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 x8664 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 32bit box). > CLUSER> (round (* 10.004939019999991d0 1e8)) > 1000493902 > 8.344650268554688d7 > > 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 BSDstyle 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 nouserinit, 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 Rhodes <csr21@ca...>  20110902 08:17:55

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 80bit precision mode.) This happens to me on x86 but not amd64; the answer from sbkernel:%unaryround is different. amd64: * (round (* 10.004939019999991d0 1d8)) 0: (SBKERNEL:%UNARYROUND 1.0004939019999992d9) 0: SBKERNEL:%UNARYROUND returned 1000493902 x86: * (round (* 10.004939019999991d0 1d8)) 0: (SBKERNEL:%UNARYROUND 1.0004939019999992d9) 0: SBKERNEL:%UNARYROUND returned 1000493901 The answer I think is that the notfixnum branch in %unaryround is completely wrong. The logic looks like a weird mix of truncate and roundtoeven, getting I think neither of them right. I think if we replace the ROUNDED binding with (if (minusp exp) (let ((fractionalbits (logand bits (lognot (ash 1 ( exp))))) (0.5bits (ash 1 ( 1 exp)))) (cond ((> fractionalbits 0.5bits) (1+ shifted)) ((< fractionalbits 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 Keller <psilord@cs...>  20110902 08:44:27

On Fri, Sep 02, 2011 at 09:17:45AM +0100, Christophe Rhodes wrote: > The answer I think is that the notfixnum branch in %unaryround is > completely wrong. The logic looks like a weird mix of truncate and > roundtoeven, getting I think neither of them right. I think if we > replace the ROUNDED binding with > (if (minusp exp) > (let ((fractionalbits (logand bits (lognot (ash 1 ( exp))))) > (0.5bits (ash 1 ( 1 exp)))) > (cond > ((> fractionalbits 0.5bits) (1+ shifted)) > ((< fractionalbits 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 Rhodes <csr21@ca...>  20110902 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 <csr21@...> writes: > The answer I think is that the notfixnum branch in %unaryround is > completely wrong. The logic looks like a weird mix of truncate and > roundtoeven, getting I think neither of them right. I think if we > replace the ROUNDED binding with > (if (minusp exp) > (let ((fractionalbits (logand bits (lognot (ash 1 ( exp))))) > (0.5bits (ash 1 ( 1 exp)))) > (cond > ((> fractionalbits 0.5bits) (1+ shifted)) > ((< fractionalbits 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 Keller <psilord@cs...>  20110902 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 <csr21@...> writes: > > > The answer I think is that the notfixnum branch in %unaryround is > > completely wrong. The logic looks like a weird mix of truncate and > > roundtoeven, getting I think neither of them right. I think if we > > replace the ROUNDED binding with > > (if (minusp exp) > > (let ((fractionalbits (logand bits (lognot (ash 1 ( exp))))) > > (0.5bits (ash 1 ( 1 exp)))) > > (cond > > ((> fractionalbits 0.5bits) (1+ shifted)) > > ((< fractionalbits 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 Rhodes <csr21@ca...>  20110902 19:39:20

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 notfixnum branch in %unaryround is >> > completely wrong. The logic looks like a weird mix of truncate and >> > roundtoeven, getting I think neither of them right. I think if we >> > replace the ROUNDED binding with >> > (if (minusp exp) >> > (let ((fractionalbits (logand bits (lognot (ash 1 ( exp))))) >> > (0.5bits (ash 1 ( 1 exp)))) >> > (cond >> > ((> fractionalbits 0.5bits) (1+ shifted)) >> > ((< fractionalbits 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 roundtoeven logic, not at all confusing. Cheers, Christophe 
From: Peter Keller <psilord@cs...>  20111007 03:42:27

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 roundtoeven logic, not at all confusing. Hello, So it turns out your fix isn't quite right: Linux > sbcl This is SBCL 1.0.51.416113d10, 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 BSDstyle licenses. See the CREDITS and COPYING files in the distribution for more information. ;; This works * (round (* 10.004939019999991d0 1e8)) 1000493902 8.344650268554688d7 ;; this works too. * (round 100000000000000.6d0) 100000000000001 0.40625d0 ;; uh oh! * (round 100000000000000.6) debugger invoked on a TYPEERROR in thread #<THREAD "initial thread" RUNNING {AAA58B1}>: The value NIL is not of type INTEGER. Type HELP for debugger help, or (SBEXT:QUIT) to exit from SBCL. restarts (invokable by number or by possiblyabbreviated name): 0: [ABORT] Exit debugger, returning to top level. (SBKERNEL:%UNARYROUND 1.e14) 0] When the 100000000000000.6 value is read as a singlefloat, 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 Rhodes <csr21@ca...>  20111007 08:12:25

Peter Keller <psilord@...> writes: > * (round 100000000000000.6) > > debugger invoked on a TYPEERROR 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 %unaryround? People shouldn't let me drive. It should read (if (minusp exp) (... complicated fractionalbits stuff ...) shifted) as any fule kno. Sorry. Juho, can I slip this in under your freeze radar? Cheers, Christophe 
From: James Cloos <cloos+<sfsbclhelp@jh...>  20110904 14:30:37

>>>>> "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.344650268554688d7 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 
From: Christophe Rhodes <csr21@ca...>  20110904 15:00:32

James Cloos <cloos+sfsbclhelp@...> 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.344650268554688d7 > > It looks to be an fpmath=387 vs fpmath=sse issue. It is not. It is a sizeofmostpositivefixnum issue, coupled with a bug. Christophe 
From: James Cloos <cloos+<sfsbclhelp@jh...>  20110906 17:32:11

>>>>> "CR" == Christophe Rhodes <csr21@...> writes: JC>> It looks to be an fpmath=387 vs fpmath=sse issue. CR> It is not. It is a sizeofmostpositivefixnum 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 