6
$\begingroup$

I have two expressions with the Gamma function that are identical:

 FullSimplify[Gamma[2*k/(2*k - 1)]/Gamma[1 + 1/(2*k)] == 
 2*k/(2*k - 1)*Gamma[1/(2*k - 1)]/Gamma[1/(2*k)]]
 (* True *)

But when using NProduct I get completely different results:

 NProduct[Gamma[2*k/(2*k - 1)]/Gamma[1 + 1/(2*k)], {k, 1, Infinity}]
 (* 1.06215 *)

 NProduct[2*k/(2*k - 1)*Gamma[1/(2*k - 1)]/Gamma[1/(2*k)], {k, 1, Infinity}]
 (* ComplexInfinity *)

The product is convergent, so this result is wrong. I don't know if this qualifies as a BUG, but it would be better to give some sort of error message than ComplexInfinity for a convergent product. This is very confusing.

If I add WorkingPrecision, the program does give error messages, but not infinity

 NProduct[Gamma[2*k/(2*k - 1)]/Gamma[1 + 1/(2*k)], {k, 1, Infinity}, WorkingPrecision -> 20]
 (* NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections... *)
 (* 1.0621509056004 *)

I can also add that Maple calculates this product to 100 decimal places without any problems.

 evalf(Product(2*k*GAMMA(1/(2*k - 1))/((2*k - 1)*GAMMA(1/(2*k))), k = 1..infinity), 120);
 1.06215090557105728069683736293809990425207955200456933340798700905379893707714082919361825368669317760219700551663195491
Improve this question
$\endgroup$
1
  • $\begingroup$ I think the difference probably comes from {Limit[Gamma[2*k/(2*k - 1)], k -> Infinity], Limit[Gamma[1 + 1/(2*k)], k -> Infinity]} vs. {Limit[Gamma[1/(2*k - 1)], k -> Infinity], Limit[Gamma[1/(2*k)], k -> Infinity]} -- The singular forms are probably (mis)handled in a different way, exposing a bug. $\endgroup$
    – Michael E2
    Commented Jun 10 at 17:08

3 Answers 3

Reset to default
5
$\begingroup$

Bug in NProduct maybe. For that second one, rote conversion to logs, then NSum, then Exp works fine.

Exp@NSum[Log[2*k] - Log[(2*k - 1)] + LogGamma[1/(2*k - 1)] - 
   LogGamma[1/(2*k)], {k, 1, Infinity}]

(* Out[196]= 1.06215 *)
Improve this answer
$\endgroup$
4
$\begingroup$

I'd call it a bug. Here's a way to get 100 digits. If correct, it shows that the Maple result is not accurate to all digits shown. Note that NSum[] on the logarithm of the factors apparently has the same bug as the OP's NProduct[]. (Update, to clarify the difference with @Daniel's answer: I took the log but didn't expand, which I thought to do later to respond to the Maple result: NSum[2*k/(2*k - 1)*Gamma[1/(2*k - 1)]/Gamma[1/(2*k)] // Log, {k, 1, Infinity}] --> ComplexInfinity seems a bug, too. Machine precision Exp@NSum[... // Log // PowerExpand, {k, 1, Infinity}] works correctly, even without the log-gamma replacement.)

opts = Options@NIntegrate;

SetOptions[NIntegrate, MaxRecursion -> 100, 
  Method -> {"GaussKronrodRule", "Points" -> 21}, MinRecursion -> 6];
prod100 = Exp@N[Sum[
    2*k/(2*k - 1)*Gamma[1/(2*k - 1)]/Gamma[1/(2*k)] // Log // 
      PowerExpand // ReplaceAll[Log[Gamma[x_]] :> LogGamma[x]], {k, 1,
      Infinity}], 100]
SetOptions[NIntegrate, opts];
(*
1.06215090557105728069683736293741464202343947470947209825644524241677123960958822156592210084494367818
*)

SetOptions[NIntegrate, MaxRecursion -> 100, 
  Method -> {"GaussKronrodRule", "Points" -> 41}, MinRecursion -> 9];
prod150 = 
 Exp@N[Sum[
    2*k/(2*k - 1)*Gamma[1/(2*k - 1)]/Gamma[1/(2*k)] // Log // 
      PowerExpand // ReplaceAll[Log[Gamma[x_]] :> LogGamma[x]], {k, 1,
      Infinity}], 150]
SetOptions[NIntegrate, opts];
(*
1.0621509055710572806968373629374146420234394747094720982564452424167712396095882215659221008449436781834466909882566217898444265692749546388048063861684
*)

Error:

prod100 - prod150

(*  0.*10^-102  -- 102 digit agreement *)

Maple:

maple = 1.06215090557105728069683736293809990425207955200456933340798700905379893707714082919361825368669317760219700551663195491;
maple - prod150 // N

(*  6.85262*10^-31  -- 30 digit agreement *)
Improve this answer
$\endgroup$
4
  • $\begingroup$ This is very interesting! We have reliably determined digits 1.06215090557105728069683736293... I did the same experiment in Maple, and if I compare the results there to 100 and 400 decimal places, I also get the first 100 digits the same (but different from Mathematica). $\endgroup$
    – Vaclav Kotesovec
    Commented Jun 10 at 21:13
  • $\begingroup$ Maple result is correct (all digits!). I used PARI default(realprecision, 200); exp(sumpos(k=1, log(2*k) + log(gamma(1/(2*k-1))) - log(2*k-1) - log(gamma(1/(2*k))) )) with exactly the same result. $\endgroup$
    – Vaclav Kotesovec
    Commented Jun 10 at 21:59
  • $\begingroup$ Can your Mathematica program be improved so that more decimals are correct? $\endgroup$
    – Vaclav Kotesovec
    Commented Jun 10 at 22:29
  • $\begingroup$ It seems Mathematica and Maple are converging to a different numbers as the precision increases (I took it out to 400, too, and the digits always agree up to the lower precision). I don't know why or how to fix it, if Maple and PARI are correct. "My" program is just using N[] on Sum[], which calls Mathematica's program for approximating sums to arbitrary precision, plus adding power to NIntegrate. You could report it to Wolfram and see if they agree that their program is in error. $\endgroup$
    – Michael E2
    Commented Jun 10 at 23:39
1
$\begingroup$

Based on @Michael E2's answer, I modified some options and obtained results that are consistent with Maple's.

opts = Options@NIntegrate;
SetOptions[NIntegrate, MaxRecursion -> 100, MinRecursion -> 15, 
  PrecisionGoal -> 100, AccuracyGoal -> 100, 
  WorkingPrecision -> 200];
res = NProduct[
  Gamma[2*k/(2*k - 1)]/Gamma[1 + 1/(2*k)], {k, 1, Infinity}, 
  Method -> "EulerMaclaurin", PrecisionGoal -> 100, 
  AccuracyGoal -> 100, WorkingPrecision -> 200, 
  NProductFactors -> 50]
SetOptions[NIntegrate, opts];
maple = 1.\
0621509055710572806968373629380999042520795520045693334079870090537989\
3707714082919361825368669317760219700551663195491
res - maple
(*0``100.03958308889237*)

In fact, if you specify Method as "EulerMaclaurin" and do not modify NProductFactors, Mathematica will give a result with 30-digit precision and the message "Euler-Maclaurin sum failed to converge to requested error tolerance."

Improve this answer
$\endgroup$
1
  • $\begingroup$ Great! Too bad Mathematica doesn't tune the parameters itself. In the case of Maple or PARI, I don't have to specify anything like that, just the desired precision. Thank you! $\endgroup$
    – Vaclav Kotesovec
    Commented Jun 11 at 16:38

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.