# Why do I get different results for the products of two identical expressions?

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

• 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. Commented Jun 10 at 17:08

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 *)


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 *)

• 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). Commented Jun 10 at 21:13
• 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. Commented Jun 10 at 21:59
• Can your Mathematica program be improved so that more decimals are correct? Commented Jun 10 at 22:29
• 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. Commented Jun 10 at 23:39

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
(*0100.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."

• 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! Commented Jun 11 at 16:38