Message Boards Message Boards

Convert Crandall's Mma MRB Constant code to C, C++ or Turbo C++

The late Richard Crandall worked to develop the following robust Mathematica code for calculating many digits of the MRB constant. Just before he died he did something with it to increase it's speed many fold! as seen in Try to Beat these MRB Constant Records! . After much consideration I decided he converted it to a low level language like C, C++ or Turbo C++. Unfortunately he died before he could share it with enough people. Can anyone here translate the code into one or more of those languages for the world?

It has two flavors: first plain and then multi threaded.

First:

(*Newer loop with Newton interior.*)prec = 5000;(*Number of required \
decimals.*)expM[pr_] := Module[{a, d, s, k, b, c}, n = Floor[1.32 pr];
  Print["Iterations required: ", n];
  d = N[(3 + Sqrt[8])^n, pr + 10];
  d = Round[1/2 (d + 1/d)];
  {b, c, s} = {-1, -d, 0};
  T0 = SessionTime[];
  Do[c = b - c;
   x = N[E^(Log[k + 1]/(k + 1)), iprec = Ceiling[prec/128]];
   pc = iprec;
   Do[nprec = Min[2 pc, pr];
    x = SetPrecision[x, nprec];(*Adjust precision*)
    x = N[x - x/(k + 1) + 1/x^k, nprec];
    pc *= 2;
    If[nprec >= pr, Break[]], {ct, 1, 19}];
   s += c*(x - 1);
   b *= 2 (k + n) (k - n)/((k + 1) (2 k + 1));
   If[Mod[k, 1000] == 0, 
    Print["Iterations: ", k, "    Cumulative time (sec): ", 
      SessionTime[] - T0];], {k, 0, n - 1}];
  N[-s/d, pr]];

MRBtest2 = expM[prec]

Then for the adventures,

(*Fastest (at RC's end) as of 30 Nov 2012.*)prec = 500000;(*Number of \
required decimals.*)ClearSystemCache[];
T0 = SessionTime[];
expM[pre_] := 
  Module[{a, d, s, k, bb, c, n, end, iprec, xvals, x, pc, cores = 4, 
    tsize = 2^7, chunksize, start = 1, ll, ctab, 
    pr = Floor[1.02 pre]}, chunksize = cores*tsize;
   n = Floor[1.32 pr];
   end = Ceiling[n/chunksize];
   Print["Iterations required: ", n];
   Print["end ", end];
   Print[end*chunksize];
   d = N[(3 + Sqrt[8])^n, pr + 10];
   d = Round[1/2 (d + 1/d)];
   {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, pc = Min[3 pc, pr];
         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]*)x, {l, 0, tsize - 1}], {j, 0, cores - 1}, 
       Method -> "EvaluationsPerKernel" -> 1]];
    ctab = Table[c = b - c;
      ll = start + l - 2;
      b *= 2 (ll + n) (ll - n)/((ll + 1) (2 ll + 1));
      c, {l, chunksize}];
    s += ctab.(xvals - 1);
    start += chunksize;
    Print["done iter ", k*chunksize, " ", SessionTime[] - T0];, {k, 0,
      end - 1}];
   N[-s/d, pr]];

t2 = Timing[MRBtest2 = expM[prec];];
MRBtest2 - MRBtest3

MRBtest3 is any known decimal expansion of the MRB constant. You can use the above Mathematica code to calculate as many trial digits for MRBtest3 as you wish. I have absolute confidence that the code gives accurate results for anywhere from 5000 to 2,000,000 digits.

As far as interpeting the code is concerned I noted that x = N[E^(Log[ll]/(ll)), iprec]; Gives k^(1/k) to only 1 decimal place; they are either 1.0, 1.1, 1.2, 1.3 or 1.4 (usually 1.1 or 1.0). On the other hand,

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

takes the short precision x and gives it the necessary precision and accuracy for k^(1/k) (k Is ll there.) It actually computes k^(1/k). Then he remarks, "(N[Exp[Log[ll]/ll], pr])."

After finding a fast way to compute k^(1/k) to necessary precision he uses Cohen's algorithm 1, found at Henri Cohen, Fernando Rodriguez Villegas, and Don Zagier,s Convergence Acceleration of Alternating Series, to accelerate convergence of Sum[(-1)^k*(k^(1/k) - 1), {k, 1, Infinity}].

Looking at Cohen's paper:

enter image description here enter image description here enter image description here

we see the $a_n$'s are multiplied over the c's. There Crandall craftily took advantage of Mathematica's dot product advantages.

So s += ctab.(xvals - 1);is in place of enter image description here.

POSTED BY: Marvin Ray Burns
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract