I found a strange behavior in Mathematica when trying to evaluate the integral

$$f(n) = \int_{1}^2 \frac{\Gamma(n)\Gamma(x)}{\Gamma(x+n)}{\rm d}x$$

I evaluate this using F[n_] := NIntegrate[(Gamma[x] Gamma[n])/Gamma[x + n], {x, 1, 2}] For large values of $n$ the two terms $\Gamma(n)$ and $\Gamma(x+n)$ will be huge, but the integrand itself is monotonely decreasing and contained in $[0,\frac{1}{n}]$.

Mathematica ( is able to evaluate this integral for all values $n\in[1,1000]$ however for $n=171$ the integration breaks down with NIntegrate::izero and gives $0$ as the result:

F[170]  -> 0.00103244
F[171]  -> 0.0 (NIntegrate::izero)
F[172]  -> 0.00105424

If I try non-integer $n$ then I find that it breaks down for all $n \in [170.6,171.6]$ Increasing MinRecursion or AccuracyGoal does not help. The strange part is that Mathematica is able to evaluate the function perfectly well for $x\in[1,2]$ and the plot of the function over this interval looks just like it should. The plot below is $n\cdot\frac{\Gamma(n)\Gamma(x)}{\Gamma(x+n)}$ for $n=171$:

$~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~$enter image description here

My questions are as follows:

  • What is going on here?
  • Is there any settings that I can turn on to fix this?
  • $\begingroup$ Psst, Beta[] is built-in... $\endgroup$ Jun 3, 2016 at 17:12
  • $\begingroup$ @J.M. Yes that does works as a solution to this particular problem. I was wondering more in general how to solve such problems) and I still want to know what is going on here. It's really strange. $\endgroup$
    – Winther
    Jun 3, 2016 at 17:15
  • 5
    $\begingroup$ It's usually a bad idea to compute ratios of gamma functions (or factorials, for that matter) in this way for large arguments. Either you re-express everything in terms of LogGamma[] and exponentiate at the very end, or figure out an alternative expression using any of Pochhammer[], Binomial[], FactorialPower[], Beta[]... $\endgroup$ Jun 3, 2016 at 17:17
  • $\begingroup$ @J.M. Thanks I'll keep that in mind. $\endgroup$
    – Winther
    Jun 3, 2016 at 17:19

2 Answers 2


Fixing it ...

Is there any settings that I can turn on to fix this?

Does using larger working precision produce results you expect?

F[n_] := 
 NIntegrate[(Gamma[x] Gamma[n])/Gamma[x + n], {x, 1, 2}, 
  WorkingPrecision -> 60, MinRecursion -> 6, MaxRecursion -> 20]

In[10]:= F[170]

Out[10]= 0.00106892313424594316124137429342057472552302576793438553998338

In[11]:= F[171]

Out[11]= 0.00106153418172547954471366095844979404173963460222797906864894

In[12]:= F[172]

Out[12]= 0.00105424004229628467250748022777662650942896871502716627259853


Following J.M.'s advice of using Beta we can see that using larger working precision fixes the issue when using Gamma:

FB[n_] := 
 NIntegrate[Beta[x, n], {x, 1, 2}, WorkingPrecision -> 60, 
  MinRecursion -> 6, MaxRecursion -> 20]

Grid[{#, FB[#]} & /@ Range[170, 172]]

enter image description here

What is going on ...

Following george2079's answer and using the option "IntegrationMonitor" we can see what is going on.

With[{n = 171}, 
 Reap@NIntegrate[(Gamma[x] Gamma[n])/Gamma[x + n], {x, 1, 2}, 
   IntegrationMonitor -> (Sow[
       Map[{#1@"Integrand", #1@"Boundaries", #1@"Integral", #1@
           "Error", #1@"GetValues"} &, #1]] &)]

enter image description here

f[n_?NumericQ, x_?NumericQ] := Gamma[x] Gamma[n]/Gamma[x + n];
With[{n = 171},
 Reap@NIntegrate[f[n, x], {x, 1, 2}, 
   IntegrationMonitor -> (Sow[
       Map[{#1@"Integrand", #1@"Boundaries", #1@"Integral", #1@
           "Error", #1@"GetValues"} &, #1]] &)]

enter image description here

  • $\begingroup$ It's easy to check: use Beta[x, n] instead. $\endgroup$ Jun 3, 2016 at 17:46
  • $\begingroup$ Thanks for your answer and the interesting features you mention that I was not aware of like IntegrationMonitor. WorkingPrecision does fix it, should have though of that! $\endgroup$
    – Winther
    Jun 4, 2016 at 1:03
  • $\begingroup$ @Winther Great then! $\endgroup$ Jun 4, 2016 at 1:06

A bit speculative, but I think we can see why 171 is a magic number, if we factor the Gamma[n] from the integral:

 Table[ {n, NIntegrate[Gamma[x] /Gamma[x + n], {x, 1, 2}]} , {n, 168, 173}]

{{168, 7.20924*10^-304}, {169, 4.26119*10^-306}, {170, 2.41873*10^-308}, {171, 0.}, {172, 0.}, {173, 0.}}

171 is where the integral fails as you can see the result is smaller than $MinMachineNumber. I'd guess NIntegrate factors the constant, but stops factoring when the factored value is too large.

note if we hide the structure the problem ( at this n ) goes away:

f[n_?NumericQ, x_?NumericQ] := Gamma[x] Gamma[n]/Gamma[x + n]
Table[{n, NIntegrate[f[n, x], {x, 1, 2}]}, {n, 168, 175}]

{{168, 0.00108399}, {169, 0.00107641}, {170, 0.00106892}, {171, 0.00106153}, {172, 0.00105424}, {173, 0.00104704}, {174, 0.00103993}, {175, 0.00103291}}

This actually (surprisingly) works even for much larger n.

  • $\begingroup$ Thank you for your answer. That sounds like a plausible reason. Also nice to learn about the use of NumericQ in the variable definition, was unaware of that you could do that. $\endgroup$
    – Winther
    Jun 4, 2016 at 1:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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