Message Boards Message Boards

1
|
9275 Views
|
8 Replies
|
9 Total Likes
View groups...
Share
Share this post:

Record breaking direct summation of MRB constant

The MRB constant is approximated by

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

.

This method is pitifully inefficient in its direct summation of MRB constant ie. -1^(1/1)+2^(1/2)-3^(1/3)+... . Here is how much error there is after 10,000 terms:

 m - (NSum[(-1)^k (k^(1/k) - 1), {k, 1, 10^4}, 
   WorkingPrecision -> 200, Method -> "AlternatingSigns"])

(* -0.
0004607086148833843903798558880939839952579916174598636084830590133354\
3604452568023408196838217734393902621089862492299952504355243485749500\
9513789757123430499002402147185148514618052628511533818725*)

.

Due to the work of Daniel Lichtblau, me, and the late Richard Crandall, found in the middle of my most successful post,

http://community.wolfram.com/groups/-/m/t/366628?ppauth=zGnyw9gs

, in the second reply that starts with the phrase "Daniel Lichtblau and others," without the quotes; I think I found the optimum direct summation method of the m= -1^(1/1)+f1 +2^(1/2)+f2 -3^(1/3)+f3+ f(k)(x)+ ...- sum[f[k]]. That is 

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

. Then subtract Sum[(-1)^k Log[(k - Log[k])/k] . Here is an example showing the much smaller error in adding the first 10^4 terms:

     f[k_] := Log[(k - Log[k])/
        k]; m - (NSum[(-1)^k (k^(1/k) - 1 + f[k]), {k, 1, 10^4}, 
         WorkingPrecision -> 200, Method -> "AlternatingSigns"] - 
        NSum[(-1)^k *f[k], {k, 1, Infinity}, WorkingPrecision -> 200, 
         Method -> "AlternatingSigns"])

(* 6.
5176019395153886286218356168146776861946664225985899368189404704176944\
5768612992269488431506197821023123569647840336109019435539071457039552\
26445450414770330729207675435060098684928265988*10^-11*)

Why don't see if you can find a better algorithm!
POSTED BY: Marvin Ray Burns
8 Replies
Sort By:

Henrik et al,

I have a big job ahead of me if I try to use Shanks to break my record of calculated digits of the MRB constant!! The terms needed to break my record of 3,014,991 digits seems impractical to list! Maybe I just need a new way of computing the Shanks' terms (someway of using known acceleration methods to accumulate them).

POSTED BY: Marvin Ray Burns

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
POSTED BY: Henrik Schachner

Dear Henrik,

that's a beautiful piece of mathemagic.

Thanks for sharing,

Marco
POSTED BY: Marco Thiel

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):

enter image description here

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):

enter image description here

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
POSTED BY: Henrik Schachner

Thank you Henrik! I think the Shanks transformation will help me. It seems, however, Mathematica has an issue with 10's and 100's of iterations of Shanks.
POSTED BY: Marvin Ray Burns

Again. 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, and then it appears to be the case that the g[10 n]'s are rational factors of Log[10]^3, to better and better precision as n gets large:

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



 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"])]



 Table[
 RootApproximant[MantissaExponent[Log[10]^3/g[10 n]][[1]]], {n, 2, 20}]

(* {3/20, 4/9, 3/16, 24/25, 5/9, 120/343, 15/64, 40/243, 3/25, \ 1200/1331, 25/36, 1200/2197, 150/343, 16/45, 75/256, 1200/4913, \ 50/243, 1200/6859, 3/20} *)
POSTED BY: Marvin Ray Burns

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]

enter image description here

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:

enter image description here

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:

enter image description here

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
POSTED BY: Henrik Schachner

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
POSTED BY: Marvin Ray Burns
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract