I've been working on tweaking the late Richard Crandall's code for calculating many digits of the MRB constant since 2012.

I don't expect anyone else to understand the algorithms as much as I do. (My hat's off to you if you do!!!!!! and I would be thankful for your or anyone else's help.)

I started to pay a fortune for coding help, but I thought I would give the community a chance first. It would cost $800.00 just for them to work 4 hours, probably doing little more than understanding the algorithms used.

In short, it uses my shortcut of Cohen's Algorithm 1 found here and three different algorithms to calculate Exp[Log[x]/x] i.e. x^(1/x). The first one is simply Mathematica's command. The last two are optimized for calculating more and more digits. They are faster than Newton's method for many digits! There is an undocumented use of ParallelTable (the second occurrence) which I found to speed up the code significantly!

For one thing, I was hoping someone could look at my looping (Do, the tables, While, etc) and see if you can make the program for 50,000 then 100,000+ digits run faster that way? There are two separate tables and for 50,000+ digits, my computer spends around 40-45% of the time idle!!!! I think if both can be put into one table it will spend less time idle.

For another thing, there is a second algorithm in the above Cohen link that converges with less error per iteration that I plan on trying. (I haven't worked too much on it yet because I don't understand it, but I am still looking into it.) Any help here would be great.

I will give you a great amount of credit for any help that I can use!!!!!!!!!!!!

This is code is the essence of all versions I use to break MRB constant world records!

(I only repeat parts for calculating "N[Exp[Log[ll]/ll],pr]" a number of times for a minimum computation time. Lately, I also have began changing to options for the ParallelTable commands.)

If you make an improvement in this code, I will give you a fair amount of credit in all my future computations with it, like I still do for Crandall.

For checking purposes, this program uses "m3M." More than 3 million digits of m3M, (2 million + of which have been verified), are found at http://marvinrayburns.com/3M.txt . You can import the file into Mathematica, remove the string-quotes (if they appear before and after the text) and evaluate the cell to get m3M.

The denominator for Cohen's algorithm and the needed 65536 iterations are done in 295.6 seconds on my swiftest computer. While troubleshooting, I've found it useful to reduce the numer of desired digits to 10000. Also when running these parallel tables from a fresh kernel initialization, it is good to reduce the digits to 1000, and then increase it to your desired level, because it takes 1 to a few seconds to initialize all the kernels.

If you are using the Wolfram Cloud, you might have to remove the remarks.

I suggest that in one notebook you enter the

Import["http://marvinrayburns.com/3M.txt"]

and in another notebook enter the code:

  Print["Start time is ", ds = DateString[], "."];
prec = 50000;
(**Number of required decimals.*.*)ClearSystemCache[];
T0 = SessionTime[];
expM[pre_] :=

  Block[

   {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.003 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/108];

   With[{$t = tsize - 1, $c = cores - 1, $ch = chunksize},

    Do[


     xvals = Flatten[ParallelTable[Table[ll = start + j*tsize + l;
         x = N[ll^(1/ll), iprec];
         pc = iprec;



         While[pc < pr/4, pc = Min[3 pc, pr/4];
          x = SetPrecision[x, pc];
          SetPrecision[y = x^ll - ll;
           x = x (1 - 2 y/((ll + 1) y + 2 ll ll));, pc]];(**N[Exp[Log[
         ll]/ll],pr/4]**)






         x = SetPrecision[x, pr];
         SetPrecision[xll = x^ll; z = (ll - xll)/xll;
          t = 2 ll - 1; t2 = t^2;

          x *= (1 + 45/10 (ll - 1)/t2 + (ll + 1) z/(2 ll t) - 
             135/10 ll (ll - 1) 1/(3 ll t2 + t^3 z)), pr];(**N[Exp[
         Log[ll]/ll],pr]**)







         x, {l, 0, $t}], {j, 0, $c}, 
        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, $ch}], 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["Denominator and ", 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], "."], 
      Print["Denominator for acceleration method calculated in ", st, 
       " s."]];(*I want to keep all progress stats, 
     but making them faster is fine!*)


     , {k, 0, end - 1}](*End of Do*)];


   N[-s/d, pr] 

   ];(* End of Block*)




t2 = Timing[MRBtest2 = expM[prec];]; Print["Finished on ", 
 DateString[], ". Processor time and absolute time were ", 
 t2[[1]], " and ", st, "s, respectively."];
(*Print[MRBtest2];*)Print["Enter MRBtest2 to print ", 
 Floor[Precision[
   MRBtest2]], " digits"]; Print["If you saved m3M, the difference \
between this and 3,000,000+ known digits is ", N[MRBtest2 - m3M, 10]]