Numbers egg does not compute the same result as other schemes

[Christian Kellermann]

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.

[Christian Kellermann]

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

[sjamaan]

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.

[sjamaan]

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

[felix winkelmann]

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

[sjamaan]

I'm not sure. Probably...

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

[felix winkelmann]

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

[felix winkelmann]

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

[sjamaan]

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

[felix winkelmann]

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

[felix winkelmann]

Milestone 4.8.0 deleted