[math-fun] I just ported my ancient Macsyma Euler-Maclaurin fcn to Mma.
(With exact remainder and optional phase parameter.) EulMac[fn,var lo,hi,trms,phase] Here are three terms of Stirling's formula: In[305]:= Assuming[n > 1, EulMac[Log[z], z, 1, n, 3]] Out[305]= 1 - n + n*Log[n] - Inactive[Sum][Log[z], {z, 2, n}] == (-(1/24))*Inactive[Sum][ ConditionalExpression[ (-1 + 3*¢*(1 + ¢)* (-1 + 10*¢*(1 + ¢)* (1 + 12*¢*(1 + ¢))) + 180*¢^3*(1 + ¢)^3* (1 + 2*¢)*(Log[¢] - Log[1 + ¢]))/(15*¢^3*(1 + ¢)^3), Re[¢] > 0 || Re[¢] < -1 || NotElement[¢, Reals]], {¢, 1, -1 + n}] + Inactive[Sum][((-1)^¢*BernoulliB[1 + ¢]* Derivative[¢][-Log[1 + #1] + Log[n + #1] & ][0])/ (1 + ¢)!, {¢, 0, 3}] The self-evidently superfluous ConditionalExpression surrounds the remainder summand. It is typically textually cumbersome but numerically small. The Sums are inactive to discourage Mathematica from struggling to "simplify" them. Activating gives Log[Pochhammer] instead of Log[Factoria] if you don't call n a positive integer. In[311]:= Assuming[n > 1 && Element[n,Integers], FullSimplify /@ Activate[Out[305]]]; In[314]:= Out[311] /. ConditionalExpression[ x_, __] :> FullSimplify[x]; In[315]:= Solve[%, Log[n!]] Out[315]= {{Log[n!] -> (1/(360*n^3))* (-1 + 30*n^2 + 331*n^3 - 360*n^4 + 180*n^3*Log[n] + 360*n^4*Log[n] + 15*n^3*Sum[(1/(15*¢^3* (1 + ¢)^3))*(-1 + 3*¢*(1 + ¢)*(-1 + 10*¢*(1 + ¢)* (1 + 12*¢*(1 + ¢))) + 180*¢^3*(1 + ¢)^3* (1 + 2*¢)*(Log[¢] - Log[1 + ¢])), {¢, 1, -1 + n}])}} In[318]:= Expand /@ %315[[1, 1]] Out[318]= Log[n!] -> 331/360 - 1/(360*n^3) + 1/(12*n) - n + Log[n]/2 + n*Log[n] + (1/24)* Sum[(1/(15*¢^3*(1 + ¢)^3))* (-1 + 3*¢*(1 + ¢)* (-1 + 10*¢*(1 + ¢)* (1 + 12*¢*(1 + ¢))) + 180*¢^3*(1 + ¢)^3* (1 + 2*¢)*(Log[¢] - Log[1 + ¢])), {¢, 1, -1 + n}] Note 331/360 vs ln (rt2 pi): In[44]:= N[{Log[Sqrt[2*Pi]], 331/360}] Out[44]= {0.918938533, 0.919444444} Testing on 8! In[319]:= %318 /. n -> 8; In[320]:= Expand /@ % Out[320]= Log[40320] -> Log[2] + Log[3] + Log[4] + Log[5] + Log[6] + Log[7] + Log[8] (because the remainder sum is exact). Punting it, In[321]:= %318 /. Sum[__] -> 0 Out[321]= Log[n!] -> 331/360 - 1/(360 n^3) + 1/(12 n) - n + Log[n]/2 + n Log[n] Testing with n=8, In[322]:= % /. n -> 8 Out[322]= Log[40320] -> -(1303169/184320) + (17 Log[8])/2 In[323]:= N[%] Out[323]= 10.6046029 -> 10.6051088 Repeating the analysis starting with Assuming[n > 1, EulMac[Log[z], z, 1, n, 5, 1/2]] (i.e., one more Bernoulli term and phase = 1/2) Out[46]= Log[n!] -> 628991/612360 - 31/(1260 (1 + 2 n)^5) + 7/(360 (1 + 2 n)^3) - 1/( 12 + 24 n) + Log[2] - n (1 + Log[2]) - (3 Log[3])/ 2 + (1/2 + n) Log[1 + 2 n] In[47]:= % /. n -> 8. Out[47]= 10.6046029 -> 10.6046256 with the ^-5 term contributing only In[48]:= %%[[2, 2]] /. {{n -> 8}, {n -> 8.}} Out[48]= {-31/1789019820, -1.73279243*10^-8} so most of the .00002 error resides in the (discarded) remainder term: Out[49]= -433349316889/62104545180 - 8 Log[2] - 5 Log[3]/2 - Log[4] - Log[5] - Log[6] - Log[7] - Log[8] + 17 Log[17]/2 In[50]:= N[%] Out[50]= 0.0000226652532 The function is EulMac[xp_, v_, lo_, hi_, deg_, y_: 0] := Integrate[xp /. v -> y + v, {v, lo, hi}] - Inactive[Sum][xp, {v, lo + 1, hi}] == (-1)^deg*Inactive[Sum][ Integrate[(deg + 1)*v^deg* D[xp /. v -> v*y + ¢ + 1, {v, deg + 1}] - BernoulliB[deg + 1, y + v]* D[xp /. v :> y + v + ¢, {v, deg + 1}], {v, 0, 1}], {¢, lo, hi - 1}]/(deg + 1)! + Inactive[Sum][(-1)^¢* Derivative[¢][Evaluate[(xp /. v -> # + hi) - (xp /. v -> # + lo)] &][y]* BernoulliB[¢ + 1, y]/(¢ + 1)!, {¢, 0, deg}] The motivation for porting was NSum's seeming inability to get so much as four digits of Sum 2 to oo of 1/Log[n!] - 1/Log[(2 n)!] - 1/Log[(-1 + 2 n)!] apparently due to NSum's internal Euler-Maclaurin failing to try Method -> "DoubleExponential" on the corresponding NIntegrate. I now believe the sum ~ 1.4424630245831935 . This is still not accurate enough for PSLQ or the like, and ISC humorously proposes 1/Log[2 + 1/(300 Khinchin Gamma[1/6])]. --rwg
On Sep 15, 2015, at 5:04 PM, Bill Gosper <billgosper@gmail.com> wrote:
Out[305]= 1 - n + n*Log[n] - Inactive[Sum][Log[z], {z, 2, n}] == (-(1/24))*Inactive[Sum][ ConditionalExpression[ (-1 + 3*¢*(1 + ¢)* (-1 + 10*¢*(1 + ¢)* (1 + 12*¢*(1 + ¢))) + 180*¢^3*(1 + ¢)^3* (1 + 2*¢)*(Log[¢] - Log[1 + ¢]))/(15*¢^3*(1 + ¢)^3), Re[¢] > 0 || Re[¢] < -1 || NotElement[¢, Reals]], {¢, 1, -1 + n}] + Inactive[Sum][((-1)^¢*BernoulliB[1 + ¢]* Derivative[¢][-Log[1 + #1] + Log[n + #1] & ][0])/ (1 + ¢)!, {¢, 0, 3}] I see you have gone beyond asserting “my two cents” by summing to ¢=3.
-Veit
On 2015-09-15 14:25, Veit Elser wrote:
On Sep 15, 2015, at 5:04 PM, Bill Gosper <billgosper@gmail.com> wrote:
Out[305]= 1 - n + n*Log[n] - Inactive[Sum][Log[z], {z, 2, n}] == (-(1/24))*Inactive[Sum][ ConditionalExpression[ (-1 + 3*¢*(1 + ¢)* (-1 + 10*¢*(1 + ¢)* (1 + 12*¢*(1 + ¢))) + 180*¢^3*(1 + ¢)^3* (1 + 2*¢)*(Log[¢] - Log[1 + ¢]))/(15*¢^3*(1 + ¢)^3), Re[¢] > 0 || Re[¢] < -1 || NotElement[¢, Reals]], {¢, 1, -1 + n}] + Inactive[Sum][((-1)^¢*BernoulliB[1 + ¢]* Derivative[¢][-Log[1 + #1] + Log[n + #1] & ][0])/ (1 + ¢)!, {¢, 0, 3}] I see you have gone beyond asserting “my two cents” by summing to ¢=3.
-Veit
At least I kept it real. In the 2¢ plane you can run into copyright problems. [big chop] The function is
EulMac[xp_, v_, lo_, hi_, deg_, y_: 0] := Integrate[xp /. v -> y + v, {v, lo, hi}] - Inactive[Sum][xp, {v, lo + 1, hi}] == (-1)^deg*Inactive[Sum][ Integrate[(deg + 1)*v^deg* D[xp /. v -> v*y + ¢ + 1, {v, deg + 1}] - BernoulliB[deg + 1, y + v]* D[xp /. v :> y + v + ¢, {v, deg + 1}], {v, 0, 1}], {¢, lo, hi - 1}]/(deg + 1)! + Inactive[Sum][(-1)^¢* Derivative[¢][Evaluate[(xp /. v -> # + hi) - (xp /. v -> # + lo)] &][y]* BernoulliB[¢ + 1, y]/(¢ + 1)!, {¢, 0, deg}]
The motivation for porting was NSum's seeming inability to get so much as four digits of Sum 2 to oo of
1/Log[n!] - 1/Log[(2 n)!] - 1/Log[(-1 + 2 n)!]
apparently due to NSum's internal Euler-Maclaurin failing to try Method -> "DoubleExponential" on the corresponding NIntegrate. I now believe the sum ~ 1.4424630245831935 .
It may be merely 1/ln2 ~ 1.44269504, depending on how gradually f(n) can decrease before Limit[-f(n+2)-f(n+3)-...-f(2n+1),n->∞] = Sum(n>0)(f(n+1) - f(2n) - f(2n+1)) is no longer 0. Can this limit vanish w/o (non-oscillatory) f converging? E.g., can we construct an f where Limit[-f(n+2)-f(n+3)-...-f(2n+1),n->∞] =0 but Limit[-f(n+2)-f(n+3)-...-f(3n+1),n->∞] >0 ? Empirically, the sum seems to vanish for f := 1/n/ln(n+a), even though Sum f diverges. Or maybe it's just very small: < 6e-15 for a=1. Can we lower this bound? For the problem at hand, f(n):=1/Log[(n+1)!], which converges only a little slower than 1/n/ln(n+1), but the smallest magnitude EulMac and NSum have found is -.0002, with a supposedly more careful result of -.00035. Can anyone help with these sums? --rwg
This is still not accurate enough for PSLQ [n]or the like, and ISC humorously proposes 1/Log[2 + 1/(300 Khinchin Gamma[1/6])]. --rwg
On Wed, Sep 16, 2015 at 1:55 PM, Bill Gosper <billgosper@gmail.com> wrote:
On 2015-09-15 14:25, Veit Elser wrote:
On Sep 15, 2015, at 5:04 PM, Bill Gosper <billgosper@gmail.com> wrote:
Out[305]= 1 - n + n*Log[n] - Inactive[Sum][Log[z], {z, 2, n}] == (-(1/24))*Inactive[Sum][ ConditionalExpression[ (-1 + 3*¢*(1 + ¢)* (-1 + 10*¢*(1 + ¢)* (1 + 12*¢*(1 + ¢))) + 180*¢^3*(1 + ¢)^3* (1 + 2*¢)*(Log[¢] - Log[1 + ¢]))/(15*¢^3*(1 + ¢)^3), Re[¢] > 0 || Re[¢] < -1 || NotElement[¢, Reals]], {¢, 1, -1 + n}] + Inactive[Sum][((-1)^¢*BernoulliB[1 + ¢]* Derivative[¢][-Log[1 + #1] + Log[n + #1] & ][0])/ (1 + ¢)!, {¢, 0, 3}] I see you have gone beyond asserting “my two cents” by summing to ¢=3.
-Veit
At least I kept it real. In the 2¢ plane you can run into copyright problems.
[big chop]
The function is
EulMac[xp_, v_, lo_, hi_, deg_, y_: 0] := Integrate[xp /. v -> y + v, {v, lo, hi}] - Inactive[Sum][xp, {v, lo + 1, hi}] == (-1)^deg*Inactive[Sum][ Integrate[(deg + 1)*v^deg* D[xp /. v -> v*y + ¢ + 1, {v, deg + 1}] - BernoulliB[deg + 1, y + v]* D[xp /. v :> y + v + ¢, {v, deg + 1}], {v, 0, 1}], {¢, lo, hi - 1}]/(deg + 1)! + Inactive[Sum][(-1)^¢* Derivative[¢][Evaluate[(xp /. v -> # + hi) - (xp /. v -> # + lo)] &][y]* BernoulliB[¢ + 1, y]/(¢ + 1)!, {¢, 0, deg}]
The motivation for porting was NSum's seeming inability to get so much as four digits of Sum 2 to oo of
1/Log[n!] - 1/Log[(2 n)!] - 1/Log[(-1 + 2 n)!]
apparently due to NSum's internal Euler-Maclaurin failing to try Method -> "DoubleExponential" on the corresponding NIntegrate. I now believe the sum ~ 1.4424630245831935 .
It may be merely 1/ln2 ~ 1.44269504, depending on how gradually f(n) can decrease before Limit[-f(n+2)-f(n+3)-...-f(2n+1),n->∞] = Sum(n>0)(f(n+1) - f(2n) - f(2n+1)) is no longer 0.
Julian to the rescue. With permission: JJZH>It can decrease at any rate faster than 1/n. Specifically, if f(n)=o(1/n) then the limit is zero, while if 1/n=O(f(n)) then the limit is non-zero. So yes, it is 1/ln2
Can this limit vanish w/o (non-oscillatory) f converging? E.g., can we construct an f where Limit[-f(n+2)-f(n+3)-...-f(2n+1),n->∞] =0 but Limit[-f(n+2)-f(n+3)-...-f(3n+1),n->∞] >0 ?
I'm not sure what the more general question is. It's easy to get Limit[-f(n+2)-f(n+3)-...-f(2n+1),n->∞] =0 with Sum[f[n],{n,∞}]=∞ — for f[n]=1/n/log(n+a), the limand is ~ln2/log n. You can't get the requested Limit[-f(n+2)-f(n+3)-...-f(2n+1),n->∞] =0 but Limit[-f(n+2)-f(n+3)-...-f(3n+1),n->∞] >0 with f always positive, since the sum of f over [n,3n] is covered by the sums over [n,2n] and [2n,4n], but f=1/n log n has the sum over [n,2n] tending to zero but the sum over [n,n^2] tending to ln2. With more slowly decreasing functions you could probably make [n,2n] zero but [n,n log n] positive. Julian Empirically, the sum seems to vanish for f := 1/n/ln(n+a), even though Sum f diverges. Or maybe it's just very small: < 6e-15 for a=1. Can we lower this bound? For the problem at hand, f(n):=1/Log[(n+1)!], which converges only a little slower than 1/n/ln(n+1), but the smallest magnitude EulMac and NSum have found is -.0002, with a supposedly more careful result of -.00035. Can anyone help with these sums? --rwg Tnx, Julian! --rwg
This is still not accurate enough for PSLQ [n]or the like, and ISC humorously proposes 1/Log[2 + 1/(300 Khinchin Gamma[1/6])]. --rwg
participants (2)
-
Bill Gosper -
Veit Elser