Possible std.math.exp() bug

Don Clugston dac at nospam.com.au
Fri Apr 7 06:20:05 PDT 2006

```anthony.difranco at 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?

```