Wolfram Cloud Document

In[]:=

Needs["SubKernels`LocalKernels`"]Block[{$mathkernel=$mathkernel<>" -threadpriority=2"},LaunchKernels[]]

Out[]=

{KernelObject[1,local],KernelObject[2,local],KernelObject[3,local],KernelObject[4,local],KernelObject[5,local],KernelObject[6,local],KernelObject[7,local],KernelObject[8,local],KernelObject[9,local],KernelObject[10,local],KernelObject[11,local],KernelObject[12,local],KernelObject[13,local],KernelObject[14,local],KernelObject[15,local],KernelObject[16,local]}

In[]:=

ClearAll[RootPade6,expM,MRB1,MRB2,err,goodDigits];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),prec/8/6/4];pc=Precision[x];While[pc<prec,pc=Min[6pc,prec];(*6xprecisioneachstep*)x=SetPrecision[x,pc];t=x^N0;z=(N0-t)/t;(*residualforx^n=n*)(*6th-orderiteration:r_{i+1}=r_iR_N(z)*)x=x*(1+A1z+A2z^2+A3z^3)/(1+B1z+B2z^2);];N[x,prec]];pr=10000;steps=Ceiling[pr/4000];expM[pre_Integer,alpha_?NumericQ]:=Module[{d,bb,cc,s,n,pr2,block,start,stop,rng,xvals,ctab,m,dot},pr2=Floor[1.005pre];n=Floor[alphapr2];block=Ceiling[n/steps];Print["Iterations required: ",n];d=N[ChebyshevT[n,3],pr2+50];bb=SetPrecision[-1,pr2+50];cc=-d;s=SetPrecision[0,pr2+50];start=1;While[start<=n,stop=Min[start+block-1,n];Print["Starting block ",start," to ",stop];rng=Range[start,stop];(*RootPade6replacementforExp[Log[rng]/rng]-1*)xvals=ParallelTable[N[RootPade6[k,pr2]-1,pr2],{k,start,stop}];ctab=Table[cc=bb-cc;m=start+j-2;bb*=2(m+n)(m-n)/((m+1)(2m+1));cc,{j,1,stop-start+1}];dot=ctab.xvals;If[Dimensions[dot]=!={},Print["ERROR: dot is not scalar."];Print["Head[dot] = ",Head[dot]];Print["Dimensions[dot] = ",Dimensions[dot]];Abort[];];s+=dot;Print[stop," iterations done."];start=stop+1;];N[-s/d,pre]];Print["Computing MRB1..."];Print["The first run took this many seconds:",AbsoluteTiming[MRB1=expM[pr,1.32]][[1]]];Print["Computing MRB2..."];Print["the second run took this many seconds:",AbsoluteTiming[MRB2=expM[pr,1.34]][[1]]];If[ValueQ[MRB1]&&ValueQ[MRB2],err=N[Abs[MRB2-MRB1],50];Print["Error estimate = ",ScientificForm[err,20]];goodDigits=If[NumericQ[err]&&err>0,Floor[-Log10[err]],"at least "<>ToString[pr]];Print["Estimated good digits = ",goodDigits],Print["MRB1 or MRB2 failed; no valid error estimate."]];Print[pr," digits are ",MRB2];

Computing MRB1...

Iterations required: 13264

Starting block 1 to 4422

4422 iterations done.

Starting block 4423 to 8844

8844 iterations done.

Starting block 8845 to 13264

13264 iterations done.

The first run took this many seconds:3.19352

Computing MRB2...

Iterations required: 13465

Starting block 1 to 4489

4489 iterations done.

Starting block 4490 to 8978

8978 iterations done.

Starting block 8979 to 13465

13465 iterations done.

the second run took this many seconds:3.21942

Error estimate = 0.×

-10000

10

Estimated good digits = at least 10000

In[]:=

ClearAll[RootPade6,expM,MRB1,MRB2,err,goodDigits];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),prec/8/6/4];pc=Precision[x];While[pc<prec,pc=Min[6pc,prec];(*6xprecisioneachstep*)x=SetPrecision[x,pc];t=x^N0;z=(N0-t)/t;(*residualforx^n=n*)(*6th-orderiteration:r_{i+1}=r_iR_N(z)*)x=x*(1+A1z+A2z^2+A3z^3)/(1+B1z+B2z^2);];N[x,prec]];pr=20000;steps=Ceiling[pr/4000];expM[pre_Integer,alpha_?NumericQ]:=Module[{d,bb,cc,s,n,pr2,block,start,stop,rng,xvals,ctab,m,dot},pr2=Floor[1.005pre];n=Floor[alphapr2];block=Ceiling[n/steps];Print["Iterations required: ",n];d=N[ChebyshevT[n,3],pr2+50];bb=SetPrecision[-1,pr2+50];cc=-d;s=SetPrecision[0,pr2+50];start=1;While[start<=n,stop=Min[start+block-1,n];Print["Starting block ",start," to ",stop];rng=Range[start,stop];(*RootPade6replacementforExp[Log[rng]/rng]-1*)xvals=ParallelTable[N[RootPade6[k,pr2]-1,pr2],{k,start,stop}];ctab=Table[cc=bb-cc;m=start+j-2;bb*=2(m+n)(m-n)/((m+1)(2m+1));cc,{j,1,stop-start+1}];dot=ctab.xvals;If[Dimensions[dot]=!={},Print["ERROR: dot is not scalar."];Print["Head[dot] = ",Head[dot]];Print["Dimensions[dot] = ",Dimensions[dot]];Abort[];];s+=dot;Print[stop," iterations done."];start=stop+1;];N[-s/d,pre]];Print["Computing MRB1..."];Print["The first run took this many seconds:",AbsoluteTiming[MRB1=expM[pr,1.32]][[1]]];Print["Computing MRB2..."];Print["the second run took this many seconds:",AbsoluteTiming[MRB2=expM[pr,1.34]][[1]]];If[ValueQ[MRB1]&&ValueQ[MRB2],err=N[Abs[MRB2-MRB1],50];Print["Error estimate = ",ScientificForm[err,20]];goodDigits=If[NumericQ[err]&&err>0,Floor[-Log10[err]],"at least "<>ToString[pr]];Print["Estimated good digits = ",goodDigits],Print["MRB1 or MRB2 failed; no valid error estimate."]];Print[pr," digits are ",MRB2];