Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

Earlier, I asked a question on MathOverflow regarding how one might analytically approximate a function of the form: $f(n) = \prod_{i=1}^{n-1} (1-ai)$ for $a \ge 0$, $(ai) < 1$, and $n > 10^5$ or $10^6$. Robert Israel answered the question in a very nice way, observing that:

$f(n) = \dfrac{a^{n-1} \Gamma(1/a)}{\Gamma(1-n+1/a)}$

...and that one could use Stirling's series for an asymptotic approximation of $\ln(\Gamma(1-n+ 1/a))$.


My question here is: how might one actually compute values for this function (using, for example, Mathematica) for very large $n$ ($n > 10^9$ or so) and very small $a$ ($a < 10^{-30}$), without overflows/underflows, provided that the known output for $f(n)$ falls in the range of, say, $10^{-12} \le f(n) \le 1$? Are there any simple strategies for achieving this?

share|improve this question
add comment

1 Answer

up vote 1 down vote accepted

Since this is a common problem many software implemented a lngamma function.
In the case of Mathematica it is named LogGamma[].

If this fails you may use an expansion of $\ln(\Gamma(x))$. You may find these at Wikipedia.

Of course once the LogGamma obtained you shouldn't compute directly the exponential ! Instead evaluate the $\log$ of the complete expression and take the exponential of the (simplified) result.

share|improve this answer
    
Thanks for the suggestion... however, even the LogGamma function fails to avoid overflows/underflows. –  Roger S. Aug 16 '12 at 15:07
    
@RogerS.: yes of course but in a much enlarged space ! Which is the value in trouble ? –  Raymond Manzoni Aug 16 '12 at 15:10
    
@Roger, do you have an example? –  J. M. Aug 16 '12 at 15:12
2  
@Roger: You're using LogGamma[] wrong. Recall the identity $\log(p/q)=\log\,p-\log\,q$... –  J. M. Aug 16 '12 at 15:51
1  
@Roger, With[{n = 10^6, a = 10^-30}, Block[{$MaxExtraPrecision = 60}, N[(n - 1) Log[a] + LogGamma[1/a] - LogGamma[1 - n + 1/a], 25]]] –  J. M. Aug 16 '12 at 16:23
show 8 more comments

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.