Message Boards Message Boards

Try to beat these MRB constant records!

GROUPS:

POSTED BY: Marvin Ray Burns .

I think this important point got buried near the end.

When it comes to mine and a few more educated people's passion to calculate many digits and the dislike possessed by a few more educated people; it is all a matter telling us that the human mind is multi faceted in giving passion, to person a, for one task and to person b for another task!

The MRB constant is defined below. See http://mathworld.wolfram.com/MRBConstant.html

enter image description here enter image description here

Here are some record computations. If you know of any others let me know..

  1. On or about Dec 31, 1998 I computed 1 digit of the (additive inverse of the) MRB constant with my TI-92's, by adding 1-sqrt(2)+3^(1/3)-4^(1/4) as far as I could. That first digit by the way is just 0.

  2. On Jan 11, 1999 I computed 3 digits of the MRB constant with the Inverse Symbolic Calculator.

  3. In Jan of 1999 I computed 4 correct digits of the MRB constant using Mathcad 3.1 on a 50 MHz 80486 IBM 486 personal computer operating on Windows 95.

  4. Shortly afterwards I computed 9 correct digits of the MRB constant using Mathcad 7 professional on the Pentium II mentioned below.

  5. On Jan 23, 1999 I computed 500 digits of the MRB constant with the online tool called Sigma.

  6. In September of 1999, I computed the first 5,000 digits of the MRB Constant on a 350 MHz Pentium II with 64 Mb of ram using the simple PARI commands \p 5000;sumalt(n=1,((-1)^n*(n^(1/n)-1))), after allocating enough memory.

  7. On June 10-11, 2003 over a period, of 10 hours, on a 450mh P3 with an available 512mb RAM: I computed 6,995 accurate digits of the MRB constant.

  8. Using a Sony Vaio P4 2.66 GHz laptop computer with 960 MB of available RAM, on 2:04 PM 3/25/2004, I finished computing 8000 digits of the MRB constant.

  9. On March 01, 2006 with a 3GH PD with 2GBRAM available, I computed the first 11,000 digits of the MRB Constant.

  10. On Nov 24, 2006 I computed 40, 000 digits of the MRB Constant in 33hours and 26min via my own program in written in Mathematica 5.2. The computation was run on a 32-bit Windows 3GH PD desktop computer using 3.25 GB of Ram.

  11. Finishing on July 29, 2007 at 11:57 PM EST, I computed 60,000 digits of MRB Constant. Computed in 50.51 hours on a 2.6 GH AMD Athlon with 64 bit Windows XP. Max memory used was 4.0 GB of RAM.

  12. Finishing on Aug 3 , 2007 at 12:40 AM EST, I computed 65,000 digits of MRB Constant. Computed in only 50.50 hours on a 2.66GH Core2Duo using 64 bit Windows XP. Max memory used was 5.0 GB of RAM.

  13. Finishing on Aug 12, 2007 at 8:00 PM EST, I computed 100,000 digits of MRB Constant. They were computed in 170 hours on a 2.66GH Core2Duo using 64 bit Windows XP. Max memory used was 11.3 GB of RAM. Median (typical) daily record of memory used was 8.5 GB of RAM.

  14. Finishing on Sep 23, 2007 at 11:00 AM EST, I computed 150,000 digits of MRB Constant. They were computed in 330 hours on a 2.66GH Core2Duo using 64 bit Windows XP. Max memory used was 22 GB of RAM. Median (typical) daily record of memory used was 17 GB of RAM.

  15. Finishing on March 16, 2008 at 3:00 PM EST, I computed 200,000 digits of MRB Constant using Mathematica 5.2. They were computed in 845 hours on a 2.66GH Core2Duo using 64 bit Windows XP. Max memory used was 47 GB of RAM. Median (typical) daily record of memory used was 28 GB of RAM.

  16. Washed away by Hurricane Ike -- on September 13, 2008 sometime between 2:00PM - 8:00PM EST an almost complete computation of 300,000 digits of the MRB Constant was destroyed. Computed for a long 4015. Hours (23.899 weeks or 1.4454*10^7 seconds) on a 2.66GH Core2Duo using 64 bit Windows XP. Max memory used was 91 GB of RAM. The Mathematica 6.0 code used follows:

    Block[{$MaxExtraPrecision = 300000 + 8, a, b = -1, c = -1 - d, 
     d = (3 + Sqrt[8])^n, n = 131 Ceiling[300000/100], s = 0}, a[0] = 1; 
     d = (d + 1/d)/2; For[m = 1, m < n, a[m] = (1 + m)^(1/(1 + m)); m++]; 
     For[k = 0, k < n, c = b - c; 
      b = b (k + n) (k - n)/((k + 1/2) (k + 1)); s = s + c*a[k]; k++]; 
     N[1/2 - s/d, 300000]]
    
  17. On September 18, 2008 a computation of 225,000 digits of MRB Constant was started with a 2.66GH Core2Duo using 64 bit Windows XP. It was completed in 1072 hours. Memory usage is recorded in the attachment pt 225000.xls, near the bottom of this post. .

  18. 250,000 digits was attempted but failed to be completed to a serious internal error which restarted the machine. The error occurred sometime on December 24, 2008 between 9:00 AM and 9:00 PM. The computation began on November 16, 2008 at 10:03 PM EST. Like the 300,000 digit computation this one was almost complete when it failed. The Max memory used was 60.5 GB.

  19. On Jan 29, 2009, 1:26:19 pm (UTC-0500) EST, I finished computing 250,000 digits of the MRB constant. with a multiple step Mathematica command running on a dedicated 64bit XP using 4Gb DDR2 Ram on board and 36 GB virtual. The computation took only 333.102 hours. The digits are at http://marvinrayburns.com/250KMRB.txt . The computation is completely documented in the attached 250000.pd at bottom of this post.

  20. On Sun 28 Mar 2010 21:44:50 (UTC-0500) EST, I started a computation of 300000 digits of the MRB constant using an i7 with 8.0 GB of DDR3 Ram on board.; But it failed due to hardware problems.

  21. I computed 299,998 Digits of the MRB constant. The computation began Fri 13 Aug 2010 10:16:20 pm EDT and ended 2.23199*10^6 seconds later | Wednesday, September 8, 2010. I used Mathematica 6.0 for Microsoft Windows (64-bit) (June 19, 2007) That is an average of 7.44 seconds per digit.. I used my Dell Studio XPS 8100 i7 860 @ 2.80 GH 2.80 GH with 8GB physical DDR3 RAM. Windows 7 reserved an additional 48.929 GB virtual Ram.

  22. I computed exactly 300,000 digits to the right of the decimal point of the MRB constant from Sat 8 Oct 2011 23:50:40 to Sat 5 Nov 2011 19:53:42 (2.405*10^6 seconds later). This run was 0.5766 seconds per digit slower than the 299,998 digit computation even though it used 16GB physical DDR3 RAM on the same machine. The working precision and accuracy goal combination were maximized for exactly 300,000 digits, and the result was automatically saved as a file instead of just being displayed on the front end. Windows reserved a total of 63 GB of working memory of which at 52 GB were recorded being used. The 300,000 digits came from the Mathematica 7.0 command

    Quit; DateString[]
    digits = 300000; str = OpenWrite[]; SetOptions[str, 
    PageWidth -> 1000]; time = SessionTime[]; Write[str, 
    NSum[(-1)^n*(n^(1/n) - 1), {n, \[Infinity]}, 
    WorkingPrecision -> digits + 3, AccuracyGoal -> digits, 
    Method -> "AlternatingSigns"]]; timeused = 
    SessionTime[] - time; here = Close[str]
    DateString[]
    
  23. 314159 digits of the constant took 3 tries do to hardware failure. Finishing on September 18, 2012 I computed 314159 digits, taking 59 GB of RAM. The digits are came from the Mathematica 8.0.4 code

    DateString[]
    NSum[(-1)^n*(n^(1/n) - 1), {n, \[Infinity]}, 
    WorkingPrecision -> 314169, Method -> "AlternatingSigns"] // Timing
    DateString[]
    

Where I have 10 digits to round off. (The command NSum[(-1)^n*(n^(1/n) - 1), {n, [Infinity]}, WorkingPrecision -> big number, Method -> "AlternatingSigns"] tends to give about 3 digits of error to the right.)

