www.digitalmars.com         C & C++   DMDScript  

digitalmars.D.bugs - Possible std.math.exp() bug

reply anthony.difranco yale.edu writes:
Ideas?  Help much appreciated; this is a huge obstacle right now.  Have
reproduced this on two different x86 linux systems.  Nothing special about 3.0
in reproducing this.  Have verified that the corresponding C library exp()
works.

real tgamma(real x);


$ cat test.d
import std.math;

int main(char[][] args)
{
printf("f(3.0) %f\n", sqrt(3.0));
printf("g(3.0) %f\n", exp(3.0));
printf("h(3.0) %f\n", tgamma(3.0));
return 1;
}

$ dmd test.d
gcc test.o -o test -lphobos -lpthread -lm

$ ./test
f(3.0) 1.732051
g(3.0) -0.000000
h(3.0) -0.000000

$ dmd
Digital Mars D Compiler v0.150
Copyright (c) 1999-2006 by Digital Mars written by Walter Bright
Documentation: www.digitalmars.com/d/index.html
Apr 03 2006
parent reply Don Clugston <dac nospam.com.au> writes:
anthony.difranco yale.edu wrote:

 Ideas?  Help much appreciated; this is a huge obstacle right now.  Have
 reproduced this on two different x86 linux systems.  Nothing special about 3.0
 in reproducing this.  Have verified that the corresponding C library exp()
 works.


The problem is that you're using %f with printf, but those functions 
return reals.
Try "%Lf" instead, or use writef.

However, there *is* an accuracy issue in tgamma, which occurred when it 
was translated from D back to C (!) for the DMC compiler. You can find 
the original, accurate tgamma() in the mathextra project at www.dsource.org.


 real tgamma(real x);

 
 $ cat test.d
 import std.math;
 
 int main(char[][] args)
 {
 printf("f(3.0) %f\n", sqrt(3.0));
 printf("g(3.0) %f\n", exp(3.0));
 printf("h(3.0) %f\n", tgamma(3.0));
 return 1;
 }
 
 $ dmd test.d
 gcc test.o -o test -lphobos -lpthread -lm
 
 $ ./test
 f(3.0) 1.732051
 g(3.0) -0.000000
 h(3.0) -0.000000
 
 $ dmd
 Digital Mars D Compiler v0.150
 Copyright (c) 1999-2006 by Digital Mars written by Walter Bright
 Documentation: www.digitalmars.com/d/index.html
Apr 03 2006
next sibling parent reply Walter Bright <newshound digitalmars.com> writes:
Don Clugston wrote:

 However, there *is* an accuracy issue in tgamma, which occurred when it 
 was translated from D back to C (!) for the DMC compiler.


I thought that was fixed? Can you redownload snn.lib?
Apr 04 2006
parent reply Don Clugston <dac nospam.com.au> writes:
Walter Bright wrote:

 Don Clugston wrote:

 However, there *is* an accuracy issue in tgamma, which occurred when 
 it was translated from D back to C (!) for the DMC compiler.


 
 I thought that was fixed? Can you redownload snn.lib?


When it was first included, it had a major bug in lgamma and a minor 
error in tgamma. The bug in lgamma() was fixed in the subsequent 
release. Last time I checked, tgamma was faulty (although it's subtle -- 
about the fifth digit is wrong).

I just tried to run my test program, but mathextra now gives an internal 
error (ztc.type.c 308) -- new regression in DMD 0.151, I'll track it down.
Apr 04 2006
parent Walter Bright <newshound digitalmars.com> writes:
Don Clugston wrote:

 Walter Bright wrote:

 Don Clugston wrote:

 However, there *is* an accuracy issue in tgamma, which occurred when 
 it was translated from D back to C (!) for the DMC compiler.


 I thought that was fixed? Can you redownload snn.lib?


 
 When it was first included, it had a major bug in lgamma and a minor 
 error in tgamma. The bug in lgamma() was fixed in the subsequent 
 release. Last time I checked, tgamma was faulty (although it's subtle -- 
 about the fifth digit is wrong).


I found the problem - forgot the 'L' suffix on the table entries.
Apr 07 2006
prev sibling parent reply anthony.difranco yale.edu writes:
Thanks.  I was actually using your tgamma from dsource.  I was also looking for
documentation for format 
strings for printf but found none in std.c.printf or std.string sections - where
is that covered?

Also, is there a better way to do the ratio of two gammas than subtracting
lgammas and exp at the end?


The problem is that you're using %f with printf, but those functions 
return reals.
Try "%Lf" instead, or use writef.

