> restart;
> Digits := 40;
> with(linalg);
> f1 := proc (n) local i, j, v; v := matrix(n+1, n+1); for i to n+1 do v[i, 1] := 1; for j from 2 to n+1 do v[i, j] := 1/(i+j-2) end do end do; RETURN(v) end proc;
> z := unapply(factorial(i-1)*exp(i+g-1/2)/((i+g-1/2)^(i-1/2)*sqrt(2*Pi)), i, g);
> f2 := proc (n, y) local i, v; v := matrix(n+1, 1); for i to n+1 do v[i, 1] := z(i, y) end do; RETURN(v) end proc;
> n := 7;
> g := 5;
> f := f1(n);
> g1 := inverse(f);
> b := f2(n, g);
> gg := multiply(g1, b);
> evalf(gg[1, 1]-1, 4*rhs(op(2, op(1, gg))[1])+4);
#output -9.3304068958898243168 e-12
> evalf(gg[rhs(op(2, op(1, gg))[1]), 1], 4*rhs(op(2, op(1, gg))[1])+4);
#output 0.000002394534913406850300243
> for i to rhs(op(2, op(1, gg))[1]) do evalf(gg[i, 1], 4*rhs(op(2, op(1, gg))[1])+2*i) end do;
#output
# 0.99999999999066959310411017566
# 76.1800917308668800992939735513657
# -86.5053203963967652303139216042357
# 24.01409899435588324058704081692617
# -1.231742921450034278218762025469663
# 0.00121555828611639057422612294568442
# -0.0000120262591451567138253281149246451
# 0.000002394534913406850300240770820124155
> f3 := evalf(gg[1, 1], 4*rhs(op(2, op(1, gg))[1])+1);
> for i from 2 to rhs(op(2, op(1, gg))[1]) do f3 := f3+evalf(gg[i, 1], 4*rhs(op(2, op(1, gg))[1])+i)/(x+i-1) end do;
> f4 := unapply((x+g+1/2)^(x+1/2)*exp(-x-g-1/2)*sqrt(2*Pi)*f3, x);
> evalf(f4(1/2));
#output
# 0.8862269254527589674644560128649165283764
> 4*(%*%)-evalf(Pi);
#output
# 6.762374918539830701792077 e-15
now compute the coefficients I cited.
> restart;
> Digits := 40;
> with(linalg);
> f1 := proc (n) local i, j, v; v := matrix(n+1, n+1); for i to n+1 do v[i, 1] := 1; for j from 2 to n+1 do v[i, j] := 1/(i+j-2) end do end do; RETURN(v) end proc;
> z := unapply(factorial(i-1)*exp(i+g-1/2)/((i+g-1/2)^(i-1/2)*sqrt(2*Pi)), i, g);
> f2 := proc (n, y) local i, v; v := matrix(n+1, 1); for i to n+1 do v[i, 1] := z(i, y) end do; RETURN(v) end proc;
> n := 14;
> g := 607/128;
> f := f1(n);
> g1 := inverse(f);
> b := f2(n, g);
> gg := multiply(g1, b);
> evalf(gg[1, 1]-1, 4*rhs(op(2, op(1, gg))[1])+4);
#output
# -2.90817953577301957618515490814322962380380 e-15
> evalf(gg[rhs(op(2, op(1, gg))[1]), 1], 4*rhs(op(2, op(1, gg))[1])+4);
#output
# 0.00000368991826595316227036759674456787362014192
> for i to rhs(op(2, op(1, gg))[1]) do evalf(gg[i, 1], 4*rhs(op(2, op(1, gg))[1])+2*i) end do;
#output
# 0.999999999999997091820464226980423814845091856770376328
# 57.156235665862923516579393485977444273254442892929107852
# -59.597960355475491248142266131602836001910100856036329460
# 14.1360979747417471738634195408767260156591641388059836898
# -0.4919138160976201997828400285303078441634827850503476487
# 0.00003399464998481188869891934155234155285736767075962214
# 0.000046523628927048575665230224965822889358967609346662905
# -0.00009837447530487956467653837063470730301913302874861670310
# 0.000158088703224912488836072413443663958858739583487234932827
# -0.00021026444172410488319269928300571190544166498661591073776001
# 0.0002174396181152126431961446496038346923763748760822190719415235
# -0.000164318106536763890217069562280184132420847097097457818988636042
# 0.00008441822398385274329281181534545498133890999303130609738370698039
# -0.0000261908384015814086696650362479549018414325910258739195513825336890
# 0.0000036899182659531622703675967445678736201430201881942403488492031328990
> f3 := evalf(gg[1, 1], 4*rhs(op(2, op(1, gg))[1])+1);
> for i from 2 to rhs(op(2, op(1, gg))[1]) do f3 := f3+evalf(gg[i, 1], 4*rhs(op(2, op(1, gg))[1])+i)/(x+i-1) end do;
> f4 := unapply((x+g+1/2)^(x+1/2)*exp(-x-g-1/2)*sqrt(2*Pi)*f3, x);
> evalf(f4(1/2));
#output
# 0.8862269254527580134320931844229516344461
# the % character is shorthand for last result
> 4*(%*%)-evalf(Pi);
#output
# -1.538422995214718383218 e-18
Bookmarks