The following records are due to the work of Richard Crandall found here.

  1. Sam Noble of Apple computed 1,000,000 digits of the MRB constant in 18 days 9 hours 11 minutes 34.253417 seconds

  2. Finishing on Dec 11, 2012 Ricard Crandall, an Apple scientist, computed 1,048,576 digits in a lighting fast 76.4 hours. That's on a 2.93 Ghz 8-core Nehalem

  3. I computed a little over 1,200,000 digits of the MRB constant in 11 days, 21 hours, 17 minutes, and 41 seconds,( finishing on on March 31 2013). I used a six core Intel(R) Core(TM) i7-3930K CPU @ 3.20 GHz 3.20 GHz.

  4. On May 17, 2013 I finished a 2,000,000 or more digit computation of the MRB constant, using only around 10GB of RAM. It took 37 days 5 hours 6 minutes 47.1870579 seconds. I used a six core Intel(R) Core(TM) i7-3930K CPU @ 3.20 GHz 3.20 GHz.

  5. Finally, I would like to announce a new unofficial world record computation of the MRB constant that was finished on Sun 21 Sep 2014 18:35:06. It took 1 month 27 days 2 hours 45 minutes 15 seconds. I computed 3,014,991 digits of the MRB constant with Mathematica 10.0. I Used my new version of Richard Crandall's code, below, optimized for my platform and large computations. I also used a six core Intel(R) Core(TM) i7-3930K CPU @ 3.20 GHz 3.20 GHz with 64 GB of RAM of which only 16 GB was used. Can you beat it (in more number of digits, less memory used, or less time taken)? This confirms that my previous "2,000,000 or more digit computation" was actually accurate to 2,009,993 digits. (They were used as MRBtest2M.)

        (**Fastest (at MRB's end) as of 25 Jul 2014*.*) 
       DateString[]
       prec = 3000000; 
       (**Number of required decimals.*.*)ClearSystemCache[]; 
       T0 = SessionTime[];
         expM[pre_] :=
        Module[{a, d, s, k, bb, c, n, end, iprec, xvals, x, pc, cores = 12, 
        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["end ", end];
      Print[end*chunksize]; 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, 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" -> 4]];
       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" -> 2];
      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];]; DateString[] 
     Print[MRBtest2]  
    MRBtest2 - MRBtest2M
    

t2 From the computation was {1.961004112059*10^6, Null}.

Here are a couple of graphs of my record computations in max digits/ year:

enter image description hereenter image description here

POSTED BY: Marvin Ray Burns
Answer
3 years ago

The following might help anyone serious about breaking my record.

Richard Crandall wrote about couple of new formulae for computing the MRB constant faster on pp 28 and 29 here. a true jewel of an expression Interestingly, as pointed out by J. Borwein,

By the way, I wrote Crandall about formula (44) -- it seems it should have a minus sign in front of it -- and he wrote back: enter image description here

POSTED BY: Marvin Ray Burns
Answer
3 years ago

The following email Crandall sent me before he died might be useful for anyone checking their results: enter image description here

POSTED BY: Marvin Ray Burns
Answer
3 years ago

Perhaps some of these speed records will be easier to beat

digits      seconds                  sec./digit
500          0.062                     0.000124

1000         0.1230                  0.000123

3000         0.8640                  0.000288

5000        2.4801                    0.00049602

10,000      9.8515                   0.00098515

20,000      47.6900668           0.0023845

30,000      121.3201698        0.00404401

40,000      233.3583473        0.00583396

50,193      363.7078029        0.00724619

100,493    2647.6114347      0.0263462

301,492    32075.94123     .  0.106391

The ultimate speed record to beat, though, is Crandall's 1,048,576 digits in a blistering fast 76.4 hours.

Sam Noble of Apple computed 1,000,000 digits of the MRB constant in 18 days 9 hours 11 minutes 34.253417 seconds

POSTED BY: Marvin Ray Burns
Answer
3 years ago

I found the following program Crandall wrote to check his code. (it gives the MRB constant to the desired level of precision.) It's hard to follow! But here it is:

t1 = 1/2 (2 EulerGamma Log[2] - Log[2]^2); t2 = 
 1/24 (\[Pi]^2 Log[2]^2 - 
    2 \[Pi]^2 Log[
      2] (EulerGamma + Log[2] - 12 Log[Glaisher] + Log[\[Pi]]) - 
    6 (Zeta^\[Prime]\[Prime])[2]);
FLAG = 1;

prec = 200;
$MaxPrecision = Infinity;
expM[pr_] := Module[{a, d, s, k, b, c}, n = Floor[1.3 pr]; Print[n];
   d = N[(3 + Sqrt[8])^n, prec]; d = 1/2 (d + 1/d);

   {b, c, s} = {-1, -d, 0};
   T0 = T00 = SessionTime[];
   Do[c = b - c;
     s = s + 
      c*N[E^(Log[k + 1]/(k + 1)) - 1 - 
         FLAG * (Log[k + 1]/(k + 1) + 1/2! Log[k + 1]^2/(k + 1)^2), 
        pr]; 
    (* s=s+c*N[(k+1)^(1/(k+1))-1,pr]; *)
    b = (k + n) (k - n) b/((k + 1) (k + 1/2));
    If[Mod[k, 1000] == 0, Print[{k, SessionTime[] - T0}];
     T0 = SessionTime[]], {k, 0, n - 1}];
   N[-s/d, pr]];

(*Begin test to given precision.*)
(*t1=Timing[MRBtrue=NSum[(-1)^n*(n^(1/n)-1),{n,1,Infinity},\
WorkingPrecision\[Rule]prec,Method\[Rule]"AlternatingSigns"];];*)

t3 = Timing[MRBtest = expM[prec];]
SessionTime[] - T00
MRBtest +  FLAG * (t1 + t2)
POSTED BY: Marvin Ray Burns
Answer
3 years ago

The MRB constantenter image description here

It is easy to compute the MRB constant using the LHS,

but how about computing it from the RHS, enter image description here?

How many correct digits can you compute using the using it? Here's how one can compute 173 accurate digits using it.

First compute known correct digits. Then compute using the LHS using NSumTerms -> 2^22 and compare.

   c1 = NSum[(-1)^n*(n^(1/n) - 1), {n, 1, Infinity}, 
      Method -> "AlternatingSigns", WorkingPrecision -> 200, 
      PrecisionGoal -> 200];

   c2 = NSum[(2 k)^(1/(2 k)) - (2 k - 1)^(1/(2 k - 1)), {k, 1, Infinity},
       WorkingPrecision -> 600, PrecisionGoal -> 600, NSumTerms -> 2^22];
   c1 - c2

gives -6.6444030376101061743388604*10^-173

While, if you try it naively,

NSum[(-1)^n*(n^(1/n) - 1), {n, 1, \[Infinity]}, 
  Method -> "AlternatingSigns", WorkingPrecision -> 100, 
  PrecisionGoal -> 100] - 
 N[Sum[(2 k)^(1/(2 k)) - (2 k - 1)^(1/(2 k - 1)), {k, 
    1, \[Infinity]}], 100]

you get 8.65025879128817600921262272409653679780852599347660397129681915033728\ 91*10^-28,

with only 28 digits of accuracy!

Attachments:
POSTED BY: Marvin Ray Burns
Answer
3 years ago

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.

POSTED BY: Daniel Lichtblau
Answer
3 years ago

For c1 I used a method for alternating series to get a very precise sum approximation: That is from the "AlternatingSigns" option.

POSTED BY: Marvin Ray Burns
Answer
3 years ago

C1, the approximation to enter image description hereand a special case of enter image description here, is computed correctly by Henri Cohen, Fernando Rodriguez Villegas, and Don Zagier,s Convergence Acceleration of Alternating Series as found next.

enter image description hereenter image description here

Below it is computed by that method to at least 1000 decimals.

  (* Newer loop with Newton interior. *)
 prec = 1000;(*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]

During evaluation of In[17]:= Iterations required: 1320

During evaluation of In[17]:= Iterations: 0 Cumulative time (sec): 0.*10^-8

During evaluation of In[17]:= Iterations: 1000 Cumulative time (sec): 0.1872004

Out[18]= 0. 1878596424620671202485179340542732300559030949001387861720046840894772 3156466021370329665443310749690384234585625801906123137009475922663043 8929348896184120837336626081613602738126379373435283212552763962171489 3217020762820621715167154084126804483635416719985197680252759893899391 4457983505561350964852107120784442309586812949768852694956420425558648 3670441042527952471060666092633974834103115781678641668915460034222258 8380025455396892947114212218910509832871227730802003644521539053639505 5332203470627551159812828039510219264914673176293516190659816018664245 8249506972033819929584209355151625143993576007645932912814517090824249 1588320416906640933443591480670556469280678700702811500938060693813938 5953360657987405562062348704329360737819564603104763950664893061360645 5280675151935082808373767192968663981030949496374962773830498463245634 7931157530028921252329181619562697369707486576547607117801719578736830 0965902260668753656305516567361288150201438756136686552210674305370591 0397357561914891

                 c1 = NSum[(-1)^n*(n^(1/n) - 1), {n, 1, \[Infinity]}, 
                   Method -> "AlternatingSigns", WorkingPrecision -> 100, 
                   PrecisionGoal -> 100]

0.18785964246206712024851793405427323005590309490013878617200468408947 723156466021370329665443310750

POSTED BY: Marvin Ray Burns
Answer
3 years ago

Above where Cohen said,

enter image description here

the coefficients of the numerators are furthered below. I know the pattern for the first and last coefficient of each line, but can you figure out a simple pattern for the rest of them? My hat's off to you if you can!

2
16, 8
98, 80, 32
576, 544, 384, 128
3362, 3312, 2912, 1792, 512
19600, 19528, 18688, 15104, 8192, 2048
114242, 114144, 112576, 103168, 76288, 36864, 8192
665856, 665728, 663040, 641536, 557056, 376832, 163840, 32768
3880898, 3880736, 3876416, 3832064, 3603968, 2945024, 1826816, 720896, 131072 
POSTED BY: Marvin Ray Burns
Answer
2 years ago

One more record, Richard Crandall was the first one to compute 319,000 digits of the MRB constant in one day*, as early as Nov 24, 2012. Here is a portion of his email:

enter image description here

*He may have been counting CPU time alone (from the Timing[] comand which is an under estimate of the actual time taken).. My personal records are in actual time taken.

POSTED BY: Marvin Ray Burns
Answer
2 years ago

OK, you might think I'm real vain for saying the following! Nonetheless, here you go. The ultimate record to break: Richard Mathar and Richard Crandall both wrote 1 scholarly article that talks about the MRB constant to a good extent. Can you break that record and publish a full article or two about the MRB constant? (You might have to come up with some new ideas.) Mathar's and Crandall's papers can be found on Google Scholar: https://scholar.google.com/scholar?q=%22MRB+constant%22 Steven R Finch, Eric W.Weisstein,and i have written some about it in other genres.

POSTED BY: Marvin Ray Burns
Answer
2 years ago

What Richard Crandall and maybe others did to come up with that method is really good and somewhat mysterious. I still don't really understand the inner workings, and I had shown him how to parallelize it. So the best I can say is that it's really hard to compete against magic. (I don't want to discourage others, I'm just explaining why I myself would be reluctant to tackle this. Someone less familiar might actually have a better chance of breaking new ground.)

In a way this should be good news. Should it ever become "easy" to compute, the MRB number would lose what is perhaps its biggest point of interest. It just happens to be on that cusp of tantalizingly "close" to easily computable (perhaps as sums of zeta function and derivatives thereof), yet still hard enough that it takes a sophisticated scheme to get more than a few dozen digits.

POSTED BY: Daniel Lichtblau
Answer
2 years ago

Daniel Lichtblau,
Which one do you not understand? how he derived the eta' formula:enter image description here from enter image description here?

Or is it how and why he added the extra loop to Cohen's method:`

ll = start + l - 2; b *= 2 (ll + n) (ll - n)/((ll + 1) (2 ll + 1)); ?`

I've been asking for a proof of the eta formula derivation. I think it could become the foundation to a great paper.

If the extra loop changes Cohen's method for enter image description here to a method for enter image description hereI think the code in that loop could help prove the derivation.

POSTED BY: Marvin Ray Burns
Answer
2 years ago

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..

POSTED BY: Daniel Lichtblau
Answer
2 years ago

Daniel Lichtblau and others, Richard Crandall did intend to explian his work on the MRB constant and his program to compute it. When I wrote him with a possible small improvement to his program he said, "It's worth observing when we write it up." See screenshot: enter image description here

POSTED BY: Marvin Ray Burns
Answer
2 years ago

Crandall is not using his eta formulas directly!!!!!!! He computes Sum[(-1)^k*(k^(1/k) - 1), {k, 1, Infinity}] directly!

Going back to Crandall's code:

(*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

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 (See a screenshot in a previous post.) to accelerate convergence of Sum[(-1)^k*(k^(1/k) - 1), {k, 1, Infinity}]. That is his secret!!

As I mentioned in a previous post the "MRBtest2 - MRBtest3" is for checking with a known-to-be accurate approximation to the MRB constant, MRBtest3

I'm just excited that I figured it out! as you can tell.

POSTED BY: Marvin Ray Burns
Answer
2 years ago

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

POSTED BY: Daniel Lichtblau
Answer
2 years ago

Richard Crandall might of had some help in developing his method. He wrote one time:

"Marvin I am working on a highly efficient method for your constant, and I've been in touch with other mathematics scholars.

Please be patient...

rec

Sent from my iPhone."

POSTED BY: Marvin Ray Burns
Answer
2 years ago

The result of k^(1/k) from the iteration in Crandall's code grows like the following example of 123^(1/123)..

Timing[ll = 123; x = 1 + 1/10; 
 Table[Block[{$MaxExtraPrecision = Infinity}, 
   x = SetPrecision[x, 20 + a^6]; y = x^ll - ll; 
   x = x (1 - 2 y/((ll + 1) y + 2 ll ll)); 
   N[ll^(1/ll) - x, a + 20]], {a, 1, 13}]]

That gives the following timing and accuracies per iteration.

  {8.829657, {-0.04239426367665370783, -0.02519017014329854597401, \
  -0.0097174389047505410455275, -0.000933400633967433388167272, \
  -9.456045144810702667437515*10^-7, \
  -9.8570302819631224677862339*10^-16, \
  -1.11649387558952636447488240*10^-42, \
  -1.622508976729293063326569833*10^-123, \
  -4.9794273187141078980478269830*10^-366, \
  -1.43931637321260733551982880202*10^-1093, \
  -3.476056847637409626605825363818*10^-3276, \
  -4.8964202322309104109089811118496*10^-9824, \
  -1.36852937628135955635109182098985*10^-29467}}

Here is a more extreme example from my big computer where we get about 64,000,000 digits of accuracy of 123^(1/123) in only 20 iterations!

In[47]:= Timing[ll = 123; x = 1 + 1/10;
 Table[Block[{$MaxExtraPrecision = Infinity}, 
   x = SetPrecision[x, 20 + a^6]; y = x^ll - ll;
   x = x (1 - 2 y/((ll + 1) y + 2 ll ll));
   N[ll^(1/ll) - x, a + 20]], {a, 1, 20}]]

Out[47]= {210.180147, {-0.04239426367665370783, \
-0.02519017014329854597401, -0.0097174389047505410455275, \
-0.000933400633967433388167272, -9.456045144810702667437515*10^-7, \
-9.8570302819631224677862339*10^-16, \
-1.11649387558952636447488240*10^-42, \
-1.622508976729293063326569833*10^-123, \
-4.9794273187141078980478269830*10^-366, \
-1.43931637321260733551982880202*10^-1093, \
-3.476056847637409626605825363818*10^-3276, \
-4.8964202322309104109089811118496*10^-9824, \
-1.36852937628135955635109182098985*10^-29467, \
-2.987998984519890677800153644718129*10^-88398, \
-3.1099929855547029560727524646233719*10^-265190, \
-3.50668133180579791324881778853196092*10^-795566, \
-5.026977910888395698424738167604170974*10^-2386694, \
-1.4809452514436132712292881458571004047*10^-7160077, \
-3.78647491015226246115034702337662410645*10^-21480228, 
  0.*10^-64000020}}

Furthermore this process can compute 6,379,358,787 digits of the Sqrt[2], from an initial guess of machine precision, in 18 iterations:

In[90]:= Timing[ll = 2; x = 2.^(1/2);
 Table[Block[{$MaxExtraPrecision = Infinity}, 
   x = SetPrecision[x, 20 + a^7.5]; y = x^ll - ll;
   x = x (1 - 2 y/((ll + 1) y + 2 ll ll));
   N[ll^(1/ll) - x, a + 20]], {a, 1, 18}]]

Out[90]= {1913.882668, {0.*10^-21, -1.800462644283455490473*10^-148, \
-7.2956225729359623776786*10^-445, \
-4.85397000370369097837099*10^-1334, \
-1.429556345255675471801415*10^-4001, \
-3.6518576944392011682818847*10^-12004, \
-6.08767627470915271802076570*10^-36012, \
-2.820100999650988078782491560*10^-108035, \
-2.8035222068968058313783728204*10^-324105, \
-2.75436831997714080172674040062*10^-972315, \
-2.612017346393607427831858236332*10^-2916945, \
-2.2276049962448401256762035120120*10^-8750835, \
-1.38173437723345455709717486529180*10^-26252505, \
-3.297491628267321222864512567028780*10^-78757516, \
-4.4818892212256882776648121052875660*10^-236272548, \
-1.12536490316593766234907240580868063*10^-708817643, 
  0.*10^-1073741813, 0.*10^-1073741813}}

The above pattern is 3 times the previous number of accurate digits per iteration.. `708817643 * 3^2= 6379358787`

I believe it will only take 191+16 = 207 iterations to compute a google digits of the Sqrt[2]:

In[55]:= Log[10.^100/708817643]/Log[3]

Out[55]= 191.04
POSTED BY: Marvin Ray Burns
Answer
2 years ago

How about computing the MRB constant from Crandall's eta derivative formulas? They are mentioned in a previous post but here they are again:

enter image description here

I computed and checked 500 digits, using the first eta derivative formula in 38.6 seconds. How well can you do? Can you improve my program? (It is a 51.4% improvement of one of Crandall's programs.) I would like a little competition in some of these records! (That formula takes just 225 summands, compared to 10^501 summands using -1^(1/1)+2^(1/2)-3^(1/3)+.... See http://arxiv.org/pdf/0912.3844v3.pdf for more summation requirements for other summation methods.)

In[37]:= mm = 
  0.187859642462067120248517934054273230055903094900138786172004684089\
4772315646602137032966544331074969038423458562580190612313700947592266\
3043892934889618412083733662608161360273812637937343528321255276396217\
1489321702076282062171516715408412680448363541671998519768025275989389\
9391445798350556135096485210712078444230958681294976885269495642042555\
8648367044104252795247106066609263397483410311578167864166891546003422\
2258838002545539689294711421221891050983287122773080200364452153905363\
9505533220347062755115981282803951021926491467317629351619065981601866\
4245824950697203381992958420935515162514399357600764593291281451709082\
4249158832041690664093344359148067055646928067870070281150093806069381\
3938595336065798740556206234870432936073781956460310476395066489306136\
0645528067515193508280837376719296866398103094949637496277383049846324\
5634793115753002892125232918161956269736970748657654760711780171957873\
6830096590226066875365630551656736128815020143875613668655221067430537\
0591039735756191489093690777983203551193362404637253494105428363699717\
0244185516548372793588220081344809610588020306478196195969537562878348\
1233497638586301014072725292301472333336250918584024803704048881967676\
7601198581116791693527968520441600270861372286889451015102919988536905\
7286592870868754254925337943953475897035633134403826388879866561959807\
3351473990256577813317226107612797585272274277730898577492230597096257\
2562718836755752978879253616876739403543214513627725492293131262764357\
3214462161877863771542054231282234462953965329033221714798202807598422\
1065564890048536858707083268874877377635047689160983185536281667159108\
4121934201643860002585084265564350069548328301205461932`1661.\
273833491444;

In[30]:= Timing[
 etaMM[m_, pr_] := 
  Module[{a, d, s, k, b, c}, a[j_] := Log[j + 1]^m/(j + 1)^m;
   n = Floor[1.32 pr];
   d = Cos[n ArcCos[3]];
   {b, c, s} = {-1, -d, 0};
   Do[c = b - c;
    s = s + c a[k];
    b = (k + n) (k - n) b/((k + 1) (k + 1/2)), {k, 0, n - 1}];
   N[s/d, pr] (-1)^m];
 eta[s_] := (1 - 2^(1 - s)) Zeta[s];
 eta1 = Limit[D[eta[s], s], s -> 1];
 MRBtrue = mm;
 prec = 500;
 MRBtest = 
  eta1 - Sum[(-1)^m etaMM[m, prec]/m!, {m, 2, Floor[.45 prec]}];
 MRBtest - MRBtrue]

Out[30]= {36.831836, 0.*10^-502}

Here is a short table of computation times with that program:

Digits      Seconds

500        36.831836
1000       717.308198
1500       2989.759165
2000       3752.354453

I just now retweaked the program. It is now

Timing[etaMM[m_, pr_] := 
  Module[{a, d, s, k, b, c}, 
   a[j_] := N[(-PolyLog[1, -j]/(j + 1))^m, pr];
   n = Floor[1.32 pr];
   d = Cos[n ArcCos[3]];
   {b, c, s} = {-1, -d, 0};
   Do[c = b - c;
    s = s + c a[k];
    b = N[(k + n) (k - n) b/((k + 1) (k + 1/2)), pr], {k, 0, n - 1}];
   Return[N[s/d, pr] (-1)^m]];
 eta[s_] := (1 - 2^(1 - s)) Zeta[s];
 eta1 = Limit[D[eta[s], s], s -> 1];
 MRBtrue = mm;
 prec = 1500;
 MRBtest = 
  eta1 - Sum[(-1)^m etaMM[m, prec]/Gamma[m + 1], {m, 2, 
     Floor[.45 prec]}, Method -> "Procedural"];
 MRBtest - MRBtrue]

Here are my best eta derivative records:

Digits        Seconds
 500          9.874863
 1000        62.587601
 1500        219.41540
 2000       1008.842867
 2500       2659.208646
 3000       5552.902395
 3500       10233.821601

That is using V10.0.2.0 Kernel. Here is a sample

Timing[etaMM[m_, pr_] := 
          Module[{a, d, s, k, b, c}, 
           a[j_] := N[(-PolyLog[1, -j]/(j + 1))^m, pr];
           n = Floor[1.32 pr];
           d = Cos[n ArcCos[3]];
           {b, c, s} = {-1, -d, 0};
           Do[c = b - c;
            s = s + c a[k];
            b = N[(k + n) (k - n) b/((k + 1) (k + 1/2)), pr], {k, 0, n - 1}];
           Return[N[s/d, pr] (-1)^m]];
         eta[s_] := (1 - 2^(1 - s)) Zeta[s];
         eta1 = Limit[D[eta[s], s], s -> 1];
         MRBtrue = mm;
         prec = 500;
         MRBtest = 
          eta1 - Sum[(-1)^m etaMM[m, prec]/Gamma[m + 1], {m, 2, 
             Floor[.45 prec]}];
        ]
         N::meprec: Internal precision limit $MaxExtraPrecision = 50. reached while evaluating 
             -Cos[660 ArcCos[3]].

         N::meprec: Internal precision limit $MaxExtraPrecision = 50. reached while evaluating 
             -Cos[660 ArcCos[3]].

         N::meprec: Internal precision limit $MaxExtraPrecision = 50. reached while evaluating 
             -Cos[660 ArcCos[3]].

         General::stop: Further output of N::meprec will be suppressed during this calculation.

         Out[1]= {9.874863, Null}
POSTED BY: Marvin Ray Burns
Answer
2 years ago

OOPS! I did believe for a little while that I figured out how to rapidly compute AND CHECK a computation of the MRB constant! but it wasn't quite a good check! (The timing given is in processor time [for computing and checking] only. T0 can be used with another SessionTime[] call at the end to figure out all time expired during running of the program.) I used both of Crandall's methods for computing it and used for a check, the nontrivial identityenter image description here ,where gamma is the Euler constant and M is the MRB constant.

Below is my first version of the code with results. If nothing else, I thought, the code pits Crandall's 2 methods against each other to show if one is wrong they both are wrong. (So it is not a real proof.) But these are two totally different methods! (the first of which has been proven by Henry Cohen to be theoretically correct here). For a second check mm is a known approximation to the constant; over 3 million non checked (as of now) digits are found in the attached file 3M.nb. (You will have to change the Format/Style to Input to use the digits.)

I at least did hope that I did more than to only show enter image description here== enter image description here, but I'm afraid I didn't!

In[15]:= (*MRB constant computation with verification! The constant's \
decimal approximation is saved as MRBtest*)prec = 5000;(*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];
  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;;, {k, 0, end - 1}];
  etaMs = N[-s/d - (EulerGamma Log[2] - Log[2]^2/2), pr]]; t2 = 
 Timing[MRBtest2 = expM[prec];];
Print["The MRB constant was computed and checked to ", prec, " digits \
in ", t1 = t2[[1]] + Timing[eta[s_] := (1 - 2^(1 - s)) Zeta[s];
     eta1 = Limit[D[eta[s], s], s -> 1];
     MRBtrue = mm;
     MRBtest = eta1 + etaMs;
     check = MRBtest - MRBtrue][[1]], " seconds"]; check

During evaluation of In[15]:= The MRB constant was computed and checked to 5000 digits in 2.12161 seconds

Out[18]= 0.*10^-5000

In[19]:= MRBtest - mm

Out[19]= 0.*10^-5000
Attachments:
POSTED BY: Marvin Ray Burns
Answer
2 years ago

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.

POSTED BY: Daniel Lichtblau
Answer
2 years ago

Daniel Lichtblau and others, I just deciphered an Identity Crandall used for checking computations of the MRB constant just before he died. It is used in a previous post about checking, where I said it was hard to follow. The MRB constant is B here. B=`enter image description here In input form that is

   B= 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}] + 
     1/24 (\[Pi]^2 Log[2]^2 - 
        2 \[Pi]^2 Log[
          2] (EulerGamma + Log[2] - 12 Log[Glaisher] + Log[\[Pi]]) - 
        6 (Zeta^\[Prime]\[Prime])[2]) + 
     1/2 (2 EulerGamma Log[2] - Log[2]^2)

For 3000 digit numeric approximation, it is

B=NSum[((-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}, 
  Method -> "AlternatingSigns", WorkingPrecision -> 3000] + 
 1/24 (\[Pi]^2 Log[2]^2 - 
    2 \[Pi]^2 Log[
      2] (EulerGamma + Log[2] - 12 Log[Glaisher] + Log[\[Pi]]) - 
    6 (Zeta^\[Prime]\[Prime])[2]) + 
 1/2 (2 EulerGamma Log[2] - Log[2]^2)

It is anylitaclly straight forward too because

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

gives

1/24 (-\[Pi]^2 (Log[2]^2 + EulerGamma Log[4] - 
      24 Log[2] Log[Glaisher] + Log[4] Log[\[Pi]]) - 
   6 (Zeta^\[Prime]\[Prime])[2])

That is enter image description here I wonder why he chose it?

POSTED BY: Marvin Ray Burns
Answer
2 years ago

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.

POSTED BY: Daniel Lichtblau
Answer
2 years ago

If you copy and paste the code for Zeta''[2] you get a non evaluable input!

I get the following results comparing two methods:

In[655]:= 
Timing[MRBtrue = 
   NSum[(-1)^n*(n^(1/n) - 1), {n, 1, Infinity}, 
    WorkingPrecision -> prec, Method -> "AlternatingSigns"]][[1]]

Out[655]= 8.049652

In[656]:= 
Timing[B = 
   NSum[((-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}, 
     Method -> "AlternatingSigns", WorkingPrecision -> 3000] + 
    1/24 (\[Pi]^2 Log[2]^2 - 
       2 \[Pi]^2 Log[
         2] (EulerGamma + Log[2] - 12 Log[Glaisher] + Log[\[Pi]]) - 
       6 (Zeta''[2])) + 1/2 (2 EulerGamma Log[2] - Log[2]^2)][[1]]

Out[656]= 12.277279
POSTED BY: Marvin Ray Burns
Answer
2 years ago

Daniel Lichtblau,

Notice Log[k]/(k) and Log[k]^2/(2k^2) are of the form Log[k]^n/(nk^n). We can take that summing of Log[k]/(k) + Log[k]^2/(2*k^2)+... to its end and write B=enter image description here or compute

NSum[(-1)^(k)*(k^(1/k) - 1 + 
     Sum[Log[k]^n/(n*k^n), {n, 1, Infinity}]), {k, 1, Infinity}, 
  WorkingPrecision -> 1000, Method -> "AlternatingSigns"] - 
 NSum[(-1)^(1 + k) Log[(k - Log[k])/k], {k, 1, Infinity}, 
  WorkingPrecision -> 1000, Method -> "AlternatingSigns"]

Simplified,.B=enter image description here,

What does that do to our "better convergence from brute force summation?"

P.S. I only used 2 terms in my sample data, so I might not have the correct pattern for Log[k]/(k) + Log[k]^2/(2*k^2)+... .

The reason that I suspect that I've got the wrong pattern is, in the following code the summands Log[k]/(k) and Log[k]^2/(2k^2) give 0 where the Log[k]^n/(nk^n).for n>2 does not. (etaMM[m, n] Gives a numeric approximation for the mth derivative of m.)

   etaMM[m_, pr_] := 
  Module[{a, d, s, k, b, c}, a[j_] := Log[j + 1]^m/(j + 1)^m;
   n = Floor[1.32 pr];
   d = Cos[n ArcCos[3]];
   {b, c, s} = {-1, -d, 0};
   Do[c = b - c;
    s = s + c a[k];
    b = (k + n) (k - n) b/((k + 1) (k + 1/2)), {k, 0, n - 1}];
   Return[N[s/d, pr] (-1)^m]]; Table[-Sum[(-1)^m*etaMM[m, 30]/m!, {m, 
      1, a}] - 
   Sum[(-1)^(k)*(Sum[Log[k]^n/(n*k^n), {n, 1, a}]), {k, 1, 
     Infinity}], {a, 1, 15}]
POSTED BY: Marvin Ray Burns
Answer
2 years ago

What about records of computing the integral analog of the MRB constant,: enter image description here?

Richard Mathar did a lot of work on it here.

I've gotten Matheamtica to compute 125 digits. However, they are not proven to be correct yet! They are

0.68765236892769436980931240936544016493963738490362254179507101010743\
366253478493706862729824049846818873192933433546612328629

. First we compute the real part as far as Mathematica will allow.

a1 = NIntegrate[Cos[x Pi] x^(1/x), {x, 1, Infinity}, 
  WorkingPrecision -> 100]

0.07077603931152880353952802183028200136575469620336299759658471973672\
987938741600037745028756981434374

a2 = NIntegrate[Cos[x Pi] x^(1/x), {x, 1, Infinity}, 
  WorkingPrecision -> 120]
a2 - a1

0.07077603931152880353952802183028200136575469620336302758317278266053\
31986618615110244568060496758380620699811570793175408

2.998658806292380331927444551064700651847986149432*10^-53

a3 = NIntegrate[Cos[x Pi] x^(1/x), {x, 1, Infinity}, 
  WorkingPrecision -> 150]
a3 - a2

0.07077603931152880353952802183028200136575469620336302758317278816361\
8457264385970709799491401005081151056924116255307801983594127144525095\
5653544005192

5.5030852586025244596853426853513292430889869429591759902612*10^-63

a4 = NIntegrate[Cos[x Pi] x^(1/x), {x, 1, Infinity}, 
  WorkingPrecision -> 200]
a4 - a3

0.07077603931152880353952802183028200136575469620336302758317278816361\
8457264382036580831881266177238210031756216795402920795214039271485948\
634659563768084109747493815003439875479076850383786911941519465

-3.9341289676101348278429410251678994599048811883800878730391469306948\
367511*10^-78

a5 = NIntegrate[Cos[x Pi] x^(1/x), {x, 1, Infinity}, 
  WorkingPrecision -> 250]
a5 - a4

0.07077603931152880353952802183028200136575469620336302758317278816361\
8457264382036580831881266177238209440733969109717926999044694539086929\
3857095687266500964737783523859835124762555195276023702167529617039725\
7261177753806842756198742365511173334813888

-5.9102224768568499379616934473239901924894999504143401327371546261745\
6363002821330856184541724766503*10^-103

Next we compute the imaginary part to the same precision.

b1 = NIntegrate[I Sin[x Pi] x^(1/x), {x, 1, Infinity}, 
   WorkingPrecision -> 100] - I/Pi

0.*10^-117 - 
 0.6840003894379321291827444599926611267109914826550016181302726087470\
544306934833279937664708191960468 I

b2 = NIntegrate[I Sin[x Pi] x^(1/x), {x, 1, Infinity}, 
   WorkingPrecision -> 120] - I/Pi
b2 - b1

0.*10^-137 - 
 0.6840003894379321291827444599926611267109914826549994343226304054256\
46767722886537984405858512438464223325361496951820797 I

0.*10^-117 + 
 2.1838076422033214076629705967900093606123067575826*10^-51 I

b3 = NIntegrate[I Sin[x Pi] x^(1/x), {x, 1, Infinity}, 
   WorkingPrecision -> 150] - I/Pi
b3 - b2

0.*10^-167 - 
 0.6840003894379321291827444599926611267109914826549994343226303771381\
5305812568206208637713014270949108628424796532117557865488349026349505\
4352728287677 I

0.*10^-137 + 
 2.8287493709597204475898028728369728973137041113531630645218*10^-62 I

b4 = NIntegrate[I Sin[x Pi] x^(1/x), {x, 1, Infinity}, 
   WorkingPrecision -> 200] - I/Pi
b4 - b3

0.*10^-218 - 
 0.6840003894379321291827444599926611267109914826549994343226303771381\
5305812497663815095983421272147867735031056071869477552727290571462108\
208123698276619850397331432861469605963724235550107655309644965 I

0.*10^-167 + 
 7.0542393541729592998801240893393740460248080312761058454887397227149\
1304910*10^-76 I

b5 = NIntegrate[I Sin[x Pi] x^(1/x), {x, 1, Infinity}, 
   WorkingPrecision -> 250] - I/Pi
b5 - b4

0.*10^-268 - 
 0.6840003894379321291827444599926611267109914826549994343226303771381\
5305812497663815095983421272147867223796451609148860995867828496814126\
9810848570299802270095261060286697622600207986034863822997401942304753\
4951409792726050072747412751162199808963072 I

0.*10^-218 + 
 5.1123460446272061655685946207464798122703884124663962338780532683279\
9843703703436946621273009904771*10^-102 I

b6 = NIntegrate[I Sin[x Pi] x^(1/x), {x, 1, Infinity}, 
   WorkingPrecision -> 300] - I/Pi
b6 - b5

0.*10^-318 - 
 0.6840003894379321291827444599926611267109914826549994343226303771381\
5305812497663815095983421272147867223796451609148860995867804988314557\
9408739051911924508290758754789975176921766748245229306743723292030351\
1357229649514450909272015113199881208930542548540913212596310791355732\
04151474091653439098975 I

0.*10^-268 + 
 2.3508499569040210951838787776180450230549672244567844123778963451625\
36786502744023594180143211599163475397637962318600032530*10^-127 I

Notice that WorkingPrecision->100 gave 51 consistant (correct) digits, WorkingPrecision->120 gave 62 correct digits, WorkingPrecision->150 gave 76 correct digits, WorkingPrecision->200 gave 102 correct digits, so it is not too much of a stretch to believe WorkingPrecision->250 gave 125 correct digits.

In[78]:= c = N[Abs[a5 + b5], 125]

Out[78]= 0.\
6876523689276943698093124093654401649396373849036225417950710101074336\
6253478493706862729824049846818873192933433546612328629
POSTED BY: Marvin Ray Burns
Answer
2 years ago

Going back to the Crandall equations for the MRB constant: enter image description here

In particular, this time we will work on equation 44, which should have a negative sign in front.

The naive approach is about as good as any other, so far::

Clear[a, c]; eta[s_] := (1 - 2^(1 - s)) Zeta[s];
a[i_] := Derivative[i][eta][0];
c[j_] := Sum[Binomial[j, d]*(-1)^d*d^(j - d), {d, 1, j}];
-N[Sum[c[m]/m!*a[m], {m, 1, 40}], 6]

giving, after some warning,

0.187859

Nonetheless, I found a shortcut for computing the m'th derivatives of the eta function( with an argument of 0). For an example below you can replace 70 with any other n>3.

Timing[Derivative[70][DirichletEta][0]]]

can be computed in a lot less time by

a = EulerGamma^2/2 - Pi^2/24 - (1/2)*(Log[2] + Log[Pi])^2 + 
   StieltjesGamma[1]; 
b = EulerGamma^3 + (3/2)*EulerGamma^2*Log[2*Pi] - (1/8)*Pi^2*
   Log[2*Pi] - (1/2)*Log[2*Pi]^3 + 3*EulerGamma*StieltjesGamma[1] + 
  3*Log[2*Pi]*StieltjesGamma[1] + 
      (3*StieltjesGamma[2])/2 - 
  Zeta[3]; Timing[((-1)^n*Log[2]^n + (-1)^(n + 1)*n*Log[2]^(n - 1)*
      Log[2*Pi] + (-1)^(n + 1)*n*(n - 1)*Log[2]^(n - 2)*a) + 
        (-1)^n*2*Binomial[n, 3]*Log[2]^(n - 3)*
    b - (Derivative[n][Zeta][0] + 
     Sum[(-1)^x*2*
       Binomial[n, x]*(Derivative[n - x][Zeta][0]*Log[2]^x), {x, 1, 
       n - 4}]) /. n -> 70]

. That shortcut leads to the following raw program to compute the MRB constant. (Unfortunately, it is severely limited in its precision due to Mathematica's refusal to give high precision results for zeta derivatives of order >3 with argument of 0. For example try N[Derivative[4][Zeta][0], 40] ).

Be careful reading the program. (There is a dot product in it. And it could use a lot of cleanup!)

a = EulerGamma^2/2 - Pi^2/24 - (1/2)*(Log[2] + Log[Pi])^2 + 
   StieltjesGamma[1]; 
b = EulerGamma^3 + (3/2)*EulerGamma^2*Log[2*Pi] - (1/8)*Pi^2*
    Log[2*Pi] - (1/2)*Log[2*Pi]^3 + 3*EulerGamma*StieltjesGamma[1] + 
   3*Log[2*Pi]*StieltjesGamma[1] + 
       (3*StieltjesGamma[2])/2 - Zeta[3]; 
  -N[With[{prec = 
     40}, (c = 
        CoefficientList[Series[Exp[(-x)*Exp[x]], {x, 0, prec}], 
         x])[[2 ;; 4]] . 
     Table[Derivative[n][DirichletEta][0], {n, 1, 3}] + 

    c[[5 ;; prec]] . 
     Table[((-1)^n*Log[2]^n + (-1)^(n + 1)*n*Log[2]^(n - 1)*
          Log[2*Pi] + (-1)^(n + 1)*n*(n - 1)*Log[2]^(n - 2)*a) + (-1)^
         n*2*Binomial[n, 3]*Log[2]^(n - 3)*b - 
                (N[Derivative[n][Zeta][0], prec] + 
         Sum[(-1)^x*2*
           Binomial[n, 
            x]*(N[Derivative[n - x][Zeta][0], prec]*Log[2]^x), {x, 1, 
           n - 4}]), {n, 4, prec - 1}]], 40]

Does anyone care to work on this?

The zeta derivatives of order n with argument of 0 do have exact forms not involving the zeta function.(I just haven't figured out the pattern to them.) For examples,

Derivative[6][Zeta][0]

is the same as

Limit[D[2 (2 Pi)^(s - 1) Sin[Pi*s/2] Gamma[1 - s] Zeta[1 - s], {s, 
   6}], s -> 0]

which gives the following that doesn't use zeta.

(1/2688)*(6720*EulerGamma^6 + 5040*EulerGamma^4*Pi^2 + 
   1596*EulerGamma^2*Pi^4 - 275*Pi^6 + 
      32256*EulerGamma^5*Log[2*Pi] + 
   13440*EulerGamma^3*Pi^2*Log[2*Pi] + 
   60480*EulerGamma^4*Log[2*Pi]^2 + 
      10080*EulerGamma^2*Pi^2*Log[2*Pi]^2 - 1596*Pi^4*Log[2*Pi]^2 + 
   53760*EulerGamma^3*Log[2*Pi]^3 + 
      20160*EulerGamma^2*Log[2*Pi]^4 - 1680*Pi^2*Log[2*Pi]^4 - 
   1344*Log[2*Pi]^6 - 
      53760*EulerGamma^3*PolyGamma[2, 1] - 
   80640*EulerGamma^2*Log[2*Pi]*PolyGamma[2, 1] + 
      6720*Pi^2*Log[2*Pi]*PolyGamma[2, 1] + 
   26880*Log[2*Pi]^3*PolyGamma[2, 1] - 13440*PolyGamma[2, 1]^2 + 
      8064*Log[2*Pi]*PolyGamma[4, 1] + 
   40320*EulerGamma^4*StieltjesGamma[1] + 
      20160*EulerGamma^2*Pi^2*StieltjesGamma[1] + 
   3192*Pi^4*StieltjesGamma[1] + 
      161280*EulerGamma^3*Log[2*Pi]*StieltjesGamma[1] + 
   40320*EulerGamma*Pi^2*Log[2*Pi]*StieltjesGamma[1] + 
      241920*EulerGamma^2*Log[2*Pi]^2*StieltjesGamma[1] + 
   20160*Pi^2*Log[2*Pi]^2*StieltjesGamma[1] + 
      161280*EulerGamma*Log[2*Pi]^3*StieltjesGamma[1] + 
   40320*Log[2*Pi]^4*StieltjesGamma[1] - 
      161280*EulerGamma*PolyGamma[2, 1]*StieltjesGamma[1] - 
   161280*Log[2*Pi]*PolyGamma[2, 1]*
        StieltjesGamma[1] + 80640*EulerGamma^3*StieltjesGamma[2] + 
   20160*EulerGamma*Pi^2*StieltjesGamma[2] + 
      241920*EulerGamma^2*Log[2*Pi]*StieltjesGamma[2] + 
   20160*Pi^2*Log[2*Pi]*StieltjesGamma[2] + 
      241920*EulerGamma*Log[2*Pi]^2*StieltjesGamma[2] + 
   80640*Log[2*Pi]^3*StieltjesGamma[2] - 
      80640*PolyGamma[2, 1]*StieltjesGamma[2] + 
   80640*EulerGamma^2*StieltjesGamma[3] + 
      6720*Pi^2*StieltjesGamma[3] + 
   161280*EulerGamma*Log[2*Pi]*StieltjesGamma[3] + 
      80640*Log[2*Pi]^2*StieltjesGamma[3] + 
   40320*EulerGamma*StieltjesGamma[4] + 
      40320*Log[2*Pi]*StieltjesGamma[4] + 8064*StieltjesGamma[5])

In traditional form that result is

enter image description here

Even

Derivative[8][Zeta][0]

is the same as

Limit[D[2 (2 Pi)^(s - 1) Sin[Pi*s/2] Gamma[1 - s] Zeta[1 - s], {s, 
   8}], s -> 0]

gives

-(Pi^8/4608) + (1/6)*
 ... (an exact form not involving zeta).

Derivative[n][Zeta][0] is approximately -n!

Try

In[656]:= ss = 
 Table[N[Limit[
     D[2 (2 Pi)^(s - 1) Sin[Pi*s/2] Gamma[1 - s] Zeta[1 - s], {s, n}],
      s -> 0] + n!, 20], {n, 1, 9}]

Out[656]= {0.081061466795327258220, -0.0063564559085848512101, \
-0.0047111668622544477611, 0.0028968119862920410128, \
-0.00023290755845472453599, -0.00093682513005092950428, \
0.00084982376500166915171, -0.00023243173551155958286, \
-0.00033058966361229644526},

so we can switch them in our ugly program to compute the MRB constant and get a rapidly converging rough approximation as follows.

a = EulerGamma^2/2 - Pi^2/24 - (1/2)*(Log[2] + Log[Pi])^2 + 
   StieltjesGamma[1];
b = EulerGamma^3 + (3/2)*EulerGamma^2*Log[2*Pi] - (1/8)*Pi^2*
    Log[2*Pi] - (1/2)*Log[2*Pi]^3 + 3*EulerGamma*StieltjesGamma[1] + 
   3*Log[2*Pi]*StieltjesGamma[1] + (3*StieltjesGamma[2])/2 - Zeta[3];
-N[With[{prec = 
     50}, (c = 
        CoefficientList[Series[Exp[(-x)*Exp[x]], {x, 0, prec}], 
         x])[[2 ;; 4]].Table[
      Derivative[n][DirichletEta][0], {n, 1, 3}] + 
    c[[5 ;; prec]].Table[((-1)^n*Log[2]^n + (-1)^(n + 1)*n*
          Log[2]^(n - 1)*Log[2*Pi] + (-1)^(n + 1)*n*(n - 1)*
          Log[2]^(n - 2)*a) + (-1)^n*2*Binomial[n, 3]*Log[2]^(n - 3)*
        b - (-n! + 
         Sum[(-1)^x*2*Binomial[n, x]*(-(n - x)!*Log[2]^x), {x, 1, 
           n - 4}]), {n, 4, prec - 1}]], 40]

, which unfortunately is off from the true value of the MRB constant by about 0.00112931. Does anyone have suggestions on how to adjust this new program to give a better rapidly converging approximation to the MRB constant? or how to proceed from here?

POSTED BY: Marvin Ray Burns
Answer
2 years ago

I could really use your opinions here! If your getting tired of my posts let me know.

Power towers are defined below. See http://mathworld.wolfram.com/PowerTower.html .

enter image description here ;;;enter image description here

How about breaking my record for computing the infinite power tower of the MRB constant. Here's how I have to go about computing n digits of it: I compute a decimal expansion of the MRB constant, which I save as m. Then I compute l = Log[m]; N[-ProductLog[-l]/l, n]. My record is 3 million digits, in which I used my 3,014,991 digit computation of the MRB constant (mentioned in the first post of this thread) for m.

It would be EXTREMELY helpful if anyone could find a way to compute the exact power tower without first computing the MRB constant! Such a revelation would be synonymous with finding a closed for solution to the MRB constant! I've been trying to find the exact power tower, lately, with no success. Even a high precision approximation to the power tower would generate a nearly equally precise approximation to the MRB constant. (Let ll be the infinite power tower of m, then ll^(1/ll)=m)

Attached are 3 million digits of the power tower.

Attachments:
POSTED BY: Marvin Ray Burns
Answer
2 years ago

Using the same programs to compute the MRB constant, MMA V10.1 is a little slower than 10.0.0.

Here is a sample of such programs:

(*Fastest (at MRB's end) as of 24 Dec 2014.*)
prec = 10000;(*Number of required decimals.*)ClearSystemCache[];
T0 = SessionTime[];
expM[pre_] := 
  Module[{a, d, s, k, bb, c, n, end, iprec, xvals, x, pc, cores = 12, 
    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["end ", end];
   Print[end*chunksize]; d = Cos[n ArcCos[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;
        h = Log[ll]/ll; x = N[Exp[h], 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" -> 4]];
    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" -> 2];
    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

I'm in the process of tweaking that program to perform better in V 10.1.

I tweaked the program a little. I computed 301493 digits in 31266 seconds, or 8 hours 41 minutes 6 seconds. (See attachment.) It's pretty hard to tweak without loosing accuracy!!!

Some other 10.1 records:

digits        seconds
10,000           9.7
20,000          44
30,000         119
40,000         228
50,193         360
100,493       2573
301,493      31266
Attachments:
POSTED BY: Marvin Ray Burns
Answer
2 years ago

Going back to integral analog of the MRB constant'

enter image description here:

Using formula 5 on page 3 of http://arxiv.org/pdf/0912.3844v3.pdf

.enter image description here

We can compute a great deal of digits of the integral analog of the MRB constant' (I once called it the MKB constant, named after Marsha Kell-Burns my, now ex, wife.) In the paper Mathar simply calls it M1.

Until further notice in this post when we compute the imaginary part of M1, we will be concerned with the imaginary part's absolute value only,

This time we will compute the Imaginary part first to at least 500 digits:

  a[1] = 0; For[n = 1, n < 11, 
    a[n] = N[2/Pi - 
       1/Pi*NIntegrate[
         Cos[Pi*x]*x^(1/x)*(1 - Log[x])/x^2, {x, 1, Infinity}, 
         WorkingPrecision -> 100*n], 50 n]; Print[a[n] - a[n - 1]], 
    n++]; Print[a[11]]
\

giving

0.6840003894379321291827444599926611267109914826549994343226303771381530581249766381509598342127214787

0.*10^-101

0.*10^-151

0.*10^-201

0.*10^-251

0.*10^-301

0.*10^-351

0.*10^-401

0.*10^-451

0.*10^-501

0.6840003894379321291827444599926611267109914826549994343226303771381530581249766381509598342127214786722379645160914886099586780498831455794087390519118879988351918366211827085883779918191195794251385436100844782462528597869421390620796113023053439642582325892202911183326091512210367124716901047132601108752764946385830438156754378694878046808312868541961166205744280461776232345922905313658259576212809654022016030244583148587352474339130505540080799774619683572540292971258866450201101870835703060314349396491402064932644813564545345219868887520120

. Likewise the real part:

b[1] = 0; For[n = 1, n < 11, 
 b[n] = N[-1/Pi*
    NIntegrate[Sin[Pi*x]*x^(1/x)*(1 - Log[x])/x^2, {x, 1, Infinity}, 
     WorkingPrecision -> 100*n], 50 n]; Print[b[n] - b[n - 1]], 
 n++]; Print[b[11]]

giving

0.07077603931152880353952802183028200136575469620336302758317278816361845726438203658083188126617723821

0.*10^-102

0.*10^-152

0.*10^-202

0.*10^-252

0.*10^-302

0.*10^-352

0.*10^-402

0.*10^-452

0.*10^-502

0.07077603931152880353952802183028200136575469620336302758317278816361845726438203658083188126617723820944073396910971792699904464538475364292258443860652193330471222906120205483985764336623434898438270710499897053952312269178485299032185072743545220051257328105422174249313177670295863771714489658779291185716175115405623656039914848817528200250723061535734571065031458992196831648681239079549382556509741967588147362548743205919028695774572411439927516593391029992733107982746794845130889328251307263102570083031527430861023428334369104098217022622689

. Then the magnitude:

N[Sqrt[a[11]^2 + b[11]^2], 500]

giving

0.68765236892769436980931240936544016493963738490362254179507101010743\
3662534784937068627298240498468188731929334335466123286287665409457565\
9577211580255650416284625143925097120589697986500952590195706813170472\
5387265069668971286335322245474865156721299946377659227025219748069576\
0895993932096027520027641920489863095279507385793449828250341732295653\
3809181101532087948181335825805498812728097520936901677028741356923292\
2644964771090329726483682930417491673753430878118054062296678424687465\
624513174205

. That checks with the 200 digits computed by the quadosc command in mpmath by FelisPhasma at https://github.com/FelisPhasma/MKB-Constant .The function is defined here: http://mpmath.googlecode.com/svn/trunk/doc/build/calculus/integration.html#oscillatory-quadrature-quadosc

enter image description here

P.S.

I just now finished 750 digits, (about the max with formula 5 from the paper, as far as Mathematica is concerned).

Here is the work:

a[1] = 0; For[n = 1, n < 16, 
 a[n] = N[2/Pi - 
    1/Pi*NIntegrate[
      Cos[Pi*x]*x^(1/x)*(1 - Log[x])/x^2, {x, 1, Infinity}, 
      WorkingPrecision -> 100*n], 50 n]; Print[a[n] - a[n - 1]], 
 n++]; Print[a[16]]; 
b[1] = 0; For[n = 1, n < 16, 
 b[n] = N[-1/Pi*
    NIntegrate[Sin[Pi*x]*x^(1/x)*(1 - Log[x])/x^2, {x, 1, Infinity}, 
     WorkingPrecision -> 100*n], 50 n]; Print[b[n] - b[n - 1]], 
 n++]; Print[b[16]]; Print[N[Sqrt[a[16]^2 + b[16]^2], 750]]

0.6840003894379321291827444599926611267109914826549994343226303771381530581249766381509598342127214787

0.*10^-101

0.*10^-151

0.*10^-201

0.*10^-251

0.*10^-301

0.*10^-351

0.*10^-401

0.*10^-451

0.*10^-501

0.*10^-551

0.*10^-601

3.*10^-650

-4.*10^-700

-2.6*10^-749

0.68400038943793212918274445999266112671099148265499943432263037713815\
3058124976638150959834212721478672237964516091488609958678049883145579\
4087390519118879988351918366211827085883779918191195794251385436100844\
7824625285978694213906207961130230534396425823258922029111833260915122\
1036712471690104713260110875276494638583043815675437869487804680831286\
8541961166205744280461776232345922905313658259576212809654022016030244\
5831485873524743391305055400807997746196835725402929712588664502011018\
7083570306031434939649140206493264481356454534521986888752011950353818\
1776359577265099302389566135475579468144849763261779452665955246258699\
8679271659049208654746533234375478909962633090080006358213908728990850\
5026759549928935029206442637425786005036048098598304092996753145589012\
64547453361707037686708654522699


0.07077603931152880353952802183028200136575469620336302758317278816361845726438203658083188126617723821

0.*10^-102

0.*10^-152

0.*10^-202

0.*10^-252

0.*10^-302

0.*10^-352

0.*10^-402

0.*10^-452

0.*10^-502

2.*10^-551

-1.*10^-600

1.8*10^-650

1.27*10^-699

4.34*10^-749

0.07077603931152880353952802183028200136575469620336302758317278816361\
8457264382036580831881266177238209440733969109717926999044645384753642\
9225844386065219333047122290612020548398576433662343489843827071049989\
7053952312269178485299032185072743545220051257328105422174249313177670\
2958637717144896587792911857161751154056236560399148488175282002507230\
6153573457106503145899219683164868123907954938255650974196758814736254\
8743205919028695774572411439927516593391029992733107982746794845130889\
3282513072631025700830315274308610234283343691040982170226226904594029\
7055093272952022662549075225941956559080574835998923469310063614655255\
0629713179601483134045038416878054929072981851045829413286377842843667\
5378730394247519728064887287780998671021887797977772522419765594172569\
277490031071938177749184834961300

0.687652368927694369809312409365440164939637384903622541795071010107433662534784937068627298240498468188731929334335466123286287665409457565957721158025565041628462514392509712058969798650095259019570681317047253872650696689712863353222454748651567212999463776592270252197480695760895993932096027520027641920489863095279507385793449828250341732295653380918110153208794818133582580549881272809752093690167702874135692329226449647710903297264836829304174916737534308781180540622966784246874656245131742049004832216427665542900559350289936114782223424261285828326467186036500189315374147638489679365569122714398706519530651330568884655048857998738535162606116788633540389660052822237449082894798620397228331715198160243676576563833057235963591510865254600

Using formula 7 from http://arxiv.org/pdf/0912.3844v3.pdf,

enter image description here .

(Treating it as we did formula 5), First, the imaginary part to at least 1000 digits::

a[1] = 0; For[n = 1, n < 21, 
 a[n] = N[2/Pi + 
    1/Pi^2 NIntegrate[
      Sin[x Pi] x^(1/x) (1 - 3 x + 2 (x - 1) Log[x] + Log[x]^2)/
        x^4, {x, 1, Infinity}, WorkingPrecision -> 100 n], 50 n];
 Print[a[n] - a[n - 1]], n++]; Print[a[21]]

0.6840003894379321291827444599926611267109914826549994343226303771381530581249766381509598342127214787

0.*10^-101

0.*10^-151

0.*10^-201

0.*10^-251

0.*10^-301

0.*10^-351

0.*10^-401

0.*10^-451

0.*10^-501

0.*10^-551

0.*10^-601

0.*10^-651

0.*10^-701

0.*10^-751

0.*10^-801

0.*10^-851

0.*10^-901

-2.*10^-950

5.*10^-1000

0.684000389437932129182744459992661126710991482654999434322630377138153058124976638150959834212721478672237964516091488609958678049883145579408739051911887998835191836621182708588377991819119579425138543610084478246252859786942139062079611302305343964258232589220291118332609151221036712471690104713260110875276494638583043815675437869487804680831286854196116620574428046177623234592290531365825957621280965402201603024458314858735247433913050554008079977461968357254029297125886645020110187083570306031434939649140206493264481356454534521986888752011950353818177635957726509930238956613547557946814484976326177945266595524625869986792716590492086547465332343754789099626330900800063582139087289908505026759549928935029206442637425786005036048098598304092996753145589012645474533617070376867086545228223060940434935219252885333298390272342234952870883304116640409421452765284609364941205344122569781634782508368641126766528707019957340895061936246645065753101916781254557006989818409283317145837167345971516970849116096077030635788389165381066055992688

Then the real part to at least 1000 digits:

b[1] = 0; For[n = 1, n < 21, 
 b[n] = N[1/Pi^2 - 
    1/Pi^2 NIntegrate[
      Cos[Pi x] x^(1/x) (1 - 3 x + 2 (x - 1) Log[x] + Log[x]^2)/
        x^4, {x, 1, Infinity}, WorkingPrecision -> 100 n], 50 n];
 Print[b[n] - b[n - 1]], n++]; Print[b[21]]

0.07077603931152880353952802183028200136575469620336302758317278816361845726438203658083188126617723821

0.*10^-102

0.*10^-152

0.*10^-202

0.*10^-252

0.*10^-302

0.*10^-352

0.*10^-402

0.*10^-452

0.*10^-502

0.*10^-552

0.*10^-602

0.*10^-652

0.*10^-702

0.*10^-752

0.*10^-802

0.*10^-852

-3.*10^-901

8.*10^-951

-4.6*10^-1000

0.0707760393115288035395280218302820013657546962033630275831727881636184572643820365808318812661772382094407339691097179269990446453847536429225844386065219333047122290612020548398576433662343489843827071049989705395231226917848529903218507274354522005125732810542217424931317767029586377171448965877929118571617511540562365603991484881752820025072306153573457106503145899219683164868123907954938255650974196758814736254874320591902869577457241143992751659339102999273310798274679484513088932825130726310257008303152743086102342833436910409821702262269045940297055093272952022662549075225941956559080574835998923469310063614655255062971317960148313404503841687805492907298185104582941328637784284366753787303942475197280648872877809986710218877979777725224197655941725692774900310719381777491848349627938468198411955193898347075098152638657614980900350262780319142430252921925131515239611841070722530473939496294305264627977744876814858325335947117076721493110160508928494597906728688873533031986215124467678736429981544321187124269147141804397293341613

Then the magnitude:

In[97]:= N[Sqrt[a[21]^2 + b[21]^2], 1000]

Out[97]= 0.\
6876523689276943698093124093654401649396373849036225417950710101074336\
6253478493706862729824049846818873192933433546612328628766540945756595\
7721158025565041628462514392509712058969798650095259019570681317047253\
8726506966897128633532224547486515672129994637765922702521974806957608\
9599393209602752002764192048986309527950738579344982825034173229565338\
0918110153208794818133582580549881272809752093690167702874135692329226\
4496477109032972648368293041749167375343087811805406229667842468746562\
4513174204900483221642766554290055935028993611478222342426128582832646\
7186036500189315374147638489679365569122714398706519530651330568884655\
0488579987385351626061167886335403896600528222374490828947986203972283\
3171519816024367657656383305723596359151086525460036387486837632622334\
2987257095524637683005910353149353985736118868884201748241906260834981\
7303422370398413326428269921074045506558966667483453656748906071577744\
4147548424388220133662816274116986724576330176058912438027319979840883\
05950589130911719199

PPS. I just now finished a 1500 digit computation of the integral analog of the MRB constant, but I don't have any way of checking it other than to see that it confirms smaller computations. Which thing it does.

In[99]:= aa = 
 N[2/Pi + 1/Pi^2 NIntegrate[
     Sin[x Pi] x^(1/x) (1 - 3 x + 2 (x - 1) Log[x] + Log[x]^2)/
       x^4, {x, 1, Infinity}, WorkingPrecision -> 3000], 1500]

Out[99]= 0.\
6840003894379321291827444599926611267109914826549994343226303771381530\
5812497663815095983421272147867223796451609148860995867804988314557940\
8739051911887998835191836621182708588377991819119579425138543610084478\
2462528597869421390620796113023053439642582325892202911183326091512210\
3671247169010471326011087527649463858304381567543786948780468083128685\
4196116620574428046177623234592290531365825957621280965402201603024458\
3148587352474339130505540080799774619683572540292971258866450201101870\
8357030603143493964914020649326448135645453452198688875201195035381817\
7635957726509930238956613547557946814484976326177945266595524625869986\
7927165904920865474653323437547890996263309008000635821390872899085050\
2675954992893502920644263742578600503604809859830409299675314558901264\
5474533617070376867086545228223060940434935219252885333298390272342234\
9528708833041166404094214527652846093649412053441225697816347825083686\
4112676652870701995734089506193624664506575310191678125455700698981840\
9283317145837167345971516970849116096077030635788389165381066055992708\
4284702473154303800276803908560080204997803241058414188902018357202062\
9532415382916822796942734253441520784640814155687968986766443021927163\
6249354786973717955004441549085673392105556692081075647388204227896978\
1483978754685921758294318270385312597177598977912650715548994562461701\
1553879109152932039370312241134127950112036269188660519350584627066913\
4925878278209048717316088629321353274101519307401594635990058104175474\
300641475776727955287474213040

In[98]:= bb = 
 N[1/Pi^2 - 
   1/Pi^2 NIntegrate[
     Cos[Pi x] x^(1/x) (1 - 3 x + 2 (x - 1) Log[x] + Log[x]^2)/
       x^4, {x, 1, Infinity}, WorkingPrecision -> 3000], 1500]

Out[98]= 0.\
0707760393115288035395280218302820013657546962033630275831727881636184\
5726438203658083188126617723820944073396910971792699904464538475364292\
2584438606521933304712229061202054839857643366234348984382707104998970\
5395231226917848529903218507274354522005125732810542217424931317767029\
5863771714489658779291185716175115405623656039914848817528200250723061\
5357345710650314589921968316486812390795493825565097419675881473625487\
4320591902869577457241143992751659339102999273310798274679484513088932\
8251307263102570083031527430861023428334369104098217022622690459402970\
5509327295202266254907522594195655908057483599892346931006361465525506\
2971317960148313404503841687805492907298185104582941328637784284366753\
7873039424751972806488728778099867102188779797777252241976559417256927\
7490031071938177749184834962793846819841195519389834707509815263865761\
4980900350262780319142430252921925131515239611841070722530473939496294\
3052646279777448768148583253359471170767214931101605089284945979067286\
8887353303198621512446767873642998154432118712426914714180439729334146\
8345902382977472975053271988386946291215512340931334841526712825988330\
6521193975174379922254198045615178994412133135553490942451521573377205\
4086429300485891441696490339106907723915822537813700713422515725943626\
7756749980892097547020923938358076198570370106085596863039832425037481\
4946826330552459256977035009973219582010379262683780372730214991685800\
3676611833579648850161974289307066295385292264148146789532534018500663\
1153014589399140567464592864024

In[109]:= c1500 = Sqrt[aa^2 + bb^2]

Out[109]= \
0.68765236892769436980931240936544016493963738490362254179507101010743\
3662534784937068627298240498468188731929334335466123286287665409457565\
9577211580255650416284625143925097120589697986500952590195706813170472\
5387265069668971286335322245474865156721299946377659227025219748069576\
0895993932096027520027641920489863095279507385793449828250341732295653\
3809181101532087948181335825805498812728097520936901677028741356923292\
2644964771090329726483682930417491673753430878118054062296678424687465\
6245131742049004832216427665542900559350289936114782223424261285828326\
4671860365001893153741476384896793655691227143987065195306513305688846\
5504885799873853516260611678863354038966005282223744908289479862039722\
8331715198160243676576563833057235963591510865254600363874868376326223\
3429872570955246376830059103531493539857361188688842017482419062608349\
8173034223703984133264282699210740455065589666674834536567489060715777\
4441475484243882201336628162741169867245763301760589124380273199798408\
8305950589130911719198776146941477264898934365742508503405073273852990\
3546587114217499635584514475429656959327732862489935076490012861232249\
2446704232200904844779690044774489466704342791971033325818579375177198\
9865742583276770011926585495711579480114327818546199372349313180236079\
1389248808154759564302727311223193005229640892474022665093207969297797\
9723087954832182561714039165214592519432072341006090867558444590500046\
6707963346545638317950978935794173691635274461184852166407791838662429\
40408834876470623546535579027725
POSTED BY: Marvin Ray Burns
Answer
2 years ago

Mathar gives a simple scheme to find better formulas at http://arxiv.org/pdf/0912.3844v3.pdf . I could use some help in programming it: (I keep getting erroneous results!) Does anyone get the right results here?

enter image description here

POSTED BY: Marvin Ray Burns
Answer
2 years ago

Below, where the upper limit of the following integrals shows Infinity, it is meant to be the (Ultraviolet limit of the sequence) as mentioned by Mathar here:

enter image description here

Until further notice in this post when we compute the imaginary part of M1, we will be concerned with the imaginary part's absolute value only,

I derived a new formula for computing the integral analog of the MRB constant':

f[x_]:=x^(1/x);-((2 I)/\[Pi]^3) + 1/\[Pi]^2 - (
 2 I)/\[Pi] + (I/Pi)^3*
  Integrate[(-1)^x*D[f[x], {x, 3}], {x, 1, Infinity}]

In traditional form that is M1= enter image description here

Using it I computed 2000 digits in only 10.8 minutes:

In[131]:= Timing[f[x_] = x^(1/x); 
 a = N[1/\[Pi]^2 + (1/Pi)^3*
     NIntegrate[Sin[Pi*x]*D[f[x], {x, 3}], {x, 1, Infinity}, 
      WorkingPrecision -> 4000], 2000]; 
 b = N[2/\[Pi]^3 + 
    2/\[Pi] + (1/Pi)^3*
     NIntegrate[Cos[Pi x]*D[f[x], {x, 3}], {x, 1, Infinity}, 
      WorkingPrecision -> 4000], 2000]; 
 Print[N[Sqrt[a^2 + b^2], 2000]]]

During evaluation of In[131]:= 0.68765236892769436980931240936544016493963738490362254179507101010743366253478493706862729824049846818873192933433546612328628766540945756595772115802556504162846251439250971205896979865009525901957068131704725387265069668971286335322245474865156721299946377659227025219748069576089599393209602752002764192048986309527950738579344982825034173229565338091811015320879481813358258054988127280975209369016770287413569232922644964771090329726483682930417491673753430878118054062296678424687465624513174204900483221642766554290055935028993611478222342426128582832646718603650018931537414763848967936556912271439870651953065133056888465504885799873853516260611678863354038966005282223744908289479862039722833171519816024367657656383305723596359151086525460036387486837632622334298725709552463768300591035314935398573611886888420174824190626083498173034223703984133264282699210740455065589666674834536567489060715777444147548424388220133662816274116986724576330176058912438027319979840883059505891309117191987761469414772648989343657425085034050732738529903546587114217499635584514475429656959327732862489935076490012861232249244670423220090484477969004477448946670434279197103332581857937517719898657425832767700119265854957115794801143278185461993723493131802360791389248808154759564302727311223193005229640892474022665093207969297797972308795483218256171403916521459251943207234100609086755844459050004667079633465456383179509789357941736916352744611848521664077918386624294040883487647062354653558109265769644276994369741555722263494599492834558291937955573706480722982389806312472239746286527176248883116124285469947303667188075506826507811479428582807366599407544908560990699866167233307144245764835741501174979679166078765231145175411199825822532170091858833628202128777966026600647843068442894310401343003939117236867245656732686719139206716028255819141802331701942027248337771633882445225049334329008827371320849006472846226868011129149192754883153995560921671208059671732704499253517327447921147157

Out[131]= {653.145, Null}

I am presently computing 10,000 digits using that formula. Come back here for results!

That formula didn't work out; I will try one of the following formulas.

Here are 2 more, more advanced formulas; remember f(x) is x^(1/x):

enter image description here

I did finish a 5,000 digit computation using M1=

enter image description here

in 48.11 minutes.

Here are the 5000 digits:of the magnitude:

0.68765236892769436980931240936544016493963738490362254179507101010743366253478493706862729824049846818873192933433546612328628766540945756595772115802556504162846251439250971205896979865009525901957068131704725387265069668971286335322245474865156721299946377659227025219748069576089599393209602752002764192048986309527950738579344982825034173229565338091811015320879481813358258054988127280975209369016770287413569232922644964771090329726483682930417491673753430878118054062296678424687465624513174204900483221642766554290055935028993611478222342426128582832646718603650018931537414763848967936556912271439870651953065133056888465504885799873853516260611678863354038966005282223744908289479862039722833171519816024367657656383305723596359151086525460036387486837632622334298725709552463768300591035314935398573611886888420174824190626083498173034223703984133264282699210740455065589666674834536567489060715777444147548424388220133662816274116986724576330176058912438027319979840883059505891309117191987761469414772648989343657425085034050732738529903546587114217499635584514475429656959327732862489935076490012861232249244670423220090484477969004477448946670434279197103332581857937517719898657425832767700119265854957115794801143278185461993723493131802360791389248808154759564302727311223193005229640892474022665093207969297797972308795483218256171403916521459251943207234100609086755844459050004667079633465456383179509789357941736916352744611848521664077918386624294040883487647062354653558109265769644276994369741555722263494599492834558291937955573706480722982389806312472239746286527176248883116124285469947303667188075506826507811479428582807366599407544908560990699866167233307144245764835741501174979679166078765231145175411199825822532170091858833628202128777966026600647843068442894310401343003939117236867245656732686719139206716028255819141802331701942027248337771633882445225049334329008827371320849006472846226868011129149192754883153995560921671208059671732704499253517327447529208297180672654123457301218758892278525894167935930983363218877512533994251978272092700003994136520699813263053327399132641690231179063314931546906927612775633995348209911166678724589467821767106592498663827057034363632241807121831546175498178011687284590439293322231263406301066863589072717290630291441982684113819198880100231182613587798104863611185433976009254862585527222843445901958943153561148829083242874018226480554274231391324767376148485531787767908124831873688579979114662856184612164534836370699371440464263768724668291617743681719766849740663590277737977490693183461320266666793472116774276618408124767965369796362732668987556797338128876129264558867657737417548617146808592137056879602982206609613881069490166381528825180204703315896719667069923077454352649723496033985893188309150391579573916059639453655188856334980355047281560296288150836680499821806918067869468571687709518088408966653716009356556714281694904914038988996962213833530636987279769672200413448893419914190954063100962251649102614676944333201213024711868954772741991675045198246947499574872027800654821823797116399297131866662866832215332914761325880983081211272181775518951539503852063119472301382766303820851467743266039356123495461914463960644386394228342211998370152351720235034997434035743513051754761571835043769475528640144621307760159481496713401409374957729200400650100318226988524015127382509490642900236553851499823658269458873976032051355393161653806016080446394196719312454167915154602448638624354575153334932298393406734174580316934939632892851077461038399470015366439910136971186909599331204517462262508377673477745789645309425145559198802530351403897927622891172233239135167420567162398873965477371498335087310395422796362380227536212159184529243644094285328763286873653399867593200891823468738537356817916009007206857590792983184556882143118383332812491747733056313117179696094921120670802012310012864110800437831852620698327457619035904268498030693438632685623213366864129523404256345542376567721287706234359125016588483777876970236084456277023948551334490591022594253744077631232660869593809453087749830900393202787736482133628148979992109544954840067942735030391105496026321872468122542495017023785810605820545392820104069279893067324597299043883381251767370331206913429284614563732308018369972360638019778425246546329838131639355043236388708044857300408692365733932897876809202025693305332974091411983635619038514442263783801745983300121464879550146672827072002317686396598587702487509572349422593441184802476344187280014450860069307120621758277552124841158659386176036703247124389223327008210072318671884895179305778728051888524412158486781863155034447221379906386062559915129172725833420555901857729690605950941678587057025641848365090809750870051863842805803189784976076099574956436664131457150096711473033060684065060747340764998195621425524824611657787212347497307297184843276100338110267863618974154272345482369968216663233417338501929114697679974461999040589290327155974468087040862022522065912789

I'm getting closer to 10K digits of M1: Using enter image description here , where f(x)=x^(1/x).

I got approx. 10K digits of the imaginary part, but the real part was a little garbled.

Finally using

enter image description here , where f(x) = x^(1/x) ,

.I got about 10000 digits of M1 in about 12 hours. (It showed that the 5000 digit computation was only correct to 4979 digits, though.) Here is a rough program to get it:

   f[x_]:=x^(1/x);  Print[DateString[]]; Print[T0 = SessionTime[]]; prec = 10000; 
  Timing[Print[
  a = N[Re[-(136584/Pi^10) - (34784*I)/Pi^9 + 
       5670/Pi^8 + (786*I)/Pi^7 - 90/Pi^6 - 
                 (4*I)/Pi^5 - 3/Pi^4 - (2*I)/Pi^3 + 
       1/Pi^2 - (2*I)/Pi] - 
             (1/Pi)^10*
      NIntegrate[Cos[Pi*x]*D[f[x], {x, 10}], {x, 1, Infinity}, 
       WorkingPrecision -> prec, 
                 PrecisionGoal -> prec], prec]]; 
 Print[SessionTime[] - T0, " seconds"]; 
     Print[
  N[b = -Im[-(136584/Pi^10) - (34784*I)/Pi^9 + 
        5670/Pi^8 + (786*I)/Pi^7 - 90/Pi^6 - 
                   (4*I)/Pi^5 - 3/Pi^4 - (2*I)/Pi^3 + 
        1/Pi^2 - (2*I)/Pi] + 
             (1/Pi)^10*
      NIntegrate[Sin[Pi*x]*D[f[x], {x, 10}], {x, 1, Infinity}, 
       WorkingPrecision -> prec, 
                 PrecisionGoal -> prec], prec]]]; Print[
 SessionTime[] - T0, " seconds"]; 
  Print[N[Sqrt[a^2 + b^2], prec]]; Print[DateString[]]; 

See attached 10000MKB.pdf and 10KMKB.nb for work and digits.

On May 5, I computed another 10,000 digits in 9.55 hours see attached faster10KMKB.

On May 6, I computed another 10,000 digits in a blistering fast 5.1 hours see attached fastest10KMKB.nb.

On May 9, I improved that timing to 4.8 hours (17355 seconds). Here is the code I used:

d = 15; f[x_] = x^(1/x); ClearAll[a, b, h];
h[n_] := Sum[
  StirlingS1[n, k]*Sum[(-j)^(k - j)*Binomial[k, j], {j, 0, k}], {k, 1,
    n}]; h[0] = 1; g = 
 2 I/Pi - Sum[-I^(n + 1) h[n]/Pi^(n + 1), {n, 1, d}]; Print[
 DateString[]];
Print[T0 = SessionTime[]]; prec = 10000;
Print[N[a = -Re[g] + (1/Pi)^(d + 1)*
      NIntegrate[
       Simplify[Cos[Pi*x]*D[f[x], {x, d + 1}]], {x, 1, Infinity}, 
       WorkingPrecision -> prec*(105/100), 
       PrecisionGoal -> prec*(105/100)], prec]];
Print[SessionTime[] - T0, " seconds"];
Print[N[b = 
    Im[g] - (1/Pi)^(d + 1)*
      NIntegrate[
       Simplify[Sin[Pi*x]*D[f[x], {x, d + 1}]], {x, 1, Infinity}, 
       WorkingPrecision -> prec*(105/100), 
       PrecisionGoal -> prec*(105/100)], prec]];
Print[SessionTime[] - T0, " seconds"]; Print[
 N[Sqrt[a^2 + b^2], prec]]; Print[DateString[]];
POSTED BY: Marvin Ray Burns
Answer
2 years ago

FelisPhasma has been helpful in providing me with a little competition in computing the integral analog of the MRB constant. See https://github.com/FelisPhasma/MKB-Constant.

I've never done this before. But I so much would like to see others breaking these records that I'm going to give away a program that is practically guaranteed to break my record of 10,000 digits, for the integral analog of the MRB constant in a day or so. The program could use some "clean up" if you care to go that far. (The imaginary part is given as a positive, real constant: it actually starts with a negative sign and of course ends with I.)

Here it is:

f[x_] = x^(1/x); Print[DateString[]]; Print[
 T0 = SessionTime[]]; prec = 11000; Timing[
 Print[a = 
   N[Re[(633666648 I)/\[Pi]^13 - 
       33137280/\[Pi]^12 - ((824760 I)/\[Pi]^11) - 
       136584/\[Pi]^10 - (34784 I)/\[Pi]^9 + 
       5670/\[Pi]^8 + (786 I)/\[Pi]^7 - 90/\[Pi]^6 - (4 I)/\[Pi]^5 - 
       3/\[Pi]^4 - (2 I)/\[Pi]^3 + 
       1/\[Pi]^2 - (2 I)/\[Pi]] + (1/Pi)^12*
      NIntegrate[Cos[Pi x]*D[f[x], {x, 12}], {x, 1, Infinity}, 
       WorkingPrecision -> prec, PrecisionGoal -> prec], prec]];
 Print[SessionTime[] - T0, " seconds"];
 Print[N[b = -Im[(633666648 I)/\[Pi]^13 - 
        33137280/\[Pi]^12 - ((824760 I)/\[Pi]^11) - 
        136584/\[Pi]^10 - (34784 I)/\[Pi]^9 + 
        5670/\[Pi]^8 + (786 I)/\[Pi]^7 + 90/\[Pi]^6 - (4 I)/\[Pi]^5 - 
        3/\[Pi]^4 - (2 I)/\[Pi]^3 + 
        1/\[Pi]^2 - (2 I)/\[Pi]] - (1/Pi)^13*
      NIntegrate[Cos[Pi x]*D[f[x], {x, 13}], {x, 1, Infinity}, 
       WorkingPrecision -> prec, PrecisionGoal -> prec], 
   prec]]]; Print[SessionTime[] - T0, " seconds"]; Print[
 N[Sqrt[a^2 + b^2], prec]]; Print[DateString[]];

Will anyone let me know you are running this program to break my record?

Edit: On Sat 2 May 2015 19:03:45 I started a 15,000 digit, new record computation of the real and imaginary parts and magnitude of the integral analog of the MRB constant, (where the imaginary part is given as a positive, real constant), using the following code.

 f[x_] = x^(1/x); ClearAll[a];
   h[n_] := Sum[
     StirlingS1[n, k]*Sum[(-j)^(k - j)*Binomial[k, j], {j, 0, k}], {k, 1,
       n}]; h[0] = 1; g = -2 I/Pi + 
     Sum[-I^(n + 1) h[n]/Pi^(n + 1), {n, 1, 18}]; Print[DateString[]];
   Print[T0 = SessionTime[]]; prec = 15000;
   Print[N[a = 
       Re[g] + (1/Pi)^19*
         NIntegrate[
          Simplify[Sin[Pi*x]*D[f[x], {x, 19}]], {x, 1, Infinity}, 
          WorkingPrecision -> prec*(105/100), 
          PrecisionGoal -> prec*(105/100)], prec]];
   Print[SessionTime[] - T0, " seconds"];
   Print[N[b = -Im[g] + (1/Pi)^19*
         NIntegrate[
          Simplify[Cos[Pi*x]*D[f[x], {x, 19}]], {x, 1, Infinity}, 
          WorkingPrecision -> prec*(105/100), 
          PrecisionGoal -> prec*(105/100)], prec]];
   Print[SessionTime[] - T0, " seconds"]; Print[
    N[Sqrt[a^2 + b^2], prec]]; Print[DateString[]];

The formula behind this computation is enter image description here

Edit: The program took 33.75 hours,The full run is attached in 15KMKB3.nb.

Edit May 9, 2015: I better than halved my time! I computed 15000 digits in 14.83 hours. See fastestMKB15K.nb/. The faster formula is

enter image description here

If you still want me to write out a code for more digits, for you to break that record, let me know.

POSTED BY: Marvin Ray Burns
Answer
2 years ago

Still talking about the integral analog of the MRB constant:enter image description here

Here are my speed records -- can you beat any of them?

enter image description here

Here is a graph of those speed records with a trendline:

enter image description here

The 20,000 digit run is attached as MKB20K.nb, and MKB20K.pdf,

Here is the algorithm used:

enter image description here

Here is the code:

d = 30; f[x_] = x^(1/x); ClearAll[a, b, h];
a[n_] := Sum[
  StirlingS1[n, k]*Sum[(-j)^(k - j)*Binomial[k, j], {j, 0, k}], {k, 1,
    n}]; a[0] = 1; g = 
 2 I/Pi - Sum[-I^(n + 1) a[n]/Pi^(n + 1), {n, 1, d}]; Print[
 DateString[]];
Print[T0 = SessionTime[]]; prec = 20000;
Print[N[a = -Re[g] - (1/Pi)^(d)*
      NIntegrate[
       Simplify[Cos[Pi*x]*D[f[x], {x, d}]], {x, 1, Infinity}, 
       WorkingPrecision -> prec*(105/100), 
       PrecisionGoal -> prec*(105/100)], prec]];
Print[SessionTime[] - T0, " seconds"];
Print[N[b = 
    Im[g] + (1/Pi)^(d + 1)*
      NIntegrate[
       Simplify[Cos[Pi*x]*D[f[x], {x, d + 1}]], {x, 1, Infinity}, 
       WorkingPrecision -> prec*(105/100), 
       PrecisionGoal -> prec*(105/100)], prec]];
Print[SessionTime[] - T0, " seconds"]; Print[
 N[Sqrt[a^2 + b^2], prec]]; Print[DateString[]];

I just now completed a 25,000 digit computation. It took 63.7 hours and confirmed the 20,000 digits. I updated MKB20K.nb and MKB20K.pdf. Here is the algorithm and the code I used:

enter image description here

d = 35; f[x_] = x^(1/x); ClearAll[a, b, h];
h[n_] := Sum[
  StirlingS1[n, k]*Sum[(-j)^(k - j)*Binomial[k, j], {j, 0, k}], {k, 1,
    n}]; h[0] = 1; g = 
 2 I/Pi - Sum[-I^(n + 1) h[n]/Pi^(n + 1), {n, 1, d}]; Print[
 DateString[]];
Print[T0 = SessionTime[]]; prec = 25000;
Print[N[a = -Re[g] + (1/Pi)^(d + 1)*
      NIntegrate[
       Simplify[Cos[Pi*x]*D[f[x], {x, d + 1}]], {x, 1, Infinity}, 
       WorkingPrecision -> prec*(105/100), 
       PrecisionGoal -> prec*(105/100)], prec]];
Print[SessionTime[] - T0, " seconds"];
Print[N[b = 
    Im[g] - (1/Pi)^(d + 1)*
      NIntegrate[
       Simplify[Sin[Pi*x]*D[f[x], {x, d + 1}]], {x, 1, Infinity}, 
       WorkingPrecision -> prec*(105/100), 
       PrecisionGoal -> prec*(105/100)], prec]];
Print[SessionTime[] - T0, " seconds"]; Print[
 N[Sqrt[a^2 + b^2], prec]]; Print[DateString[]];

Here is new a graph of those speed records with a trendline: enter image description here

Edit:

On Tue 26 May 2015 06:21:00, I started a 30,000 digit computation using the following code.

Does any one else want to try to break the record?

 $MaxExtraPrecision = 100; d = 43; f[x_] = x^(1/x); ClearAll[a, b, h];
h[n_] := Sum[
  StirlingS1[n, k]*Sum[(-j)^(k - j)*Binomial[k, j], {j, 0, k}], {k, 1,
    n}]; h[0] = 1; g = 
 2 I/Pi - Sum[-I^(n + 1) h[n]/Pi^(n + 1), {n, 1, d}]; Print[
 DateString[]];
Print[T0 = SessionTime[]]; prec = 30000;
Print[N[a = -Re[g] + (1/Pi)^(d + 1)*
      NIntegrate[
       Simplify[Cos[Pi*x]*D[f[x], {x, d + 1}]], {x, 1, Infinity}, 
       WorkingPrecision -> prec*(105/100), 
       PrecisionGoal -> prec*(105/100)], prec]];
Print[SessionTime[] - T0, " seconds"];
Print[N[b = 
    Im[g] - (1/Pi)^(d + 1)*
      NIntegrate[
       Simplify[Sin[Pi*x]*D[f[x], {x, d + 1}]], {x, 1, Infinity}, 
       WorkingPrecision -> prec*(105/100), 
       PrecisionGoal -> prec*(105/100)], prec]];
Print[SessionTime[] - T0, " seconds"]; Print[
 N[Sqrt[a^2 + b^2], prec]]; Print[DateString[]];

Edit: My first full 30,0000 run finished on Sun 31 May 2015 00:45:09.

Time span: {"4.767 days", "114.4 hours", "6864 minutes", "411849 seconds"} See attached MKB30K2.nb worksheet.

Here is an updated speed record plot, with trendline. (I think the 30,000 digit run can be done faster.)

enter image description here

POSTED BY: Marvin Ray Burns
Answer
2 years ago

I think I came up with a rough program that computes any "prec" digits of the integral analog of the MRB constant. It chooses, d, the best (or close to the best) order of derivative to use in Mathar's algorithm mentioned in a previous post (formula (12) at http://arxiv.org/pdf/0912.3844v3.pdf ), Then uses the appropriate code that integrates the integral analog of the constant. It shows the real and imaginary parts as postive real constants and the value the integral, and gives some timings. It could use a lot of cleanup! I hope someone can help me test it with varying values of prec. Please reply your intentions to use it and results.If no one else can clean it up I will after I tested it more.

prec = 2000; d = Ceiling[0.264086 + 0.00143657 prec]; If[
 Mod[d, 4] == 0, f[x_] = x^(1/x); ClearAll[a, b, h];
 a[n_] := 
  Sum[StirlingS1[n, k]*
    Sum[(-j)^(k - j)*Binomial[k, j], {j, 0, k}], {k, 1, n}]; a[0] = 1;
 g = 2 I/Pi - Sum[-I^(n + 1) a[n]/Pi^(n + 1), {n, 1, d}];
 Print[DateString[]];
 Print[T0 = SessionTime[]];
 Print[N[a = -Re[g] - (1/Pi)^(d + 1)*
      NIntegrate[
       Simplify[Sin[Pi*x]*D[f[x], {x, d + 1}]], {x, 1, Infinity}, 
       WorkingPrecision -> prec*(105/100), 
       PrecisionGoal -> prec*(105/100)], prec]];
 Print[SessionTime[] - T0, " seconds"];
 Print[N[b = 
    Im[g] - (1/Pi)^(d + 1)*
      NIntegrate[
       Simplify[Cos[Pi*x]*D[f[x], {x, d + 1}]], {x, 1, Infinity}, 
       WorkingPrecision -> prec*(105/100), 
       PrecisionGoal -> prec*(105/100)], prec]];
 Print[SessionTime[] - T0, " seconds"];
 Print[N[Sqrt[a^2 + b^2], prec]]; Print[DateString[]];, 
 If[Mod[d, 4] == 1, f[x_] = x^(1/x); ClearAll[a, b, h];
  h[n_] := 
   Sum[StirlingS1[n, k]*
     Sum[(-j)^(k - j)*Binomial[k, j], {j, 0, k}], {k, 1, n}];
  h[0] = 1; g = 2 I/Pi - Sum[-I^(n + 1) h[n]/Pi^(n + 1), {n, 1, d}];
  Print[DateString[]];
  Print[T0 = SessionTime[]];
  Print[N[
    a = -Re[g] - (1/Pi)^(d + 1)*
       NIntegrate[
        Simplify[Cos[Pi*x]*D[f[x], {x, d + 1}]], {x, 1, Infinity}, 
        WorkingPrecision -> prec*(105/100), 
        PrecisionGoal -> prec*(105/100)], prec]];
  Print[SessionTime[] - T0, " seconds"];
  Print[N[
    b = Im[g] + (1/Pi)^(d + 1)*
       NIntegrate[
        Simplify[Sin[Pi*x]*D[f[x], {x, d + 1}]], {x, 1, Infinity}, 
        WorkingPrecision -> prec*(105/100), 
        PrecisionGoal -> prec*(105/100)], prec]];
  Print[SessionTime[] - T0, " seconds"];
  Print[N[Sqrt[a^2 + b^2], prec]]; Print[DateString[]];, 
  If[Mod[d, 4] == 2, f[x_] = x^(1/x); ClearAll[a, b, h];
   a[n_] := 
    Sum[StirlingS1[n, k]*
      Sum[(-j)^(k - j)*Binomial[k, j], {j, 0, k}], {k, 1, n}];
   a[0] = 1; g = 2 I/Pi - Sum[-I^(n + 1) a[n]/Pi^(n + 1), {n, 1, d}];
   Print[DateString[]];
   Print[T0 = SessionTime[]];
   Print[N[
     a = -Re[g] - (1/Pi)^(d)*
        NIntegrate[
         Simplify[Cos[Pi*x]*D[f[x], {x, d}]], {x, 1, Infinity}, 
         WorkingPrecision -> prec*(105/100), 
         PrecisionGoal -> prec*(105/100)], prec]];
   Print[SessionTime[] - T0, " seconds"];
   Print[N[
     b = Im[g] + (1/Pi)^(d + 1)*
        NIntegrate[
         Simplify[Cos[Pi*x]*D[f[x], {x, d + 1}]], {x, 1, Infinity}, 
         WorkingPrecision -> prec*(105/100), 
         PrecisionGoal -> prec*(105/100)], prec]];
   Print[SessionTime[] - T0, " seconds"];
   Print[N[Sqrt[a^2 + b^2], prec]]; Print[DateString[]];, 
   If[Mod[d, 4] == 3, f[x_] = x^(1/x); ClearAll[a, b, h];
    h[n_] := 
     Sum[StirlingS1[n, k]*
       Sum[(-j)^(k - j)*Binomial[k, j], {j, 0, k}], {k, 1, n}];
    h[0] = 1; g = 2 I/Pi - Sum[-I^(n + 1) h[n]/Pi^(n + 1), {n, 1, d}];
    Print[DateString[]];
    Print[T0 = SessionTime[]];
    Print[
     N[a = -Re[g] + (1/Pi)^(d + 1)*
         NIntegrate[
          Simplify[Cos[Pi*x]*D[f[x], {x, d + 1}]], {x, 1, Infinity}, 
          WorkingPrecision -> prec*(105/100), 
          PrecisionGoal -> prec*(105/100)], prec]];
    Print[SessionTime[] - T0, " seconds"];
    Print[
     N[b = Im[g] - (1/Pi)^(d + 1)*
         NIntegrate[
          Simplify[Sin[Pi*x]*D[f[x], {x, d + 1}]], {x, 1, Infinity}, 
          WorkingPrecision -> prec*(105/100), 
          PrecisionGoal -> prec*(105/100)], prec]];
    Print[SessionTime[] - T0, " seconds"];
    Print[N[Sqrt[a^2 + b^2], prec]]; Print[DateString[]];]]]]

Here are some of my best timings to compare with the program's results:

digits  seconds
1000    38.8650545
2000    437.4906125
3000    889.473875
4000    1586.000714
5000    2802.591704
6000    4569.41586
7000    6891.057587
8000    9659.318566
9000    13491.43967
10000   17355
11000   
12000   
13000   
14000   
15000   53385.02323
16000   
17000   
18000   
19000   
20000   123876.4331
21000   
22000   
23000   
24000   
25000   229130.3088
26000   
27000   
28000   
29000   
30000   411848.6322

Edit: On Fri 5 Jun 2015 20:41:45 I started a 35,000 digit computation with the above "automated" program.

Edit: The 35,000 digit computation should be done by 10:24:38 am EDT | Sunday, June 14, 2015. In the above "automated" program I forgot to adjust the MaxExtraPrecision, but that shouldn't affect the accuracy in that program. It already computed the real part of the integral to 35,000 digits and the first 30,000 of those are the same as the real part of my previously mentioned 30,000 digit calculation. I will keep you posted.

Edit: The 35,000 digit computation finished on Sun 14 Jun 2015 06:52:29. It is attached as 35KMKB.nb. The first 30,000 digits of those are the same as the ones of my previously mentioned 30,000 digit calculation. (That shows the computation didn't take any "wild" turns because of the lack of MaxExtraPrecision.) Further it is a good check of the 30,000 digit run, as all of the bigger computations are of the smaller, because they all are calculated with distinct formulas using different orders of the derivative of x^(1/x).

727844 seconds

Attachments:
POSTED BY: Marvin Ray Burns
Answer
2 years ago

Here is an extensive record of records of computing the integral analog of the MRB constant:

enter image description here

enter image description here

Here is a graph of those records. (The progression of computed digits is so extreme, it is almost unbelievable!) enter image description here

POSTED BY: Marvin Ray Burns
Answer
2 years ago

For 2000 digits Mathematica 10. 2.0 shows some remarkable improvement over 10.1.2 with the above "automated program" for computing the digits of the integral analog of the MRB constant.

I will post some speed records that are strictly what the program produces in V 10.2.0, below, no picking and choosing of the methods by a human being. Some results will naturally be slower than my previously mentioned speed records, because I tried so very many methods and recorded only the fastest results.

digits          seconds

2000    256.3853590 
3000    794.4361122
4000       1633.5822870
5000        2858.9390025
10000      17678.7446323 
20000      121431.1895170
40000       I got error msg

to be continued

I have to change the program for 40,000 digits! I'll post the new program when I get 40,000 to work.

POSTED BY: Marvin Ray Burns
Answer
2 years ago

As of Wed 29 Jul 2015 11:40:10, one of my computers was happily and busily churning away at 40,000 digits of the integral analog of the MRB constant, using the following formula.

Edit: Mathematica crashed at 11:07 PM 8/4/2016

(I used MKB as a symbol for the integral analog because it is called the MKB constant. You can find the name MKB constant at http://www.ebyte.it/library/educards/constants/MathConstants.html .) If you can weed through my code, at the bottom of this reply, you might want to check the formula for the placement of pluses, minuses and imaginary units!!! A little hint when checking if the formula matches the code, d is 80 so Mod[d,4] =0.

f[x_] = x^(1/x) : a[n_] := 
  Sum[StirlingS1[n, k]*
    Sum[(-j)^(k - j)*Binomial[k, j], {j, 0, k}], {k, 1, n}]; a[0] = 1;
g = 2 I/Pi - Sum[-I^(n + 1) a[n]/Pi^(n + 1), {n, 1, 80}]
MKB = -g + (I/Pi)^81*
   Integrate[f[x]*D[f[x], {x, 81}], {x, 1, Infinity}]

In traditional form that looks like

enter image description here

Here is the code,cleaned up a little: This is the version from Aug 6, 2015 452 pm; for the first time the imaginary part is signed and shown to be multiplied by the imaginary unit!

Block[{$MaxExtraPrecision = 200}, prec = 4000; f[x_] = x^(1/x);
 ClearAll[a, b, h];
 Print[DateString[]];
 Print[T0 = SessionTime[]];

 If[prec > 35000, d = Ceiling[0.002 prec], 
  d = Ceiling[0.264086 + 0.00143657 prec]];

 h[n_] := 
  Sum[StirlingS1[n, k]*
    Sum[(-j)^(k - j)*Binomial[k, j], {j, 0, k}], {k, 1, n}];

 h[0] = 1;
 g = 2 I/Pi - Sum[-I^(n + 1) h[n]/Pi^(n + 1), {n, 1, d}];

 sinplus1 := 
  NIntegrate[
   Simplify[Sin[Pi*x]*D[f[x], {x, d + 1}]], {x, 1, Infinity}, 
   WorkingPrecision -> prec*(105/100), 
   PrecisionGoal -> prec*(105/100)];

 cosplus1 := 
  NIntegrate[
   Simplify[Cos[Pi*x]*D[f[x], {x, d + 1}]], {x, 1, Infinity}, 
   WorkingPrecision -> prec*(105/100), 
   PrecisionGoal -> prec*(105/100)];

 middle := Print[SessionTime[] - T0, " seconds"];

 end := Module[{}, Print[SessionTime[] - T0, " seconds"];
   Print[c = Abs[a + b]]; Print[DateString[]]];


 If[Mod[d, 4] == 0, 
  Print[N[a = -Re[g] - (1/Pi)^(d + 1)*sinplus1, prec]];
  middle;
  Print[N[b = -I (Im[g] - (1/Pi)^(d + 1)*cosplus1), prec]];
  end];


 If[Mod[d, 4] == 1, 
  Print[N[a = -Re[g] - (1/Pi)^(d + 1)*cosplus1, prec]];
  middle;
  Print[N[b = -I (Im[g] + (1/Pi)^(d + 1)*sinplus1), prec]]; end];

 If[Mod[d, 4] == 2, 
  Print[N[a = -Re[g] + (1/Pi)^(d + 1)*sinplus1, prec]];
  middle;
  Print[N[b = -I (Im[g] + (1/Pi)^(d + 1)*cosplus1), prec]];
  end];

 If[Mod[d, 4] == 3, 
  Print[N[a = -Re[g] + (1/Pi)^(d + 1)*cosplus1, prec]];
  middle;
  Print[N[b = -I (Im[g] - (1/Pi)^(d + 1)*sinplus1), prec]];
  end];]

Come back to see if I decided whether to try the 40K run again.

EDIT: It looks like I've only got one more test for the program (if it passes) before I retry the 40,000 digit calculation!

EDIT:On Thu 6 Aug 2015 17:23:18, I restarted the 40K run with Windows 10.

EDIT: My first thought was the program took up too much RAM, apparently over 115 GB! ( I have 64GB installed and a 51 GB paging file; nevertheless, Windows 10 closed the Mathematica kernel to keep from the computer from loosing data. Can someone else try the 40K run on their computer? It should take 2 weeks on a fast one. Please let me know if you try it and let me know the results, so I will know I don't have a problem with my computer., If two weeks is too great of a commitment, can you try taking note on the RAM used for two progressively larger runs, like 20K and 30 K? I will do the same, and we can compare notes. Thank You!

EDIT: I've been monitoring memory usage for smaller runs and found the program only uses minimal memory! This makes the action of Windows 10 (closing Mathematica kernel to avoid data loss) all the more a mystery! Could the 40K run really use up all of that RAM?

I know there are quite a few of you viewing this post; however, is anyone out there working on these calculations?.

POSTED BY: Marvin Ray Burns
Answer
2 years ago

I think I've reached an impasse in calculating more digits of the integral analog of the MRB constant. So I believe my best bet in setting a new record is to compute 4,000,000 digits of the MRB constant proper!

Here is the data of my previous records on my big computer, the six core Intel(R) Core(TM) i7-3930K CPU @ 3.20 GHz 3.20 GHz with 64 GB of RAM.

(Digits are rounded down, slightly.)

1.2 million in 120,360 seconds

2 million 1,870,579 second

3 million 4,15,7400 seconds

Remarkably they are very linear; so using the Fit command as shown next,

Fit[data, {1, x}, x]

I get

y, in seconds, = -2.5894022213114775*10^6 + 2.2446041393442626 x.

Substituting 4,000,000 for x gives 6,389,014;seconds or almost 2 and a half months. (I would round that to 3 months to be safe!)

Before I start, I will listen to your suggestions, or encouragement; if you have any, (I do get lonely when I don't get any feedback.)

EDIT:

At 10:00 PM Friday Sept 04, 2015 I started the 4,000,000 run of the MRB constant proper. Does anyone want to work on it with me?

Here is the code I am using:

(*Fastest (at MRB's end) as of 24 dEC 2014.*)
Block[{$MaxExtraPrecision = 50}, 
 prec = 4000000;(*Number of required decimals.*)ClearSystemCache[];
 T0 = SessionTime[];
 expM[pre_] := 
  Module[{a, d, s, k, bb, c, n, end, iprec, xvals, x, pc, cores = 16, 
    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["end ", end];
   Print[end*chunksize]; d = Cos[n ArcCos[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;
        h = Log[ll]/ll; x = N[Exp[h], 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" -> 4]];
    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" -> 2];
    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];]; 
 Print["Difference from 3014991 known digits is ", MRBtest2 - m3M];
 MRBtest2]

3014991 digits are saved as m3M.

I will post several updates as edits here as I go.

EDIT: As of Sat, Sept 5, 2015, the program is indicating that the run will take 3.3 months.

EDIT: Having to run the program for over 3 months, there is a chance that I would suffer a power outage longer than the 15 minutes or so that my battery backup will run my computer! Of course I'll keep you posted if that happens. (I've tried calculations that took longer than that before!)

EDIT I didn't know that upgrading to Windows 10 caused my Windows update option to read, automatically install updates, so my computer restarted. I started the 4,000,000 digit run again at 8:34 PM, Thursday Sept 10, 2015. This time I made sure the deffer updates option was chosen!!

EDIT As of Sept 19,2015, the program says it will take 94 more days to complete the 4 million digit run.

EDIT As of Sept 23,2015, 9:20 PM, the program says it will take 90.9563 more days to complete the 4 million digit run.

EDIT As of Sept 27,2015, 9:20 PM, the program says it will take 86.9283 more days to complete the 4 million digit run

POSTED BY: Marvin Ray Burns
Answer
2 years ago

Windows 10 restarted my machine again.The Windows updates were installed again, arrggghhhh!

I don't know if I am going to retry it again real soon. (It really ties up my best computer!)

This might give other people a better chance to break my record.

POSTED BY: Marvin Ray Burns
Answer
2 years ago

I am waiting for an improvement in Mathematica to try to break any new records, especially because the goal of 4,000,000 digits of the MRB constant requires, in the neighborhood of, 3 months worth of computing time with the best version so far! Also, I thought I would take advantage of this opportunity to introduce this discussion to any new members, that they too may learn or be encouraged by its voluminous replies, as manifold experienced members have -- and a few have added to it-- contributing to its success. Thank you all!

In this post we have listed record calculations of the MRB constant (MRB) and its integral analog (MKB); found new formulas for computing each (many, many formulas for computing MKB with each more complicated and powerful that the previous likeenter image description here

); looked at records using only special formulas for MRB; and looked at a few related constants, like the infinite power tower of MRB. I haven't decided whether to start a new post when I try my next record, or not; if i decide to start a new one I will post a link to it here. Please feel free to continue to use this post for your contributions or questions!

POSTED BY: Marvin Ray Burns
Answer
1 year ago

I've been working on Crandall's code for many digits of the MRB constant. On my big computer Windows 10 is just a little slower than Windows 7 was. This new code runs on Windows 10 at the same speed the my last mention in this blog program did in Windows 7.

Here is what it now looks like in traditional form: enter image description here enter image description here enter image description here

Here it is in input form:

Quiet[prec = 30000; ClearSystemCache[]; 
    T0 = SessionTime[]; expM[pre_] := 
      Module[{a, d, s, k, bb, c, n, end, iprec, xvals, x, 
          pc, cores = 16, 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["end ", end]; Print[end*chunksize]; 
         d = (1/2)*(3 + 2*Sqrt[2])^n*(1 + (3 + 2*Sqrt[2])^
                  (-2*n)); {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[ll^(1/ll), iprec]; pc = iprec; 
                     While[pc < pr, pc = Min[3*pc, pr]; x = 
                          SetPrecision[x, pc]; y = x^ll - ll; x = 
                          x - (2*x*y)/(y + ll*(2*ll + y))]; x, 
                    {l, 0, tsize - 1}], {j, 0, cores - 1}, 
                  Method -> "EvaluationsPerKernel" -> 4]]; 
            ctab = ParallelTable[Table[c = b - c; 
                   ll = start + l - 2; b *= (2*(ll^2 - n^2))/
                       (1 + 3*ll + 2*ll^2); c, {l, chunksize}], 
                Method -> "EvaluationsPerKernel" -> 2]; 
            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]

Also, here is my new fastest program for calculating digits of the MRB constant via the eta formula,

enter image description here This program is nearly twice as fast as the one in I gave "best eta timings" for, above..

prec = 1500;
to = SessionTime[]; 
etaMM[m_, pr_] := 
 Module[{a, d, s, k, b, c}, a[j_] := N[(Log[j + 1]/(j + 1))^m, pr];
  n = Floor[132 /100 pr];
  d = N[(1/2)*(3 + 2*Sqrt[2])^n*(1 + (3 + 2*Sqrt[2])^(-2*n)), prec];
  {b, c, s} = {-1, -d, 0};
  Do[c = b - c;
   s = s + c a[k];
   b = (k + n) (k - n) b/((k + 1) (k + 1/2)), {k, 0, n - 1}];
  Return[N[s/d, pr] (-1)^m]];
eta1 = N[EulerGamma Log[2] - Log[2]^2/2, prec];
MRBtest = 
  eta1 - ParallelSum[(Cos[Pi m]) etaMM[m, prec]/
      N[Gamma[m + 1], prec], {m, 2, Floor[.40 prec]}];
Print[MRBtest]; SessionTime[] - to

gives

(*    During evaluation of In[188]:= 0.18785964246206712024851793405427323005590309490013878617200468408947723156466021370329665443310749690384234585625801906123137009475922663043892934889618412083733662608161360273812637937343528321255276396217148932170207628206217151671540841268044836354167199851976802527598938993914457983505561350964852107120784442309586812949768852694956420425558648367044104252795247106066609263397483410311578167864166891546003422225883800254553968929471142122189105098328712277308020036445215390536395055332203470627551159812828039510219264914673176293516190659816018664245824950697203381992958420935515162514399357600764593291281451709082424915883204169066409334435914806705564692806787007028115009380606938139385953360657987405562062348704329360737819564603104763950664893061360645528067515193508280837376719296866398103094949637496277383049846324563479311575300289212523291816195626973697074865765476071178017195787368300965902260668753656305516567361288150201438756136686552210674305370591039735756191489093690777983203551193362404637253494105428363699717024418551654837279358822008134480961058802030647819619596953756287834812334976385863010140727252923014723333362509185840248037040488819676767601198581116791693527968520441600270861372286889451015102919988536905728659287086875425492533794395347589703563313440382638887986656195980733514739902565778133172261076127975852722742777308985774922305970962572562718836755752978879253616876739403543214513627725492293131262764357321446216187786377154205423128223

    Out[192]= 105.7159526*)
POSTED BY: Marvin Ray Burns
Answer
1 year ago

V 11 has a significant improvement in my new most recently mentioned fastest program for calculating digits of the MRB constant via the eta formula, Here are some timings:

Digits           seconds

1500                42.6386632

2000             127.3101969

3000             530.4442911

4000           1860.1966540

5000           3875.6978162

6000           8596.9347275



 10,000        53667.6315476

From an previous message that starts with "How about computing the MRB constant from Crandall's eta derivative formulas?" here are my first two sets of records to compare with the just mentioned ones. You can see that I increased time efficiency by 10 to 29 to even 72 fold for select computations! In the tests used in that "previous message," 4000 or more digit computations produced a seemingly indefinitely long hang-on.

enter image description here

Digits      Seconds

500        36.831836
1000       717.308198
1500       2989.759165
2000       3752.354453


Digits        Seconds
 500          9.874863
 1000        62.587601
 1500        219.41540
 2000       1008.842867
 2500       2659.208646
 3000       5552.902395
 3500       10233.821601

Comparing first of the just mentioned 2000 digit computations with the "significant improvement" one we get the following.

3752/127 ~=29.

And from the slowest to the fastest 1500 digit run we get

2989/42 ~=72,

POSTED BY: Marvin Ray Burns
Answer
11 months ago

Reply starts with UPDATE on 8/13/2017 below.

Search 8/13/2017 below for where this reply was placed.

POSTED BY: Marvin Ray Burns
Answer
7 days ago

On a recent post, I derived a bunch of approximations of the MRB constant.

It was MRB constant approximations using TranscendentalRecognize -- http://community.wolfram.com/groups/-/m/t/819837?p_p_auth=m7XRds7k.

Some of the last ones consist of arbitrary approximations to any n digits,

while still using small coefficients, compared with the thousand of digits of the target.

I wonder if anyone could show me how to extract the coefficient like

terms which come before the summations, in a list like this:

enter image description here

So someone could find a closed form consisting of those coefficient like

terms

POSTED BY: Marvin Ray Burns
Answer
1 year ago

Repeated on top of first message

When it comes to mine and a few more educated people's passion to calculate many digits and the dislike of a few more educated people; it is all a matter telling us that the human mind is multi faceted in giving passion, to person a, for one task and to person b for another task!

POSTED BY: Marvin Ray Burns
Answer
1 year ago

My plan in the previous was to find a numeric list to dot product(.) a set of sums that equal the MRB constant,

list or at least gave an improvable result. In the post

( MRB constant approximations using TranscendentalRecognize ),

I have made several numeric list and sums that are inner products,

but I haven't found any good pattern to the numeric lists!

Below is a sample of what I want to investigating with just a sample value for a, b , c..

So from

 m = NSum[(-1)^n*(n^(1/n) - 1), {n, 1, Infinity},Method -> "AlternatingSigns", WorkingPrecision -> 1000]
  N[m - (-129858773922357615372945307143544254 - 
 332618118135196201861563173048187520*3^(1/5) + 
 360989141074787715535168098417020609*2^(2/5)*3^(1/5))/
261711912538111957032871243762602971, 200]

 Out[116]= 0.*10^-99

you could have

In[112]:= 

a = {-129858773922357615372945307143544254, 
   332618118135196201861563173048187520, 
   360989141074787715535168098417020609};


In[113]:= b = {1, -3^(1/5), 2^(2/5) 3^(1/5)};

In[114]:= c = 261711912538111957032871243762602971;

In[115]:= N[m - (a.b/c), 200]

Out[115]= 0.*10^-99

If only there was a formula for the elements of a,b, and c; up to the precision given by TranscendentalRecognize in a polynomial like from,

you could derive a,b,c without TranscendentalRecognize getting arbitrary digits of m of the simple formula of a.b/c .

POSTED BY: Marvin Ray Burns
Answer
1 year ago

I just reasoned a great reason to compute the MRB constant!

As we saw above at https://web.archive.org/web/20130430193005/http://www.perfscipress.com/papers/UniversalTOC25.pdf , Richard Crandall wrote about Unified algorithms for polylogarithm, L-series, and zeta variants. Since this paper was cited 8 times in its first year (how much longer he lived), I believe this paper is suitable for an introduction to most of what is said in the subject mentioned above. In chapter 7 i.e. "Key fundamental constants." I searched for the meaning of Key fundamental constants, on google and found it referred to different constants, all depending on the subject that uses the constant. So a whole section of that chapter was called The MRB constant. So I wonder does considering the MRB constant to be a "Key fundamental constants." make to calculating of its digits any more interesting?

Even if you don't want to use all of the computer power to break a record here, If anyone can find a faster algorithm that I could use, that would bee terrific!

POSTED BY: Marvin Ray Burns
Answer
1 year ago

I think this article might be starting a second life for

enter image description here,

since the efficient analysis of equation 44 is missing:

enter image description here

http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.695.5959&rep=rep1&type=pdf <- New occasionally working reference, so subject not dead! This reference, on my site, is up a little more: http://www.marvinrayburns.com/UniversalTOC25.pdf .

We touched formula 44 but didn't dig out any of the good stuff that Crandall believed to be in it!

Here B is the MRB constant. Plus Crandall acknowledged that there should be a negative sign before first summation.

enter image description here

POSTED BY: Marvin Ray Burns
Answer
1 year ago

V. 11 is about 1.25 times faster with my newest program for calculating MKB, (the integral analog of the MRB sum). V 10. 4 calculated 20,000 digits in 121431.1895170 seconds and V 11 did it in 96979.6545388 seconds. I've got a little more testing to do, (about 1 day's worth), then I'll try 40,000 digits again, which should take about 12 days. I will post all my updates here, so you might want to save this message as a favorite so you won't loose it.

Update 1

The 40 K automatically started against my wishes on Thu 11 Aug 2016 15:42:08, (due to my pasting two codes at once). I'll keep you informed, how it goes.

Update 2

Windows 10 is pushing an update. Wednesday is the latest it will let me restart. I will restart now with all the updates I can get, Then deffer further ones and hopefully get 12 restart free days to do my 40K.

Update 3

I ran all the updates I could find, differed further ones and restarted 40K on Sun 14 Aug 2016 10:32:40.

Update 4

Widows 10 stopped the calculation! AGAIN! Can anyone else try it and see if you get anywhere? Here is my latest code:

(*Other program:For large calculations.Tested for 1000-35000 digits-- \
see post at \
http://community.wolfram.com/groups/-/m/t/366628?p_p_auth=KA7y1gD4 \
and search for "analog" to find pertinent replies.Designed to include \
40000 digits.A157852 is saved as c,the real part as a and the \
imaginary part as b.*)Block[{$MaxExtraPrecision = 200}, 
 prec = 40000(*Replace 40000 with number of desired digits.40000 \
digits should take two weeks on a 3.5 GH Pentium processor.*); 
 f[x_] = x^(1/x);
 ClearAll[a, b, h];
 Print[DateString[]];
 Print[T0 = SessionTime[]];
 If[prec > 35000, d = Ceiling[0.002 prec], 
  d = Ceiling[0.264086 + 0.00143657 prec]];
 h[n_] := 
  Sum[StirlingS1[n, k]*
    Sum[(-j)^(k - j)*Binomial[k, j], {j, 0, k}], {k, 1, n}];
 h[0] = 1;
 g = 2 I/Pi - Sum[-I^(n + 1) h[n]/Pi^(n + 1), {n, 1, d}];
 sinplus1 := 
  NIntegrate[
   Simplify[Sin[Pi*x]*D[f[x], {x, d + 1}]], {x, 1, Infinity}, 
   WorkingPrecision -> prec*(105/100), 
   PrecisionGoal -> prec*(105/100)];
 cosplus1 := 
  NIntegrate[
   Simplify[Cos[Pi*x]*D[f[x], {x, d + 1}]], {x, 1, Infinity}, 
   WorkingPrecision -> prec*(105/100), 
   PrecisionGoal -> prec*(105/100)];
 middle := Print[SessionTime[] - T0, " seconds"];
 end := Module[{}, Print[SessionTime[] - T0, " seconds"];
   Print[c = Abs[a + b]]; Print[DateString[]]];
 If[Mod[d, 4] == 0, 
  Print[N[a = -Re[g] - (1/Pi)^(d + 1)*sinplus1, prec]];
  middle;
  Print[N[b = -I (Im[g] - (1/Pi)^(d + 1)*cosplus1), prec]];
  end];
 If[Mod[d, 4] == 1, 
  Print[N[a = -Re[g] - (1/Pi)^(d + 1)*cosplus1, prec]];
  middle;
  Print[N[b = -I (Im[g] + (1/Pi)^(d + 1)*sinplus1), prec]]; end];
 If[Mod[d, 4] == 2, 
  Print[N[a = -Re[g] + (1/Pi)^(d + 1)*sinplus1, prec]];
  middle;
  Print[N[b = -I (Im[g] + (1/Pi)^(d + 1)*cosplus1), prec]];
  end];
 If[Mod[d, 4] == 3, 
  Print[N[a = -Re[g] + (1/Pi)^(d + 1)*cosplus1, prec]];
  middle;
  Print[N[b = -I (Im[g] - (1/Pi)^(d + 1)*sinplus1), prec]];
  end];] (*Marvin Ray Burns,Aug 06 2015*)
POSTED BY: Marvin Ray Burns
Answer
1 year ago

To try to get windows 10 from closing Mathematica during the computation I tried the instructions found at https://www.autoitscript.com/forum/topic/177749-stopping-windows-10-from-auto-closing-programs-to-free-up-ram/ . I will record progress in this spot as I did before.

UPDATE I followed the memory usage on my computer and it did use around 64 GB of RAM. And then Windows closed down the Mathematica kernel. I assume that If I can ever afford to maximize my RAM to its 128GB limit the computation will be successful!

Anyone have better luck?

POSTED BY: Marvin Ray Burns
Answer
9 months ago

I have two versions of MRB constant code that I presently use, a short one (non-parallel) and a long one. Both are given below.

enter image description here

Short

(*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]

Long

 (**Fastest (at MRB's end) as of 25 Jul 2014*.*)DateString[]
 prec = 7000;
 (**Number of required decimals.*.*)ClearSystemCache[];
 T0 = SessionTime[];
 expM[pre_] := 
   Module[{a, d, s, k, bb, c, n, end, iprec, xvals, x, pc, cores = 12, 
     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["end ", end];
    Print[end*chunksize]; 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, 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" -> 4]];
     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" -> 2];
     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];]; DateString[]
 Print[MRBtest2]
 MRBtest2 - MRBtest3
POSTED BY: Marvin Ray Burns
Answer
6 months ago

Try to beat these MRB constant records!

POSTED BY: Marvin Ray Burns
Answer
1 month ago

As mentioned a few times above Richard Crandall made the following observations about the MRB constant: enter image description here

His last comment is that we might find some algorithm advantages by focusing on the fixed argument of 0 for the eta derivatives. He gives Chapter 3 as a possible source to find the algorithm advantages. (Here is another link to his paper: http://marvinrayburns.com/UniversalTOC25.pdf .)

Here is what he was looking for.

c[0] = 0; 
c[1] = Log[2]; Table[
c[n] = -2*Derivative[n - 1][Zeta][0], {n, 2, 22}]; 
l = CoefficientList[Normal[Series[Exp[(-x)*E^x], {x, 0, 21}]], 
x]; N[-Sum[
l[[n + 1]]*((-1)^
n*(Sum[(-1)^x*(c[x + 1]*Binomial[n, x]*c[1]^(n - x)), {x, 1, 
n - 1}] + c[1]^n + 
c[0]) + 1/2 c[n + 1]), {n, 1, 17}]]

Mathematica still has trouble with higher zeta derivatives of 0, however. I'm going to look for a procedure for computing them better.

[Edit]

I computed the zeta derivatives of 0 by

   zeta0[p_] := (-1)^p Sum[StieltjesGamma[p + k]/k!, {k, 0, Infinity}] - 
     p!

and got the same results Mathematica got for the higher zeta derivatives of 0.

I also checked the CoefficientList of

CoefficientList[Normal[Series[Exp[(-x)*E^x], {x, 0, 21}]]

by the following, and haven't found where Mathematica is having trouble!

c[0] = 0; c[1] = Log[2]; Table[c[n] = -2*Derivative[n - 1][Zeta][0], 
{n, 2, 
18}]; N[-Sum[(Sum[(-1)^d*d^(n - d)*Binomial[n, d], {d, 1, n}]*
((-1)^
n*(Sum[(-1)^x*(c[x + 1]*Binomial[n, x]*c[1]^(n - x)), 
{x, 1, n - 1}] + c[1]^n + c[0]) + (1/2)*
c[n + 1]))/n!, {n, 1, 17}]]

I even checked the binomials by their formula. Also the StieltjesGamma and Derivative[n ][Zeta][0]] results appear smooth for discrete operations:

 Table[N[Derivative[n ][Zeta][0]], {n, 10, 25}]

and

 N[Table[StieltjesGamma[x], {x, 1, 30}]]

.

[Edit]

I found the problem as to why Mathematica doesn't increase its accuracy in the above formula after 17 iterations: Mathematica will give a numeric solution to many zeta derivatives of 0, only to machine precision!

For machine precision:

  Table[N[Derivative[n ][Zeta][0]], {n, 1, 30}]

gives the machine sized numeric to many derivatives:

 {-0.918939, -2.00636, -6.00471, -23.9971, -120., -720.001, \
-5040., -40320., -362880., -3.6288*10^6, -3.99168*10^7, \
-4.79002*10^8, -6.22702*10^9, -8.71783*10^10, -1.30767*10^12, \
-2.09228*10^13, -3.55687*10^14, -6.40237*10^15, -1.21645*10^17, \
-2.4329*10^18, -5.10909*10^19, -1.124*10^21, -2.5852*10^22, \
-6.20448*10^23, -1.55112*10^25, -4.03291*10^26, -1.08889*10^28, \
-3.04888*10^29, -8.84176*10^30, -2.65253*10^32}.

Then

   Table[N[Derivative[n][Zeta][0], 10], {n, 1, 30}]

gives the high precision numeric for only the first 3 derivatives:

{-0.91893853320467274178032973640561763986`10., \
-2.0063564559085848512101000267299604382`10., 
   -6.0047111668622544478`10., Derivative[4][Zeta][0], 
 Derivative[5][Zeta][0], Derivative[6][Zeta][0], 
   Derivative[7][Zeta][0], Derivative[8][Zeta][0], 
 Derivative[9][Zeta][0], Derivative[10][Zeta][0], 
   Derivative[11][Zeta][0], Derivative[12][Zeta][0], 
 Derivative[13][Zeta][0], Derivative[14][Zeta][0], 
   Derivative[15][Zeta][0], Derivative[16][Zeta][0], 
 Derivative[17][Zeta][0], Derivative[18][Zeta][0], 
   Derivative[19][Zeta][0], Derivative[20][Zeta][0], 
 Derivative[21][Zeta][0], Derivative[22][Zeta][0], 
   Derivative[23][Zeta][0], Derivative[24][Zeta][0], 
 Derivative[25][Zeta][0], Derivative[26][Zeta][0], 
   Derivative[27][Zeta][0], Derivative[28][Zeta][0], 
 Derivative[29][Zeta][0], Derivative[30][Zeta][0]}.

