On Sat, Sep 1, 2012 at 11:47 AM, Bill Gosper <billgosper@gmail.com> wrote:
Stephen Finch has given this "Constant of Theodorus" (good name, but not to be confused with Mathworld's Theodorus's Constant = √3) his usual superb treatment in http://www.people.fas.harvard.edu/~sfinch/csolve/th.pdf , but the rapid numerical method I give below may be new.
It seems to me this umbral notation trick should apply more generally to the Euler-Maclaurin formula (perhaps even accommodating my prescription for Bernpolys over plain Berns), streamlining both notation and computation.
For Stirling's formula, this works unreasonably well: Out[817]= Log[k!] ->(Log[2 π]-1)/2 + (-B + k + y) Log[-B + k + y] - k (The constant term coming from k! == (2^(-2 k) Sqrt[π] (2 k)!)/(k-1/2)! ) In[818]:= Normal[Series[#,{B,0,5}]&/@%817]/.B^(p_:1)->BernoulliB[p,y] Out[818]= Log[k!]->-(1/2)-k+(1/6-y+y^2)/(2 (k+y))+(y/2-(3 y^2)/2+y^3)/(6 (k+y)^2)+(-(1/30)+y^2-2 y^3+y^4)/(12 (k+y)^3)+(-(y/6)+(5 y^3)/3-(5 y^4)/2+y^5)/(20 (k+y)^4)+1/2 Log[2 π]+(-(1/2)+y) (-1-Log[k+y])+k Log[k+y]+y Log[k+y] In[819]:= TableForm[Table[N[%],{k,1,3},{y,0,1,1/3}]] Out[819]//TableForm= 0.->-0.000505911 0.->5.91272*10^-6 0.->0.0000436355 0.->-0.0000212515 0.693147->0.693126 0.693147->0.693149 0.693147->0.693151 0.693147->0.693144 1.79176->1.79176 1.79176->1.79176 1.79176->1.79176 1.79176->1.79176 In[820]:= TableForm[Table[N[%%,9],{k,3,5},{y,0,1,1/2}]]
Out[820]//TableForm=
1.79175947->1.79175644 1.79175947->1.79176085 1.79175947->1.79175873 3.17805383->3.17805309 3.17805383->3.17805423 3.17805383->3.17805358 4.78749174->4.78749150 4.78749174->4.78749189 4.78749174->4.78749164
Probably even better with y some simple function of k. --rwg