However, there *is* an accuracy issue in tgamma, which occurred when it 
was translated from D back to C (!) for the DMC compiler. You can find 
the original, accurate tgamma() in the mathextra project at www.dsource.org.
Apr 04 2006
parent reply Don Clugston <dac nospam.com.au> writes:
anthony.difranco yale.edu wrote:

 Thanks.  I was actually using your tgamma from dsource. 


Cool! I hope to make an official release of the mathextra libraries 
sometime soon. Some of the functions in std.complex have known issues, 
but the statistics ones are very well tested.

Good to have another mathematics programmer around!

  I was also looking for

 documentation for format 
 strings for printf but found none in std.c.printf or std.string sections -
where
 is that covered?


It's in the docs for the DMC standard library. There are a few things 
like this that aren't included in the downloaded docs, and which should 
at least be linked to.

http://www.digitalmars.com/rtl/stdio.html#fprintf


 
 Also, is there a better way to do the ratio of two gammas than subtracting
 lgammas and exp at the end?


Do you mean, more accurate way? If the numbers are not too large, I 
don't think it matters much, even a simple division of the tgammas is 
probably OK, but lgamma() will save you from overflow problems.
If the arguments are large, the Stirling approximation is used, so there 
may be potential for greater accuracy. I'd be amazed if it actually 
mattered, though.
Apr 04 2006
parent reply anthony.difranco yale.edu writes:
Good to have another mathematics programmer around!


I'm trying.  D has been very kind to me so far, with both reasonable semantics
and reasonable efficiency, and now a reasonable community.


 Also, is there a better way to do the ratio of two gammas than subtracting
 lgammas and exp at the end?


Do you mean, more accurate way? If the numbers are not too large, I 
don't think it matters much, even a simple division of the tgammas is 
probably OK, but lgamma() will save you from overflow problems.
If the arguments are large, the Stirling approximation is used, so there 
may be potential for greater accuracy. I'd be amazed if it actually 
mattered, though.


Actually, I should have said an as-overflow-resistant-as-possible way to do
something like permutations (basically ratio of gammas).  Or maybe some useful
scaling identity that doesn't wreck precision, though squeezing out the last few
bits is not a concern.  My application (statistics related) is making even
lgamma overflow as a matter of course.
Apr 06 2006
parent reply Don Clugston <dac nospam.com.au> writes:
anthony.difranco yale.edu wrote:

 Good to have another mathematics programmer around!


 I'm trying.  D has been very kind to me so far, with both reasonable semantics
 and reasonable efficiency, and now a reasonable community.
 

 Also, is there a better way to do the ratio of two gammas than subtracting
 lgammas and exp at the end?


 Do you mean, more accurate way? If the numbers are not too large, I 
 don't think it matters much, even a simple division of the tgammas is 
 probably OK, but lgamma() will save you from overflow problems.
 If the arguments are large, the Stirling approximation is used, so there 
 may be potential for greater accuracy. I'd be amazed if it actually 
 mattered, though.


 
 Actually, I should have said an as-overflow-resistant-as-possible way to do
 something like permutations (basically ratio of gammas).  Or maybe some useful
 scaling identity that doesn't wreck precision, though squeezing out the last
few
 bits is not a concern.  My application (statistics related) is making even
 lgamma overflow as a matter of course.



For x > 1e10, it looks as though this is all lgamma() is doing:

const real LOGSQRT2PI  =  0.91893853320467274178L;

return ( x - 0.5L ) * log(x) - x + LOGSQRT2PI;

Disturbing. This doesn't look very accurate, although Stirling's 
approximation is probably very good by then.

So for large a, b
lgamma(a) - lgamma(b)
= a*(log(a)-1) - b*(log(b)-1) - 0.5 * (log(a) - log(b));
So overflows can be avoided quite easily ... but is any accuracy left?
Apr 07 2006
parent anthony.difranco yale.edu- writes:
Thanks for the view from within, but it ended up being something else causing
the overflow, and regular lgamma is fine in the correct context.  I'll file that
away for forays into the combinatorics of ridiculous quantities.


For x > 1e10, it looks as though this is all lgamma() is doing:

const real LOGSQRT2PI  =  0.91893853320467274178L;

return ( x - 0.5L ) * log(x) - x + LOGSQRT2PI;

Disturbing. This doesn't look very accurate, although Stirling's 
approximation is probably very good by then.

So for large a, b
lgamma(a) - lgamma(b)
= a*(log(a)-1) - b*(log(b)-1) - 0.5 * (log(a) - log(b));
So overflows can be avoided quite easily ... but is any accuracy left?
Apr 12 2006