In Mathematica, the 6th degree convergence algorithm for n^(1/n), by Marvin Ray Burns
In Mathematica, the 6th degree convergence algorithm for n^(1/n), by Marvin Ray Burns
compare it with the Chen-Hash 4th degree and brute force fro a nth root of a large integer. https://www.sciencedirect.com/science/article/pii/0898122189900242
compare it with the Chen-Hash 4th degree and brute force fro a nth root of a large integer. https://www.sciencedirect.com/science/article/pii/0898122189900242
ClearSystemCache[];(*::Section::*)(*Sixth–OrderPadéRootAlgorithmforn^(1/n)*)Clear[RootPade6];RootPade6[n_Integer,prec_Integer]:=Module[{x,pc,z,t,N0=n,A1,A2,A3,B1,B2},(*Padécoefficientsfordegree[3/2]*)A1=3(2N0+1)/(5N0);A2=3(N0+1)(2N0+1)/(20N0^2);A3=(N0+1)(2N0+1)/(60N0^3);B1=2(3N0-1)/(5N0);B2=(2N0-1)(3N0-1)/(20N0^2);(*initialseed*)x=N[N0^(1/N0),100];pc=Precision[x];While[pc<prec,pc=Min[6pc,prec];(*sixtimesprecisioneachstep*)x=SetPrecision[x,pc];(*computez_i=(n-x^n)/x^n*)t=x^N0;z=(N0-t)/t;(*6th–orderiteration:r_{i+1}=r_iR_N(z)*)x=x*(1+A1z+A2z^2+A3z^3)/(1+B1z+B2z^2);];N[x,prec]];n=1234567891;(*nistheninn^(1/n)*)x=N[n^(1/n),100];(*xstartsoutasarelativelysmallprecisionapproximationton^(1/n)*)pc=Precision[x];pr=20000000;(*pristhedesiredprecisionofyourn^(1/n)*)Print[,t0=Timing[While[pc<pr,pc=Min[4pc,pr];x=SetPrecision[x,pc];y=x^n;z=(n-y)/y;t=2n-1;t2=t^2;x=x*(1+SetPrecision[4.5,pc](n-1)/t2+(n+1)z/(2nt)-SetPrecision[13.5,pc]n(n-1)/(3nt2+t^3z))];(*YougetamuchfasterversionofN[n^(1/n),pr]*)N[n-x^n,20]](*Theerror*)];t=Timing[p=RootPade6[n,pr]][[1]];Print[{"6th degree{time in sec, error}->",t,N[n-,20]}];ClearSystemCache[];n=123456789;Print["Brute force {time in sec, error}->",t1=Timing[N[n-N[n^(1/n),pr]^n,20]]]
n
p
4th degree {time in sec, error}->{50.0574,0.×}
-19999982
10
{6th degree{time in sec, error}->,28.0504,0.×}
-19999973
10
Brute force {time in sec, error}->{228.817,0.×}
-19999984
10
10,000 digits
10,000 digits
Timing for 10, 20, and 100,000 digits of the MRB constant:
Timing for 10, 20, and 100,000 digits of the MRB constant:
Print["Start time is ",ds=DateString[],"."];prec=10000;(*Numberofrequireddecimals.*)ClearSystemCache[];T0=SessionTime[];(*Sixth–orderPadékernelforn^(1/n)*)Clear[RootPade6];RootPade6[n_Integer,prec_Integer]:=Module[{x,pc,z,t,N0=n,A1,A2,A3,B1,B2},(*initialseedatmodestprecision*)x=N[N0^(1/N0),10];pc=Precision[x];While[pc<prec,pc=Min[6pc,prec];x=SetPrecision[x,pc];(*coefficientswithcurrentworkingprecisionpc*)A1=SetPrecision[3(2N0+1)/(5N0),pc];A2=SetPrecision[3(N0+1)(2N0+1)/(20N0^2),pc];A3=SetPrecision[(N0+1)(2N0+1)/(60N0^3),pc];B1=SetPrecision[2(3N0-1)/(5N0),pc];B2=SetPrecision[(2N0-1)(3N0-1)/(20N0^2),pc];(*residualz_i=(n-x^n)/x^n*)t=x^N0;z=(N0-t)/t;(*6th–orderPadéupdate*)x=x*(1+A1z+A2z^2+A3z^3)/(1+B1z+B2z^2);];N[x,prec]]expM[pre_]:=Module[{a,d,s,k,bb,c,end,xvals,x,pc,cores=16,tsize=2^7,chunksize,start=1,ll,ctab,pr=Floor[1.005pre]},chunksize=cores*tsize;n=Floor[1.32pr];end=Ceiling[n/chunksize];Print["Iterations required: ",n];Print["Will give ",end," time estimates, each more accurate than the previous."];Print["Will stop at ",end*chunksize," iterations to ensure precision of around ",pr," decimal places."];d=ChebyshevT[n,3];{b,c,s}={SetPrecision[-1,1.1*n],-d,0};Do[xvals=Flatten[ParallelTable[Table[ll=start+j*tsize+l;(*solekernelforll^(1/ll)atprecisionpr*)x=RootPade6[ll,pr];x,{l,0,tsize-1}],{j,0,cores-1},Method"EvaluationsPerKernel"32]];ctab=ParallelTable[Table[c=b-c;ll=start+l-2;b*=2(ll+n)(ll-n)/((ll+1)(2ll+1));c,{l,chunksize}],Method"EvaluationsPerKernel"16];s+=ctab.(xvals-1);start+=chunksize;st=SessionTime[]-T0;kc=k*chunksize;ti=(st)/(kc+10^-4)*(n)/(3600)/(24);If[kc>1,Print[kc," iterations done in ",N[st,4]," seconds."," Should take ",N[ti,4]," days or ",N[ti*24*3600,4],"s, finish ",DatePlus[ds,ti],"."]];,{k,0,end-1}];N[-s/d,pr]]t2=Timing[MRB=expM[prec];];Print["Finished on ",DateString[],". Processor time was ",t2[[1]]," s."];Print["error= ",N[mtest-MRB,20]];Print["Enter MRB to print ",Floor[Precision[MRB]]," digits"];
20,000 digits
20,000 digits
100,000 digits
100,000 digits