Geometrically Constructing the MRB Constant

With the Wolfram Assistant

as seen in this notebook:

While preparing for the verification of some future 7 million digits, I broke some speed records with two of my i9-14900K 6400MHz RAM (overclocked CPUs), using Mathematica 11.3 and the lightweight grid. They are generally faster than the 3-node MRB constant supercomputer with remote kernels! Unless noted by the command Timing[], i.e., the last yellow and the final two blue-highlighted columns in the fourth table, these are all absolute timings. How does your computer compare to these? What can you do with other software?

A few of the record runs are missing from the documentation below because of repeating the run, and many of the runs listed are not in the records.

The 10 M (415-day, 64GB RAM) digit time was estimated from this notebook.

The 10 M (205-day, over 128GB RAM) digit time was estimated from this notebook.

The 10M (350-day, under 128GB RAM) digit time is estimated and displayed in this older version here and in this updated one here.

The 7 M (7 million) digit time was estimated from this notebook.

The 100k and 300k parts of the documentation for the 1X then 2X 14900KS Parallel Runs is here. They were run sequentially, only "symbolically" parallel, since I don't have 10 14900KSs.

Enough of the 1 M to estimate wall time is here.

The sixth-order convergence algorithm is my own work, from Mathematica:

  PadeApproximant[x^(2/x), {x, 1, {3, 2}}]

(* (1 + 67/20 (-1 + x) + 161/60 (-1 + x)^2 - 
 23/60 (-1 + x)^3)/(1 + 27/20 (-1 + x) + 59/60 (-1 + x)^2)*)

Here are the code and the results:

(* ================================================*)(*Sixth-order \
Padé kernel for n^(1/n)*)(* ================================================*)


Clear[RootPade6];
RootPade6[n_Integer, prec_Integer] := 
  Module[{x, pc, z, t, N0 = n, A1, A2, A3, B1, 
      B2},(*initial seed at modest precision*)
  x = N[N0^(1/N0), prec/1000];
    pc = Precision[x];
    While[pc < pr/180, pc = Min[3 pc, pr/180];
      x = SetPrecision[x, pc];
      y = x^ll - ll;
      x = x (1 - 2 y/((ll + 1) y + 2 ll ll));];
    (**N[Exp[Log[ll]/ll],pr/60]**)
  While[pc < prec, pc = Min[6 pc, prec];
      x = SetPrecision[x, pc];
      (*coefficients with current working precision pc*)
      A1 = SetPrecision[3 (2 N0 + 1)/(5 N0), pc];
      A2 = SetPrecision[3 (N0 + 1) (2 N0 + 1)/(20 N0^2), pc];
      A3 = SetPrecision[(N0 + 1) (2 N0 + 1)/(60 N0^3), pc];
      B1 = SetPrecision[2 (3 N0 - 1)/(5 N0), pc];
      B2 = SetPrecision[(2 N0 - 1) (3 N0 - 1)/(20 N0^2), pc];
      (*residual z_i=(n-x^n)/x^n*)t = x^N0;
      z = (N0 - t)/t;
      (*6th-order Padé update*)
      x = x*(1 + A1 z + A2 z^2 + A3 z^3)/(1 + B1 z + B2 z^2);];
    N[x, prec]]

(* ================================================*)
(*Step 1:precompute x_n=n^(1/n) via RootPade6*)
(* ================================================*)

xvalsPrecompute[pre_] := 
  Module[{pr, n, cores = 1, tsize = 2^7, chunksize, end, Ntot}, 
    pr = Floor[1.005 pre];
    n = Floor[1.32 pr];
    chunksize = cores*tsize;
    end = Ceiling[n/chunksize];
    Ntot = end*chunksize;(*total terms actually used by expM*)
    ParallelTable[RootPade6[ll, pr], {ll, 1, Ntot}, 
      Method -> "EvaluationsPerKernel" -> 32]]

(* ================================================*)
(*Step 2:Chebyshev/Euler sum over precomputed x*)
(* ================================================*)

