Here is a summary of my speed records.
First a full quote from a message summarizing the Burns-Mathar method I derived from the work in this paper.
I made a quicker program for calculating the digits of the MKB constant (also called M1, I{2N} or MKB) in V12.1.0
Module[{$MaxExtraPrecision = 200, sinplus1, cosplus1, middle, end, a,
b, c, d, g, h}, prec = 5000; f[x_] = x^(1/x);
Print[DateString[]];
Print[T0 = SessionTime[]];
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 := Module[{},
NIntegrate[
Simplify[Sin[Pi*x]*D[f[x], {x, d + 1}]], {x, 1, Infinity},
WorkingPrecision -> prec*(105/100),
PrecisionGoal -> prec*(105/100)]];
cosplus1 := Module[{},
NIntegrate[
Simplify[Cos[Pi*x]*D[f[x], {x, d + 1}]], {x, 1, Infinity},
WorkingPrecision -> prec*(105/100),
PrecisionGoal -> prec*(105/100)]];
middle := Module[{}, Print[SessionTime[] - T0, " seconds"]];
end := Module[{}, Print[SessionTime[] - T0, " seconds"];
Print[N[Sqrt[a^2 - b^2], prec]]; 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];]
Whether it will allow me to calculate more digits is a question that
will be answered in a week or two.
Here is a comparison of timings on similar computers.
digits seconds
(Impoved code)
V 10.1.2 V10.3 v11.3 V12.0 V12.1
2000 437 256 67 67 58
3000 889 794 217 211 186
4000 1633 514 492 447
5000 2858 1005 925 854
10000 17678 8327 7748 7470
20000 121431 71000 66177
30000 411848 ? 229560
Seethe following cloud notebook for the results from my improved code.
https://www.wolframcloud.com/obj/bmmmburns/Published/2nd%2040k%20mkb%20prep.nb
Then there are my recent programs from the Abel-Plana formula on V12.0 that computes M2 which is (the MKB constant also called M1, I{2N} or MKB)+C, specifically,
Method ( A, B or C)
in seconds
digits A B C
2000 -> 23 -> 20 -> 18
3000 -> 96 -> 80 -> 36
4000 -> 165 -> 141 -> 64
5000 -> 442 -> 418 -> 386
6000 -> 623 -> 591 -> 544
10000 -> 3250 -> 3070 -> 2800
40000 -> 175, 551 -> 164, 005 -> 148,817
(So, for example, my 10,000 digit computations went from 17,678 seconds to 2,800 seconds.)
II predicted the 40000 digit run in method C would finish in (164005/1.09=150463) seconds, about 2 days. It took only 148817 seconds. Here is the work and results with a check against 100,000 computed digits from Method A:
In[32]:= N[(Timing[
M2 = Quiet[(NIntegrate[(Exp[Log[t]/t - Pi t/I]), {t, 1,
Infinity I}, WorkingPrecision -> 40000,
Method -> "Trapezoidal", MaxRecursion -> 15] + I/Pi)]]), 20]
Out[32]= {148817., 0.070776039311528803540 - 0.047380617070350786107 I}
In[33]:= M2100k - M2
Out[33]= 0.*10^-40001 + 0.*10^-40001 I
On 9.22.2021 I improved my timing for a 40,000 digit computation of M2 of the MKB constant using Method C by 1/2 an hour:
In[20]:= N[(Timing[(NIntegrate[(Exp[Log[t]/t - Pi t/I]), {t, 1,
Infinity I}, WorkingPrecision -> 40000,
Method -> "Trapezoidal", MaxRecursion -> 15] + I/Pi)])]
Out[20]= {147079., 0.070776 - 0.0473806 I}
200,000 digits of M2 are being computed using the MRB constant supercomputer in Methods A, B, and C.
Method A
g[x_] =
x^(1/x); t = (Timing[
test3 = -(I NIntegrate[(g[(1 + t I)])/(Exp[Pi t]), {t, 0,
Infinity}, WorkingPrecision -> 40000,
Method -> "Trapezoidal", MaxRecursion -> 15] + I/Pi)])[[1]];
Method B :
f[x_] = E^(I \[Pi] x) (1 - (1 + x)^(1/(1 + x))); Timing[NIntegrate[I (f[I t]), {t, 0, Infinity},
WorkingPrecision -> 40000, Method -> "Trapezoidal",
MaxRecursion -> 15]]
Method C :
N[(Timing[(NIntegrate[(Exp[Log[t]/t - Pi t/I]), {t, 1, Infinity I},
WorkingPrecision -> 40000, Method -> "Trapezoidal",
MaxRecursion -> 15] + I/Pi)])]
Where the MaxRecursion (M.R.) was no less than shown in this table.
Some of the AbelPlana records are shown here.
