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);
The Gamma function, Γ(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

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);
> The Gamma function, Γ(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
> 
> 

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?

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.

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

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.

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

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.

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?

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?