OpenEuphoria

Re: New proposal for math.e

• Posted by Matt Lewis <matthewwalkerlewis at gma?l?com> Aug 03, 2007
• 568 views

Juergen Luethje wrote:
> 
> Me wrote:
> 
> <snip>
> 
> > Using the code:
> > }}}
<eucode>
> > function log2 (atom x)
> >    return log(x)/log(2)
> > end function
> >  
> > for i = 0 to 32 do
> >    printf(1, "%d: %0.16g\n", i & log2( power( 2, i ) ))
> > end for
> > </eucode>
{{{

> > with the EXW.EXE interpreter 3.0.2, I can't manage to get imprecise results
> > anymore, neither on my old system nor on my new system.  *LOL*
> 
> <snip>
> 
> I think know I've found the explanation.
> Matt, when I had tested my log2() function, I probably had used a different
> way than you to look for imprecise results. Please test this:

<snip>

> These are the two cases that I treated specially in my proposal for
> "math.e".
> 
> So it seems that for i = 29 and for i = 31, the results of log2() are
> not equal according to the != operator, but they are (or only look)
> equal according to printf() ...

OK, I see the difference now.  Using the straight log/log method, there
is still a difference for 29 and 31.

OTOH, using the constant, the difference for 29 and 31 was in the next to 
least significant bit, making the differences twice as big.  And there 
were many more differences.  Definitely looks like rounding errors.

And by looking at what I got by parsing the constant, I can see that 
by letting the interpreter parse your constant, it loses accuracy.
Below, I've shown the bits of the double gotten in these different
ways:

1: log(2)
2: 0.6931471805599453
3: 0.6931471805599453e0 (using my parser)

1: 1111011110011100010111110111111101000010011101000110 01111111110 0 (-1)
2: 0111011110011100010111110111111101000010011101000110 01111111110 0 (-1)
3: 1111011110011100010111110111111101000010011101000110 01111111110 0 (-1)

This certainly calls into question the other constants, since they might
also be one bit off (probably ok for most applications, but still...)
For log(2), we could just do the calculation on assignment, but that 
doesn't help with the others.  The easiest solution is to use the
float64s, and call machine_func(46,...) on them (to avoid including
machine.e) and putting the human readable value in a comment:

LN2        = machine_func(47, {239,57,250,254,66,46,230,63}),   --
    0.6931471805599453e0
LN10       = machine_func(47, {21,85,181,187,177,107,2,64}),    --
    2.3025850929940456e0
E          = machine_func(47, {105,87,20,139,10,191,5,64}),     --
    2.7182818284590452e0
SQRT2      = machine_func(47, {204,59,127,102,158,160,246,63}), --
    1.4142135623730950e0
HALF_SQRT2 = machine_func(47, {204,59,127,102,158,160,230,63}), --
    0.7071067811865475e0
PI         = machine_func(47, {24,45,68,84,251,33,9,64}),       --
    3.1415926535897932e0
HALF_PI    = machine_func(47, {24,45,68,84,251,33,249,63}),     --
    1.5707963267948966e0
QUARTER_PI = machine_func(47, {24,45,68,84,251,33,233,63}),     --
    0.7853981633974483e0
TWO_PI     = machine_func(47, {24,45,68,84,251,33,25,64})       --
    6.2831853071795864e0

Using the constant for LN2 (as long as it's the *right* one) is probably
faster (I think that's a reasonable assumption), we should probably go
with your original function, but without the special cases.

For a new proposal, I know that we've covered the most common cases, but
how about a logn (or logx?) function for the general case:

global function logn( object x, atom base )

Matt