expMFromXvals[pre_, xvals_List] := 
  Module[{pr, n, d, b, c, s, Ntot, ll}, pr = Floor[1.005 pre];
    n = Floor[1.32 pr];
    d = ChebyshevT[n, 3];
    b = SetPrecision[-1, 1.1*n];
    c = -d;
    s = 0;
    Ntot = Length[xvals];
    (*global index L reproducing the original recurrence*)Do[ll = L;
      c = b - c;
      s += c*(xvals[[L + 1]] - 1);(*uses x_{L+1}*)
      b *= 2 (ll + n) (ll - n)/((ll + 1) (2 ll + 1));, {L, 0, 
    Ntot - 1}];
    N[-s/d, pr]]

(* ================================================*)
(*Example usage (assuming mtest is defined)*)
(* ================================================*)




Print[AbsoluteTiming[prec = 10000;
      xvals10k = xvalsPrecompute[prec];]];
Print[AbsoluteTiming[MRB2 = expMFromXvals[prec, xvals10k];
    N[mtest - MRB2, 20]]]

With 1X I-14900KS: With 2X I-14900KS: With 3X I-14900KS:

A few of the record runs are missing from the documentation below, and many of the documents listed are not in the records.**

That light blue highlight, red text columns, which shows processor times result here->7200 corrected timing for 30k and was tweaked in 2025.

For the partial column of two 6000MHz 14900K' with red text and yellow highlight, see speed 100 300 1M XMP tweaked.

For column "=F" (highlighted in green) see linked "10203050100" .

At the bottom, see attached "kernel priority 2 computers.nb" for column =B,

"3 fastest computers together.nb" for column =C

and linked "speed records 5 10 20 30 K"

Most of the sixth-degree results are here

For the mostly red column, including the single record, 10,114 seconds, 300,000 digit run " =E" is in the linked "3 fastest computers together 2.nb.}

For column "=J," see 574 s 100k , .106.1 sec 50k and 6897s 300k

The 27-hour million-digit computation is found here.

