Summary

Numbers egg does not compute the same result as other schemes

Metadata

Attachments

Description

Dominic Pearson reported on #chicken the following problem:

When computing the sum i = 0 to n where n = 1000 of n^n chicken seems to produce an incorrect result.

I have attached the results and procedures used to this report.

Changes and comments

[2011-05-27 11:28:22 UTC] ckeen attached numbers-error.txt (description=#f)

[2011-05-27 11:31:07 UTC] ckeen wrote:

I have tested this on a 64 bit linux machine:

Linux foo 2.6.35-28-generic #50-Ubuntu SMP Fri Mar 18 18:42:20 UTC 2011 x86_64 GNU/Linux

[2011-05-29 23:52:56 UTC] sjamaan changed status from new to assigned

[2011-05-29 23:52:56 UTC] sjamaan changed owner from sjamaan to felix

[2011-05-29 23:52:56 UTC] sjamaan changed milestone from 4.7.0 to 4.8.0

[2011-05-29 23:52:56 UTC] sjamaan changed component from unknown to core libraries

[2011-05-29 23:52:56 UTC] sjamaan changed summary

[2011-05-29 23:52:56 UTC] sjamaan wrote:

So far it appears to be a bug in Chicken core's implementation of expt; it doesn't check errno and produces fixnum results when pow() overflows.

Updating summary and reassigning to Felix.

[2011-05-30 00:09:16 UTC] sjamaan wrote:

See also http://lists.nongnu.org/archive/html/chicken-users/2011-05/msg00067.html

[2011-05-30 08:36:22 UTC] felix wrote:

Would returning +inf on overflow, -inf on underflow and +nan on EDOM be sufficient?

[2011-05-30 09:19:11 UTC] sjamaan wrote:

I'm not sure. Probably...

Maybe John Cowan has something sensible to say about this :)

[2011-05-31 08:59:54 UTC] felix changed priority from major to critical

[2011-05-31 08:59:54 UTC] felix changed keywords from numbers eggs to numbers egg

[2011-05-31 08:59:54 UTC] felix wrote:

Replying to sjamaan: > I'm not sure. Probably... > > Maybe John Cowan has something sensible to say about this :)

He often does, but what he says isn't automatically sensible...

I will add a slight variation of the C `expt` routine in runtime.c to `numbers` and change the default procedure to return +/-inf and signal an error on EDOM (which would otherwise return a complex number). The latter change should be a Change Request, I think, since old code may stop working (even if it computed incorrect results).

[2011-05-31 09:00:45 UTC] felix changed version from 4.6.x to 4.7.x

[2011-05-31 14:21:17 UTC] felix wrote:

I committed a change to `numbers/trunk` that uses a fixnum-specific version of `C_expt`, but still get an incorrect result on a 64-bit machine (well, actually I forgot what the correct result was, but it is different from Racket's result).

[2011-05-31 14:21:31 UTC] felix changed owner from felix to sjamaan

[2011-05-31 14:22:26 UTC] felix changed summary

[2011-06-01 00:06:31 UTC] sjamaan wrote:

I think the real problem is that 64-bit floating point numbers only have 52 bits at their disposal; if the top-most bits of a 64 bit number _and_ the lowest bit of that number are used, pow() does not seem to set errno:

 gosh> (expt 999 6)
 994014980014994001
 gosh> (* 999 999 999 999 999 999)
 994014980014994001
 gosh> (* 999.0 999.0 999.0 999.0 999.0 999.0)
 9.94014980014994e17
 gosh> (inexact->exact (* 999.0 999.0 999.0 999.0 999.0 999.0))
 994014980014994048

As you can see in the last call, precision is lost. This can be reproduced in Chicken, and even simplified to:

 #;1> (inexact->exact (exact->inexact 994014980014994001))
 994014980014994048

There is no way to get around this I'm afraid; on 64-bits systems integers can represent larger numbers than doubles when they have their lowest and highest bits set :(

What's really annoying is that as far as I could find, nothing indicates this truncation even happened! The comparison between m1 and r doesn't help, since the truncation already happened before converting it to a C_word, so it compares as being equal.

I ripped out the special case for fixpoint args so it always takes the "slow" route in [[23864]|] and the calculation now gives the correct result.

If you agree with this and don't know a better solution either, please close the ticket and I'll tag a new release of "numbers". I'm unsure what should happen in core; I suppose it can just return the wrong result since there's nothing sane we can do.

God this sucks

[2011-06-01 09:03:40 UTC] felix changed status from assigned to closed

[2011-06-01 09:03:40 UTC] felix set resolution to worksforme

[2011-06-01 09:03:40 UTC] felix wrote:

Yeah, I agree completely. Very unfortunate. Thanks for figuring this out.

[2012-09-24 21:47:48 UTC] felix changed milestone from 4.8.0 to 4.9.0

[2012-09-24 21:47:48 UTC] felix wrote:

Milestone 4.8.0 deleted