3
$\begingroup$

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?

1 Answers 1

1

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.

  • 0
    Thanks for the suggestion... however, even the LogGamma function fails to avoid overflows/underflows.2012-08-16
  • 0
    @RogerS.: yes of course but in a much enlarged space ! Which is the value in trouble ?2012-08-16
  • 0
    @Roger, do you have an example?2012-08-16
  • 0
    @RogerS.: btw I tried with pari and the output parameter is of order of the input parameter : this means that you nearly need out of bound parameter at the input to get a problem !2012-08-16
  • 0
    @J.M. Try n ~ 10^9 and a ~ 10^(-30). Mathematica's LogGamma function appears to choke in this range.2012-08-16
  • 0
    @RogerS.: you are probably in 'numerical mode' (machine precision) use the 'infinite precision' mode instead (I don't have Mma so..)2012-08-16
  • 0
    @Raymond Manzoni, I've tried both modes...2012-08-16
  • 0
    @Roger S.: Define "choke". You can specify how many digits of precision you want mathematica to return. Calling mathematica through wolframalpha doesn't have any obvious problems: http://www.wolframalpha.com/input/?i=%2810%5E9-1%29+Log%2810%5E%28-30%29%29+%2B+LogGamma%2810%5E30%29+-+LogGamma%281-10%5E9%2B10%5E30%292012-08-16
  • 0
    @Roger: sorry in this case it seems really to be a bug in your version of Mma...2012-08-16
  • 0
    @RaymondManzoni, I updated the question with a specific example showing a problem with computing the approximation using LogGamma versus computing a direct product with output ~1.2012-08-16
  • 2
    @Roger: You're using `LogGamma[]` wrong. Recall the identity $\log(p/q)=\log\,p-\log\,q$...2012-08-16
  • 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]]]`2012-08-16