Print["Start time is ",ds=DateString[],"."];
prec=1000000;
(**Number of required decimals.*.*)ClearSystemCache[];
T0=SessionTime[];
expM[pre_]:=Module[{a,d,s,k,bb,c,end,iprec,xvals,x,pc,cores=16(*=4*number of physical cores*),tsize=2^7,chunksize,start=1,ll,ctab,pr=Floor[1.005 pre]},chunksize=cores*tsize;\[IndentingNewLine]n=Floor[1.32 pr];\[IndentingNewLine]end=Ceiling[n/chunksize];\[IndentingNewLine]Print["Iterations required: ",n];\[IndentingNewLine]Print["Will give ",end," time estimates, each more accurate than the previous."];\[IndentingNewLine]Print["Will stop at ",end*chunksize," iterations to ensure precsion of around ",pr," decimal places."];d=ChebyshevT[n,3];\[IndentingNewLine]{b,c,s}={SetPrecision[-1,1.1*n],-d,0};\[IndentingNewLine]iprec=Ceiling[pr/186624];\[IndentingNewLine]Do[xvals=Flatten[Parallelize[Table[Table[ll=start+j*tsize+l;\[IndentingNewLine]x=N[E^(Log[ll]/(ll)),iprec];\[IndentingNewLine]pc=iprec;\[IndentingNewLine]While[pc<pr/1024,pc=Min[3 pc,pr/1024];\[IndentingNewLine]x=SetPrecision[x,pc];\[IndentingNewLine]y=x^ll-ll;\[IndentingNewLine]x=x (1-2 y/((ll+1) y+2 ll ll));];\[IndentingNewLine](**N[Exp[Log[ll]/ll],pr/1024]**)x=SetPrecision[x,pr/256];\[IndentingNewLine]xll=x^ll;z=(ll-xll)/xll;\[IndentingNewLine]t=2 ll-1;t2=t^2;\[IndentingNewLine]x=x*(1+SetPrecision[4.5,pr/256] (ll-1)/t2+(ll+1) z/(2 ll t)-SetPrecision[13.5,pr/256] ll (ll-1) 1/(3 ll t2+t^3 z));(*N[Exp[Log[ll]/ll],pr/256]*)x=SetPrecision[x,pr/64];\[IndentingNewLine]xll=x^ll;z=(ll-xll)/xll;\[IndentingNewLine]t=2 ll-1;t2=t^2;\[IndentingNewLine]x=x*(1+SetPrecision[4.5,pr/64] (ll-1)/t2+(ll+1) z/(2 ll t)-SetPrecision[13.5,pr/64] ll (ll-1) 1/(3 ll t2+t^3 z));(**N[Exp[Log[ll]/ll],pr/64]**)x=SetPrecision[x,pr/16];\[IndentingNewLine]xll=x^ll;z=(ll-xll)/xll;\[IndentingNewLine]t=2 ll-1;t2=t^2;\[IndentingNewLine]x=x*(1+SetPrecision[4.5,pr/16] (ll-1)/t2+(ll+1) z/(2 ll t)-SetPrecision[13.5,pr/16] ll (ll-1) 1/(3 ll t2+t^3 z));(**N[Exp[Log[ll]/ll],pr/16]**)x=SetPrecision[x,pr/4];\[IndentingNewLine]xll=x^ll;z=(ll-xll)/xll;\[IndentingNewLine]t=2 ll-1;t2=t^2;\[IndentingNewLine]x=x*(1+SetPrecision[4.5,pr/4] (ll-1)/t2+(ll+1) z/(2 ll t)-SetPrecision[13.5,pr/4] ll (ll-1) 1/(3 ll t2+t^3 z));(**N[Exp[Log[ll]/ll],pr/4]**)x=SetPrecision[x,pr];\[IndentingNewLine]xll=x^ll;z=(ll-xll)/xll;\[IndentingNewLine]t=2 ll-1;t2=t^2;\[IndentingNewLine]x=x*(1+SetPrecision[4.5,pr] (ll-1)/t2+(ll+1) z/(2 ll t)-SetPrecision[13.5,pr] ll (ll-1) 1/(3 ll t2+t^3 z));(*N[Exp[Log[ll]/ll],pr]*)x,{l,0,tsize-1}],{j,0,cores-1}]]];\[IndentingNewLine]ctab=ParallelTable[Table[c=b-c;\[IndentingNewLine]ll=start+l-2;\[IndentingNewLine]b*=2 (ll+n) (ll-n)/((ll+1) (2 ll+1));\[IndentingNewLine]c,{l,chunksize}],Method->"Automatic"];\[IndentingNewLine]s+=ctab.(xvals-1);\[IndentingNewLine]start+=chunksize;\[IndentingNewLine]st=SessionTime[]-T0;kc=k*chunksize;\[IndentingNewLine]ti=(st)/(kc+10^-4)*(n)/(3600)/(24);\[IndentingNewLine]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}];\[IndentingNewLine]N[-s/d,pr]];
t2=Timing[MRBtest2=expM[prec];];Print["Finished on ",DateString[],". Proccessor time was ",t2[[1]]," s."];Print["Actual time was ",st];
(*Print[*)MRBtest2(*]*)(*Remove (**) or enter MRBtest2 to print output*);Print["Enter MRBtest2 to print ",Floor[Precision[MRBtest2]]," digits"];Print["If you saved m3M, the difference between this and 3,014,991 known digits is ",N[MRBtest2-m3M,10]]

Remembering that the integrated analog of the MRB constant is

NIntegrate[(-1)^n (n^(1/n) - 1), {n, 1, Infinity  I}, 
 Method -> "Trapezoidal", WorkingPrecision -> 20]

These results are from the Timing[] command:

The 14900KS at 7200 MHz (extreme tuning!) documented here

The i9-12900KS column is documented here.

"Windows10 2024 i9-14900KS 6000 MHZ RAM" documentation here

