(* notebook for Master class given by Kay Wiese, (C) 2020 *)
(* we wish to resum the following series, using Padé-Borel resummation *)
f[g]=(∑_(n=0)^(Subscript[n, max]))(((-g))^n) (Subscript[a, n]);
(* the coefficients are *)
an[n_]:=(4n)!/((2n)!)/2^(2n)/n!

(* our own integration routine, allows to cut out the singularities *)
Clear[OwnInt];
workingPrecision=20;precisionGoal=12;
OwnInt[f_,tmax_,g_,posroots_,δ_,verbatim_:False]:=Module[{X,points},
LVPrint[x___]:=If[verbatim,Print[x]];
Clear[t];
nroots=Length[posroots];
Off[NIntegrate::slwcon];
If[verbatim,Print[nroots," positive roots located at ",posroots]];
If[nroots==0,Return[ NIntegrate[f[t g]Exp[-t],{t,0,tmax}] ]];
If[tmax< Max[posroots],Print["KayInt[] error: tmax< Max[nroots]"]];
points={posroots,tmax}//Flatten//Sort;
If[verbatim,Print["points = ",points]];
X=NIntegrate[Re[f[t g]]Exp[-t],{t,0,points[[1]]-δ},WorkingPrecision->workingPrecision,PrecisionGoal->precisionGoal,Exclusions->posroots];
LVPrint["first integral = ",X];
For[i=1,i≤nroots,i++,
LVPrint["i=",i];
Y=NIntegrate[Re[f[t g]]Exp[-t],{t,points[[i]]+δ,points[[i+1]]-δ},Exclusions->posroots,WorkingPrecision->workingPrecision,PrecisionGoal->precisionGoal];
LVPrint["Int[f(tg)(e^(-t)),{t,",points[[i]]+δ,",",points[[i+1]]-δ,"}=",Y];
X+=Y;
];
On[NIntegrate::slwcon];
Return[X];
]

Pverbatim=True;
Fverbatim=False;
nmax=6;δ=10^(-4);tmax=10;
origser=Sum[an[n]/n!(-t)^n,{n,0,nmax}];
Clear[pade,Iexact]
Iexact[g_]=Exp[1/(32g)]BesselK[1/4,1/(32 g)]/4/Sqrt[Pi g];
Pexact=Plot[Iexact[g],{g,0,1},PlotStyle->{Green,Dashed,Thick}];
For[r=0,r≤ nmax,r++,
Print["************************************************************"];
Print["Padé ",{nmax-r,r}];
Clear[pade];
pade[t_]=PadeApproximant[origser,{t,0,{nmax-r,r}}];
If[Fverbatim,Print["Padé=",pade[t]]];
denom=Denominator[pade[t]];
numer=Numerator[pade[t]];
If[Fverbatim,Print["denom = ",denom]];
roots=t/.NSolve[denom==0,t];
zeroes=t/.NSolve[numer==0,t];
If[roots==t,roots={}];
If[zeroes==t,zeroes={}];
If[Pverbatim,Show[{ListPlot[zeroes//ReIm,PlotStyle->Blue,AxesOrigin->{0,0},PlotRange->All],ListPlot[roots//ReIm,PlotStyle->Red,AxesOrigin->{0,0},PlotRange->All]},PlotRange->All]//Print];
Print["roots = ",roots];
posroots=Select[roots,tmax>Re[#]>0&&Im[#]==0&];nroots=Length[posroots];
Print["positive roots = ",posroots];
P1=Plot[pade,{t,0,5}];
tab=Table[{g,OwnInt[pade[#]&,tmax,g,posroots,δ,False]},{g,0,1,0.01}];
P2={Pexact,ListPlot[tab,PlotStyle->Red],ListLinePlot[tab,PlotStyle->Blue]};
If[Pverbatim,P2//Show//Print]
]

(* some checks on routines *)
pade[t]
posroots
Plot[pade[t],{t,0,0.3}]
OwnInt[pade[#]&,10,0.1,posroots,0.1,True]