With some confidence I can say

enter image description here

POSTED BY: Marvin Ray Burns
Answer
1 month ago

I fixed it!!!

I forced Mathematica to give up the needed digits of the zeta derivatives.

My new program does not really compete with the previously mentioned record-breakers. But it is a fairly efficient use of Crandall's formula (44). First it calculates many digits on MRB the old fashion way to use as a check. So here is my new program:

m = NSum[(-1)^n (n^(1/n) - 1), {n, 1, Infinity}, 
   Method -> "AlternatingSigns", WorkingPrecision -> 120];

Block[{$MaxExtraPrecision = 1000}, s = 100(*number of terms*);
 c[1] = Log[2];
 c[n_] := N[1 - 2*Derivative[n - 1][Zeta][0], Floor[3/2 s]] - 1;
 mtest = (-Sum[(Sum[(-1)^d*d^(n - d)*Binomial[n, d], {d, 1, 
          n}]*((-1)^
            n*(Sum[(-1)^x*(c[x + 1]*Binomial[n, x]*c[1]^(n - x)), {x, 
               1, n - 1}] + c[1]^n) + (1/2)*c[n + 1]))/n!, {n, 1, 
      s}])]; Print[s, " terms give ", mtest, ", and the error is ", 
 m - mtest]

100 terms give 0.18785964486341229249700337222, and the error is -2.40134517224848543816*10^-9

Changing s to 200 shows all correct digits:

200 terms give 0.187859642465, and the error is -3.*10^-12

POSTED BY: Marvin Ray Burns
Answer
1 month ago

I noticed that accuracy is not directly related to the number of terms used in the programs of Crandall's formula (44).

I made a slight improvement in my program that computes Crandall's eta derivative of 0 formula (44):

First compute a known value of the MRB constant:

 m = NSum[(-1)^n (n^(1/n) - 1), {n, 1, Infinity}, 
   Method -> "AlternatingSigns", WorkingPrecision -> 120];

followed by

  Block[{$MaxExtraPrecision = 1000}, s = 129(*number of terms*);
   sc[t_] := 
    StieltjesGamma[
     t];(*-Pi/(2 n+2) NIntegrate[Log[1/2+I x]^(n+1)/Cosh[Pi \
  x]^2,{x,-Infinity,Infinity},WorkingPrecision\[Rule]50];*)
   zeta0[p_] := 
    N[(-1)^p Sum[sc[p + k]/k!, {k, 0, Infinity}] - p!, 3/2 s];
   c[1] = Log[2];
   c[n_] := N[1 - 2*zeta0[n - 1], Floor[3/2 s]] - 1;
   mtest = (-Sum[(Sum[(-1)^d*d^(n - d)*Binomial[n, d], {d, 1, 
            n}]*((-1)^
              n*(Sum[(-1)^x*(c[x + 1]*Binomial[n, x]*c[1]^(n - x)), {x, 
                 1, n - 1}] + c[1]^n) + (1/2)*c[n + 1]))/n!, {n, 1, 
        s}])]; Print[s, " terms give ", mtest, ", and the error is ", 
   m - mtest]

