Dear Marvin,

thank you for your nice answer! I need to add the reference, where I found the Shanks transformation:

Carl M. Bender, Steven A. Orszag; Advanced Mathematical Methods for Scientists and Engineers; Springer, 1999

An outline of the idea goes like this: Suppose the sequence of the partial sums $A_n$ behaves like: $$A_n = A + \alpha q^n$$ This actually describes roughly the situation in case of the MRB. (If it were an exact description, we had the "exact" solution after just one single Shanks.) The wanted value of course is $A$, and because we have a three parameter model, we need three terms of the sequence to calculate $A$, e.g. $A_{n-1},A_{n},A_{n+1}$. Then one gets: $$A = \frac{A_{n+1}A_{n-1}-A_{n}^2}{A_{n+1}+A_{n-1}-2A_{n}}$$ I find the possible improvement by this simple transformation most impressive. It can be illustrated like so (using the definition of sac from my first post):

And - at last - to answer your question: I do not think that Mathematica has an issue with more Shanks iterations. But then of course one needs to use more terms and more digits:

ClearAll["Global`*"]
m = Table[N[(-1)^k (k^(1/k) - 1), 2000], {k, 1, 2000}];
(* partial sums of the series *)
am = Accumulate@m;
shanks[ac_List] := Table[(ac[[n + 1]] ac[[n - 1]] - ac[[n]]^2)/(ac[[n + 1]] + ac[[n - 1]] - 2 ac[[n]]), {n, 2, Length[ac] - 1}]
sac = NestList[shanks, am, 24];
ListLogPlot[Abs[Differences /@ sac], Joined -> True, GridLines -> Automatic, ImageSize -> Large]

giving (after two or three seconds on my old PC):

Then for an approximation of your MRB constant we get:

err = N[Abs[sac[[-1, -1]] - sac[[-1, -2]]], 5]
(*  Out:  6.2079*^-146  *)
numberOfDigits = Floor[-Log[10, err]]
(*  Out:  145  *)
mrb = N[sac[[-1, -1]], numberOfDigits]
(*  Out:  \
0.18785964246206712024851793405427323005590309490013878617200468408947\
7231564660213703296654433107496903842345856258019061231370094759226630\
43892937955325691395341832`145.   *)

Regards -- Henrik

Dear Marco,

thank you for your nice and encouraging comments! Yes, it is almost magical - and it shows a sort of "beauty of simplicity"! I learned about Shanks transformation when reading - or trying to read - Benders book (which for me is somewhat challenging) and watching his lectures (which I find excellent!). Bender shows lots of most interesting methods I had never heard of before, particularly not in my math courses (about 30 years ago ...).

Best wishes -- Henrik

Dear Marvin,

as a preliminary remark I need to say that I have not followed this obviously long ongoing discussion closely, because I am not really interested in this kind of calculations. So forgive me please if my idea was posted long time ago or if my post is completely off topic.

It is obvious from the definition that your MRB constant emerges from an alternating series. Any partial sum has an equivalent behavior (in the following I am using not tens of thousands of terms but just the first 100 terms - as a proof of concept):

(* first 100 terms of the series: *)    
m100 = Table[N[(-1)^k (k^(1/k) - 1), 200], {k, 1, 100}];
(* partial sums of the series *)
am100 = Accumulate@m100;
ListPlot[am100, GridLines -> Automatic]

This is a typical situation where a Shanks transformation seems promising:

shanks[ac_List] := 
 Table[(ac[[n + 1]] ac[[n - 1]] - ac[[n]]^2)/(ac[[n + 1]] + ac[[n - 1]] - 2 ac[[n]]), {n, 2, Length[ac] - 1}]

which already leads to a huge improvement:

The nice thing about Shanks transformation is that it can be applied iteratively, leading to further improvements, e.g.:

sac = NestList[shanks, am100, 7];

The difference between the last and the second last term can be regarded as a measure of the error. What we finally get is:

Doing these calculations is just a matter of milliseconds.

Finally I can not resist in citing the great Carl Bender in his very first lecture on mathematical physics (at about min 43):

If you need to find the sum of a series, the worst possible technique is to add the numbers of the series together one at a time ...

Regards -- Henrik

Let g[n] be the error in using f, as you use 10^n terms, as shown in the original post; using the following piece of code:

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

N::meprec: Internal precision limit $MaxExtraPrecision = 50.` reached while evaluating -E^(2 I Interval[{0,\[Pi]}])+E^(2 I Interval[{0,\[Pi]}]). >>

f[k_] := Log[(k - Log[k])/k];

ClearAll[g]; 
g[n_] := Module[{}, 
  m - (NSum[(-1)^k (k^(1/k) - 1 + f[k]), {k, 1, 10^n}, 
      WorkingPrecision -> Floor[3 n + 100], 
      Method -> "AlternatingSigns"] - 
     NSum[(-1)^k*f[k], {k, 1, Infinity}, 
      WorkingPrecision -> Floor[3 n + 100], 
      Method -> "AlternatingSigns"])]

. You find there is a pattern to g[n]

For[n = 2, n < 20, 
 Print[RootApproximant[
   MantissaExponent[g[10 n]/g[10 (n + 1)]][[1]]]]; n++]

8/27

27/64

64/125

125/216

216/343

343/512

512/729

729/1000

1000/1331

1331/1728

1728/2197

2197/2744

2744/3375

3375/4096

4096/4913

4913/5832

5832/6859

6859/8000

, and find g[10 n]/g[10 (n + 1)] ->n^3/(n+1)^3 *10 to a factor:

For[n = 2, n < 20, Print[
  RootApproximant[
    MantissaExponent[
      g[10 n]/g[10 (n + 1)]
      ][[1]]
    ] - n^3/(n + 1)^3
  ]; n++]

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0