Attachments:

3 fastest comput...nb

kernel priority ...nb

Distribution of MRB constant digits as given in all record calculations

Here is the distribution of digits within the first 6,000,000 decimal places (.187859,,,), "4" shows up more than other digits, followed by "0," "8" and "7."

Here is the distribution of digits within the first 5,000,000 decimal places (.187859,,,), "4" shows up a lot more than other digits, followed by "0," "8" and "6."

Here is a similar distribution over the first 4,000,000:

3,000,000 digits share a similar distribution:

Over the first 2 and 1 million digits "4" was not so well represented. So, the heavy representation of "4" is shown to be a growing phenomenon from 2 million to 5 million digits. However, "1,2,and 5" still made a very poor showing:

I attached more than 6,000,000 digits of the MRB constant.

I calculated 6,500,000 digits of the MRB constant!!

The MRB constant supercomputer said,

Finished on Wed 16 Mar 2022 02 : 02 : 10. Processor and actual time were 6.2662810^6 and 1.6026403541959210^7 s.respectively Enter MRB1 to print 6532491 digits. The error from a 6, 000, 000 or more digit calculation that used a different method is 0.*10^-6029992

"Processor time" 72.526 days

"Actual time" 185.491 days

For the digits see the attached 6p5millionMRB.nb. For the documentation of the computation see 2nd 6p5 million.nb.

Attachments:

2nd 6p5 million.nb

6p5millionMRB.nb

4/22/2019

Let $$M=\sum _{n=1}^{\infty } \frac{(-1)^{n+1} \eta ^n(n)}{n!}=\sum _{n=1}^{\infty } (-1)^n \left(n^{1/n}-1\right).$$ Then using what I learned about the absolute convergence of $\sum _{n=1}^{\infty } \frac{(-1)^{n+1} \eta ^n(n)}{n!}$ from https://math.stackexchange.com/questions/1673886/is-there-a-more-rigorous-way-to-show-these-two-sums-are-exactly-equal, combined with an identity from Richard Crandall: , Also using what Mathematica says:

$$\sum _{n=1}^1 \frac{\underset{m\to 1}{\text{lim}} \eta ^n(m)}{n!}=\gamma (2 \log )-\frac{2 \log ^2}{2},$$

I figured out that

$$\sum _{n=2}^{\infty } \frac{(-1)^{n+1} \eta ^n(n)}{n!}=\sum _{n=1}^{\infty } (-1)^n \left(n^{1/n}-\frac{\log (n)}{n}-1\right).$$

So I made the following major breakthrough in computing MRB from Candall's first eta formula. See attached 100 k eta 4 22 2019. Also shown below.

The time grows 10,000 times slower than the previous method!

I broke a new record, 100,000 digits: Processor and total time were 806.5 and 2606.7281972 s respectively.. See attached 2nd 100 k eta 4 22 2019.

Here is the work from 100,000 digits.