producing

129 terms give 0.187859642466672738539355017, and the error is -4.605618290837083*10^-12.

POSTED BY: Marvin Ray Burns
Answer
1 month ago

UPDATE on 8/13/2017

Search 8/13/2017 in the above for this reply's place in Que. I placed it here to make it easy to find on this late of a date.

Here is another great improvement in calculating digits of MRB though this first Crandall eta formula

enter image description here

Here is an example calculating 1500 digits in less than 23 seconds. First calculate known accurate digits:

mTrue = NSum[(-1)^n (n^(1/n) - 1), {n, 1, Infinity}, 
   Method -> "AlternatingSigns", WorkingPrecision -> 2000];

Then use this program:

prec = 1500;
to = SessionTime[];
etaMM[m_, pr_] := 
  Module[{a, s, k, b, c}, a[j_] := (Log[j + 1]/(j + 1))^m;
   {b, c, s} = {-1, -d, 0};
   Do[c = b - c;
    s = s + c a[k];
    b = (k + n) (k - n) b/((k + 1) (k + 1/2)), {k, 0, n - 1}];
   Return[N[s/d, pr] (-1)^m]];
eta1 = N[EulerGamma Log[2] - Log[2]^2/2, prec]; n = 
 Floor[132/100 prec]; d = N[ChebyshevT[n, 3], prec];
MRBtest = 
  eta1 - ParallelSum[(Cos[Pi m]) etaMM[m, prec]/
      N[Gamma[m + 1], prec], {m, 2, Floor[.365 prec]}];
Print[mTrue - MRBtest]; SessionTime[] - to

The larger prec (precission) you want to choose, the smaller number you can replace .365 with in
ParallelSum[(Cos[Pi m]) etaMM[m, prec]/N[Gamma[m + 1], prec], {m, 2, Floor[.365 prec]}. For prec=5000, change .365 to .31 (approximately). This advantage for smaller upper limits being needed in the sum is due to the Gamma (!) in the denominator. On my big computer 5000 digits now gives a timing of around 1800 seconds (1/2 an hour). See message above the previous occurrence of 8/13/2017 for previous records and records before those.

POSTED BY: Marvin Ray Burns
Answer
7 days ago

Group Abstract Group Abstract