Re: New proposal for math.e
Matt Lewis wrote:
>
> 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:
> }}}
<eucode>
> 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
> </eucode>
{{{
>
> 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:
> }}}
<eucode>
> global function logn( object x, atom base )
> </eucode>
{{{
>
> Matt
Your solution for the constants is certainly the most rigorous one. Let humans
read material for humans and computers read material for computers, without
interference.
As for deg_to_rad's name: too many keystrokes as usual, but I'll pretend I don't
care, since I don't think I ever had to write it before this post.
logn(x,base)=log(x)/log(base)? I have nothing against it, but the doc should
emphasize that there may be tiny discrepancies, for instance between LN2 and
logn(2,E), because of the just discussed rounding errors.
Since no one replied to my first proposal, I'll repeat it: add a gamma()
function and an erf() function, since they re very commonly used in analysis and
statistics. Note that an incomplete gamma function could be an even better idea
than erf(), as erf() is a special, very frequent occurrence of
incomplete_gamma(x,z) after a trivial transform. The general incomplete_gamma()
arises in physics rather. The value of the Euler constant (=-gamma'(1)) would be
a logical addition then too.
1/sqrt(2*PI) is a very common normalisation constant as well, both in statistics
and in Fourier series calculations, so it may not hurt to have it on board as
well.
The most difficult thing here could be to port some Fortran code to Eu to do the
above, assuming the C compiler don't have this as standard. There is ample supply
of such free code.
CChris
|
Not Categorized, Please Help
|
|