Print["Start time is ", ds = DateString[], "."];
prec = 100000;
(**Number of required decimals.*.*)ClearSystemCache[];
T0 = SessionTime[];
expM[pre_] := 
  Module[{a, d, s, k, bb, c, end, iprec, xvals, x, pc, cores = 16(*=4*
    number of physical cores*), tsize = 2^7, chunksize, start = 1, ll,
     ctab, pr = Floor[1.005 pre]}, chunksize = cores*tsize;
   n = Floor[1.32 pr];
   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 precsion of around ", pr, 
    " decimal places."]; d = ChebyshevT[n, 3];
   {b, c, s} = {SetPrecision[-1, 1.1*n], -d, 0};
   iprec = Ceiling[pr/27];
   Do[xvals = Flatten[ParallelTable[Table[ll = start + j*tsize + l;
        x = N[E^(Log[ll]/(ll)), iprec];
        pc = iprec;
        While[pc < pr/4, pc = Min[3 pc, pr/4];
         x = SetPrecision[x, pc];
         y = x^ll - ll;
         x = x (1 - 2 y/((ll + 1) y + 2 ll ll));];(**N[Exp[Log[ll]/
        ll],pr/4]**)x = SetPrecision[x, pr];
        xll = x^ll; z = (ll - xll)/xll;
        t = 2 ll - 1; t2 = t^2;
        x = 
         x*(1 + SetPrecision[4.5, pr] (ll - 1)/
              t2 + (ll + 1) z/(2 ll t) - 
            SetPrecision[13.5, pr] ll (ll - 1) 1/(3 ll t2 + t^3 z));(**
        N[Exp[Log[ll]/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) (2 ll + 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[], ". Proccessor time was ", t2[[1]], " s."];
Print["Enter MRBtest2 to print ", Floor[Precision[MRBtest2]], 
  " digits"];


 (Start time is )^2Tue 23 Apr 2019 06:49:31.

 Iterations required: 132026

 Will give 65 time estimates, each more accurate than the previous.

 Will stop at 133120 iterations to ensure precsion of around 100020 decimal places.

 Denominator computed in  17.2324041s.

...

129024 iterations done in 1011. seconds. Should take 0.01203 days or 1040.s, finish Mon 22 Apr 
2019 12:59:16.

131072 iterations done in 1026. seconds. Should take 0.01202 days or 1038.s, finish Mon 22 Apr 
2019 12:59:15.

Finished on Mon 22 Apr 2019 12:59:03. Processor time was 786.797 s.

 Print["Start time is " "Start time is ", ds = DateString[], "."];
 prec = 100000;
 (**Number of required decimals.*.*)ClearSystemCache[];
 T0 = SessionTime[];
 expM[pre_] := 
   Module[{lg, a, d, s, k, bb, c, end, iprec, xvals, x, pc, cores = 16(*=
     4*number of physical cores*), tsize = 2^7, chunksize, start = 1, 
     ll, ctab, pr = Floor[1.0002 pre]}, chunksize = cores*tsize;
    n = Floor[1.32 pr];
    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 precsion of around ", pr, 
     " decimal places."]; d = ChebyshevT[n, 3];
    {b, c, s} = {SetPrecision[-1, 1.1*n], -d, 0};
    iprec = pr/2^6;
    Do[xvals = Flatten[ParallelTable[Table[ll = start + j*tsize + l;
         lg = Log[ll]/(ll); x = N[E^(lg), iprec];
         pc = iprec;
         While[pc < pr, pc = Min[4 pc, pr];
          x = SetPrecision[x, pc];
          xll = x^ll; z = (ll - xll)/xll;
          t = 2 ll - 1; t2 = t^2;
          x = 
           x*(1 + SetPrecision[4.5, pc] (ll - 1)/
                t2 + (ll + 1) z/(2 ll t) - 
              SetPrecision[13.5, 2 pc] ll (ll - 1)/(3 ll t2 + t^3 z))];
          x - lg, {l, 0, tsize - 1}], {j, 0, cores - 1}, 
        Method -> "EvaluationsPerKernel" -> 16]];
     ctab = ParallelTable[Table[c = b - c;
        ll = start + l - 2;
        b *= 2 (ll + n) (ll - n)/((ll + 1) (2 ll + 1));
        c, {l, chunksize}], Method -> "EvaluationsPerKernel" -> 16];
     s += ctab.(xvals - 1);
     start += chunksize;
     st = SessionTime[] - T0; kc = k*chunksize;
     ti = (st)/(kc + 10^-10)*(n)/(3600)/(24);
     If[kc > 1, 
      Print[kc, " iterations done in ", N[st - stt, 4], " seconds.", 
       " Should take ", N[ti, 4], " days or ", ti*3600*24, 
       "s, finish ", DatePlus[ds, ti], "."], 
      Print["Denominator computed in  ", stt = st, "s."]];, {k, 0, 
      end - 1}];
    N[-s/d, pr]];
 t2 = Timing[MRBeta2toinf = expM[prec];]; Print["Finished on ", 
  DateString[], ". Processor and total time were ", 
  t2[[1]], " and ", st, " s respectively."];

Start time is  Tue 23 Apr 2019 06:49:31.

Iterations required: 132026

Will give 65 time estimates, each more accurate than the previous.

Will stop at 133120 iterations to ensure precision of around 100020 decimal places.

Denominator computed in  17.2324041s.

...

131072 iterations done in 2589. seconds. Should take 0.03039 days or 2625.7011182s, finish Tue 23 Apr 2019 07:33:16.

Finished on Tue 23 Apr 2019 07:32:58. Processor and total time were 806.5 and 2606.7281972 s respectively.

 MRBeta1 = EulerGamma Log[2] - 1/2 Log[2]^2

 EulerGamma Log[2] - Log[2]^2/2

   N[MRBeta2toinf + MRBeta1 - MRB, 10]

   1.307089967*10^-99742

Attachments:

2nd 100 k eta 4 ...nb

50 k eta 4 22 2019.nb

The new sum is this.

Sum[(-1)^(k + 1)*(-1 + (1 + k)^(1/(1 + k)) - Log[1 + k]/(1 + k) - 
         Log[1 + k]^2/(2*(1 + k)^2)), {k, 0, Infinity}] 

That appears to be the same as for MRB except now we subtract two terms from the series expansion at the origin of k^(1/k). For each k these terms are Log[k]/k + 1/2*(Log[k]/k)^2. Accounting for the signs (-1)^k and summing, as I did earlier for just that first term, we get something recognizable.

Sum[(-1)^(k)*(Log[k]/(k) + Log[k]^2/(2*k^2)), {k, 1, Infinity}]

(* Out[21]= 1/24 (24 EulerGamma Log[2] - 2 EulerGamma \[Pi]^2 Log[2] - 
   12 Log[2]^2 - \[Pi]^2 Log[2]^2 + 24 \[Pi]^2 Log[2] Log[Glaisher] - 
   2 \[Pi]^2 Log[2] Log[\[Pi]] - 6 (Zeta^\[Prime]\[Prime])[2]) *)

So what does this buy us? For one thing, we get even better convergence from brute force summation, because now our largest terms are O((logk/k)^3) and alternating (which means if we sum in pairs it's actually O~(1/k^4) with O~ denoting the "soft-oh" wherein one drops polylogarithmic factors).

How helpful is this? Certainly it cannot hurt. But even with 1/k^4 size terms, it takes a long time to get even 40 digits, let alone thousands. So there is more going on in that Crandall approach.

The identity in question is straightforward. Write n^(1/n) as Exp[Log[n]/n], take a series expansion at 0, and subtract the first term from all summands. That means subtracting off Log[n]/n in each summand. This gives your left hand side. We know it must be M - the sum of the terms we subtracted off. Now add all of them up, accounting for signs.

Expand[Sum[(-1)^n*Log[n]/n, {n, 1, Infinity}]]

(* Out[74]= EulerGamma Log[2] - Log[2]^2/2 *)

So we recover the right hand side.

I have not understood whether this identity helps with Crandall's iteration. One advantage it confers, a good one in general, is that it converts a conditionally convergent alternating series into one that is absolutely convergent. From a numerical computation point of view this is always good.

Nice work. Worth a bit of excitement, I' d say.

I can't say I understand either. My guess is the Eta stuff comes from summing (-1)^k*(Log[k]/k)^n over k, as those are the terms that appear in the double sum you get from expanding k^(1/k)-1 in powers of Log[k]/k (use k^(1/k)=Exp[Log[k]/k] and the power series for Exp). Even if it does come from this the details remain elusive..

It is hard to be certain that c1 and c2 are correct to 77 digits even though they agree to that extent. I'm not saying that they are incorrect and presumably you have verified this. Just claiming that whatever methods NSum may be using to accelerate convergence, there is really no guarantee that they apply to this particular computation. So c1 aand c2 could agree to that many places because they are computed in a similar manner without all digits actually being correct.