Universal CMRB code (Optimised for 10,000 -x?,000,000 digits)

In[]:=
(*================================================*)(*Sixth-order\​Padékernelforn^(1/n)*)(*================================================*)​​Clear[RootPade6];​​(*::Section::*)​​(*3rd-OrderHalleyRootAlgorithmforn^(1/n)*)​​​​Clear[RootHalley3];​​​​RootHalley3[n_Integer,prec_Integer]:=Module[{x,pc,t,f,fp,fpp,N0=n},(*initialseedatmodestprecision*)x=N[Exp[Log[N0]/N0],prec/8/6/4];​​pc=Precision[x];​​(*precisionrampwithcubic(order-3)refinement*)While[pc<prec,pc=Min[3pc,prec];(*3xprecisioneachstep*)x=SetPrecision[x,pc];​​t=x^N0;(*t=x^n*)f=t-N0;(*f(x)=x^n-n*)fp=N0t/x;(*f'(x)=nx^(n-1)=nt/x*)fpp=N0(N0-1)t/x^2;(*f''(x)=n(n-1)x^(n-2)*)(*Halleyupdate:x_{k+1}=x-(2ff')/(2f'^2-ff'')*)x=x-(2ffp)/(2fp^2-ffpp);];​​N[x,prec]];​​​​(*================================================*)​​(*Step1:precomputex_n=n^(1/n)viaRootPade6*)​​(*================================================*)​​​​xvalsPrecompute[pre_]:=​​Module[{pr,n,cores=16,tsize=2^7,chunksize,end,Ntot},​​pr=Floor[1.005pre];​​n=Floor[1.32pr];​​chunksize=cores*tsize;​​end=Ceiling[n/chunksize];​​Ntot=end*chunksize;(*totaltermsactuallyusedbyexpM*)​​ParallelTable[RootHalley3[ll,pr],{ll,1,Ntot},​​Method->"EvaluationsPerKernel"->32]]​​​​(*================================================*)​​(*Step2:Chebyshev/Eulersumoverprecomputedx*)​​(*================================================*)​​​​expMFromXvals[pre_,xvals_List]:=​​Module[{pr,n,d,b,c,s,Ntot,ll},pr=Floor[1.005pre];​​n=Floor[1.32pr];​​d=ChebyshevT[n,3];​​b=SetPrecision[-1,1.1*n];​​c=-d;​​s=0;​​Ntot=Length[xvals];​​(*globalindexLreproducingtheoriginalrecurrence*)Do[ll=L;​​c=b-c;​​s+=c*(xvals[[L+1]]-1);(*usesx_{L+1}*)​​b*=2(ll+n)(ll-n)/((ll+1)(2ll+1));,{L,0,Ntot-1}];​​N[-s/d,pr]]​​​​(*================================================*)​​(*Exampleusage(assumingmtestisdefined)*)​​(*================================================*)​​​​​​
In[]:=
AbsoluteTiming[prec=3000;​​xvals10k=xvalsPrecompute[prec]][[1]]
Out[]=
0.52963
In[]:=
AbsoluteTiming[prec=3000;​​xvals10k=xvalsPrecompute[prec]][[1]]
Out[]=
0.566462
In[]:=
Table[AbsoluteTiming[prec=10000;​​xvals10k=xvalsPrecompute[prec]][[1]],{a,1,10}]
Out[]=
{3.1669399`,3.1690454`,3.1932415`,3.1799816`,3.1598764`,3.1597313`
,
3.1669737`,3.1642072`,3.1555066`,3.1610988`}
x
In[]:=
Table[AbsoluteTiming[prec=20000;​​xvals10k=xvalsPrecompute[prec]][[1]],{a,1,10}]
Out[]=
{12.4992333`,12.4838926`,12.5221884`,12.5230455`,13.2662231`,13.5918186`,13.5466749`,13.5550278`,13.5244558`,13.5338026`}
In[]:=
Table[AbsoluteTiming[prec=30000;​​xvals10k=xvalsPrecompute[prec]][[1]],{a,1,10}]
Out[]=
{33.1996783`,34.4798806`,36.1601118`,36.1088642`,36.1456314`,36.1653597`,36.1581085`,36.1798604`,36.2166899`,36.277533`}
In[]:=
For[a=1,a<10,Print[AbsoluteTiming[prec=5000;​​xvals10k=xvalsPrecompute[prec]][[1]]];​​Print[AbsoluteTiming[MRB2=expMFromXvals[prec,xvals10k];​​N[mtest-MRB2,20]]];a++]
1.58153
{0.161944,0.×
-5021
10
}
1.49568
{0.151282,0.×
-5021
10
}
1.50504
{0.183655,0.×
-5021
10
}
1.5504
{0.201134,0.×
-5021
10
}
1.46703
{0.167608,0.×
-5021
10
}
1.52977
{0.189694,0.×
-5021
10
}
1.56019
{0.195513,0.×
-5021
10
}
1.57797
{0.198669,0.×
-5021
10
}
1.46532
{0.142556,0.×
-5021
10
}
In[]:=
For[a=1,a<10,Print[AbsoluteTiming[prec=5000;​​xvals10k=xvalsPrecompute[prec]][[1]]];​​Print[AbsoluteTiming[MRB2=expMFromXvals[prec,xvals10k];​​N[mtest-MRB2,20]]];a++]
1.12778
{0.0721416,0.×
-5021
10
}
1.10371
{0.0938062,0.×
-5021
10
}
1.08441
{0.097715,0.×
-5021
10
}
1.13601
{0.104932,0.×
-5021
10
}
1.10353
{0.132836,0.×
-5021
10
}
1.09938
{0.10498,0.×
-5021
10
}
1.11757
{0.133261,0.×
-5021
10
}
1.15132
{0.132982,0.×
-5021
10
}
1.20617
{0.107208,0.×
-5021
10
}