(* SetDirectory["~/math/screened"]; *)
ndegree=2; (* degree of algebraic approximant *)
nm=103 (*103 - max*);
ndig=128;
gmin=1/100;
gmax=3;
gstep=1/100;
xmin=0;
xmax=3;
xstep=1/100;
nen=2; (* only the last nen approximants will be calculated *)
<<series.dat;
os=OpenWrite["sum.dat"];
eng[1]=-1/2;
Do[eng[n]=en[n-1],{n,2,nm}];
Do[
Do[en[n]=eng[n]/.g->lambda,{n,nm}];
en0=en[1];
ndegall={ndegree};
ndegm=Length[ndegall];
      Clear[cf,ene];
e1=Table[Null,{n,2,nm}];
a0=en[1];
Do[e1[[n-1]]=en[n],{n,2,nm}];
  Clear[en];
km0=2;
kmm=ndegall[[ndegm]]+1;
ek=Table[Null,{k,1,kmm},{n,1,nm}];
c=Table[Null,{k,1,kmm},{n,1,nm}];
r=Table[Null,{k,1,kmm},{k,1,kmm}];
   Do[ndegree=ndegall[[ndegn]];
Print["--------- degree = ",ndegree];
km=ndegree+1;
ek[[1,1]]=1;
Do[ek[[1,n]]=0, {n,2,nm}];
Do[ek[[2,n]]=e1[[n-1]], {n,2,nm}];
If[km<3,Goto[1]];
Do[Print["Evaluating E^",k-1];
   Do[ek[[k,n]]=Sum[e1[[n-l]]*ek[[k-1,l]],{l,k-1,n-1}],
      {n,k,nm}],{k,km0,km}];
Label[1];km0=km+1;
Do[Do[r[[k,n]]=(-a0)^(n-k)*Binomial[n-1,k-1],
      {n,1,km}],{k,1,km}];
Do[Do[c[[k,n]]=ek[[k,n]],{n,k,nm}],{k,1,km}];
Do[nf=ka-km;
   nf2=nf+2;
   Do[n=nm-n5+nf2;
      c[[1,n]]=c[[1,n-1]],{n5,nf2,nm}];
   c[[1,nf+1]]=1;
   Do[n0=nf+k;
      a=c[[1,n0]]/c[[k,n0]];
      c[[k,nf+1]]=-a;
      cf[k,nf+1]=N[c[[k,nf+1]],ndig];
      Do[c[[1,n]]=c[[1,n]]-a*c[[k,n]], {n,n0,nm}], {k,2,km}];
   Do[a=c[[1,n]];
      Do[c[[k-1,n]]=c[[k,n]], {k,2,km}];
   c[[km,n]]=a, {n,nf2,nm}], {ka,km,nm}]
,{ndegn,1,ndegm}];
Clear[e1,ek,c,r];
acc=3; (* accuracy limit *)
roots=Table[Null,{n,1,nen},{nd,1,ndegree}];
r=Table[Null,{k,0,ndegree},{k,0,ndegree}];
energy=Table[Null,{n,1,nen}];
delta=Table[Null,{n,1,nen}];
km=ndegree+1;
Clear[p,t];
rs=Roots[Sum[p[n]*t^n,{n,0,ndegree}]==0,t,Cubics->True,Quartics->True];
Do[
Do[Do[r[[k,n]]=(-a0)^(n-k)*Binomial[n-1,k-1],
      {n,1,km}],{k,1,km}];
x=SetPrecision[Delta,ndig];
Do[n=ka;
   na=n-km+1;
   cf[1,na]=1;
   Do[a=x*cf[1,na]*r[[k,1]];
      Do[a=a+cf[m,na]*r[[k,m]],{m,2,km}];
      Do[r[[k,m-1]]=r[[k,m]],{m,2,km}];
      p[k-1]=r[[k,km]]=a, {k,1,km}];
   If[n<nm-nen+1,Goto[1]];
   rs0=N[rs,ndig];
   s0=Solve[rs0,t];
   s=N[s0];
   Do[roots[[nm-n+1,k]]=Replace[t,s[[k]][[1]]],{k,1,ndegree}];
   Label[1], {ka,km,nm}];
d=10^66;
Do[Do[d1=Abs[roots[[1,k1]]-roots[[2,k2]]];
   If[d1>d,Goto[3]];
   d=d1;
   k10=k1;k20=k2;
   Label[3],{k1,1,ndegree}],{k2,1,ndegree}];
y0=roots[[1,k10]];
Do[d=10^66;
  Do[d1=Abs[roots[[n,k]]-y0];
     If[d1>d,Goto[4]];
     k0=k;
     d=d1;
     Label[4],{k,1,ndegree}];
   energy[[n]]=roots[[n,k0]];
   delta[[n]]=d,{n,1,nen}];
d=0;
Do[d1=delta[[n]];
   If[d1<d,Goto[5]];
   d=d1;
   Label[5],{n,1,nen}];
(*If[d>acc,Exit[]];*)
yre=SetPrecision[Re[y0],12];
yim0=Abs[SetPrecision[Im[y0],12]];
yim=0;
If[yim0>1.1*10^-5,yim=yim0];
d=SetPrecision[If[d==0,99.,-Log[10,d]],3];
ene[z]=yre+I yim;
Write[os,SetPrecision[lambda,3]//FortranForm,OutputForm[" \t"],
SetPrecision[Delta,3]//FortranForm,OutputForm[" \t"],
yre//FortranForm,OutputForm[" \t"],yim//FortranForm,
   OutputForm[" \t"],d//FortranForm];
Print[SetPrecision[lambda,3],OutputForm[" \t"],
SetPrecision[Delta,3],OutputForm[" \t"],
yre,OutputForm[" \t"],yim,
   OutputForm[" \t"],d],
{Delta,xmin,xmax,xstep}];
Clear[p,pol,t,s0,s,c,r,rt],
{lambda,gmin,gmax,gstep}];
Close[os];