Π day

[robve]

(03-14-2022 11:06 PM)ttw Wrote:  One can often apply various sequence transformations to speed up these series. A couple of new methods are given here.

https://arxiv.org/pdf/1702.07199.pdf

https://people.mpim-bonn.mpg.de/zagier/f...lltext.pdf


In most record-setting attempts, the desired result is a fraction with really big integers, even bigger than 355/113. An (to me, anyway) amusing feature is that most of the time is taken by the final long division. Third order, first order, fourth order, etc., methods all take about the same amount of time as the time for arithmetic grows pretty fast.

I recall there was also a HP Forum post on accelerating alternating sums with Albert Chan suggesting to use Aitken extrapolation correction. In addition, I observed that summing the terms in reverse order may sometimes improve accuracy.

Here is the code I wrote a year ago (for my collection):

' SUMALT - summation of infinite alternating series
' Euler transformation:
' https://en.wikipedia.org/wiki/Series_acceleration
' https://mathworld.wolfram.com/EulersSeri...ation.html
' Adapted from:
' https://albillo.hpcalc.org/programs/HP%2...Series.pdf
' Combined with Aitken extrapolation correction applied to series S2:
' https://en.wikipedia.org/wiki/Aitken%27s...ed_process
' Combined with reverse term summation
' Performs N+D+2 function evaluations and D*(D+1)/2+2*(D+1) table operations
' The term y(i) to sum is defined in line 200 using I as index and returning value Y

' VARIABLES
' N number of terms
' D order of differences, e.g. 7 should suffice
' S S1 and final sum_i=0^inf (-1)^i*y(i)
' A() auto array with difference terms t[0] to t[d] stored in A(27) to A(27+D)
' I sum iterator
' J iterator
' X S2
' Y value of y(i)
' Z scratch

' driver program
10 N=10,D=7
20 INPUT "N=";N,"D=";D
30 GOSUB 100: PRINT S: END

' init differences t[0] to t[d]
100 FOR I=N+1 TO N+D+1: GOSUB 200: A(26+I-N)=Y: NEXT I
' apply Euler transformation to compute differences t[0] to t[d] and store in A(27) to A(27+D)
110 FOR I=1 TO D: FOR J=D TO I STEP -1: A(27+J)=A(27+J)-A(26+J): NEXT J: NEXT I
120 Z=2
130 FOR I=0 TO D: A(27+I)=A(27+I)/Z,Z=-Z-Z: NEXT I
' apply Aitken extrapolation as a correction to initialize S2
140 S=0,X=-A(27+D)*A(27+D)/(A(27+D)-A(26+D)),Z=1
' compute S = S1 = sum_0^n (-1)^i*y(i)
150 FOR I=0 TO N: GOSUB 200: S=S+Z*Y,Z=-Z: NEXT I
' compute X = S2 = sum_0^d t[i] in reverse order to retain accuracy
160 FOR I=D TO 0 STEP -1: X=X+A(27+I): NEXT I
' S = S1+S2 if N is odd, S = S1-S2 otherwise, where Z=1 if N is odd or -1
170 S=S+Z*X
180 RETURN
' y(i)
200 Y=1/(3^I*(2*I+1)): RETURN

Line 200 is the alternating \( \displaystyle \frac{1}{3^i\,(2i+1)} \) term.

This method only takes N=5 to produce 8 decimal places of \( \pi \) with order of differences D=3:

RUN
N=5
D=3
         9.068996707E-01
9.068996707E-01*6/√3
3.141592614

Not bad!

- Rob

Edit note/warning: there is no reason to push this method hard with N=10 and D=7, which will require N+D+2=18 function evaluations and D*(D+1)/2+2(D+1)=44 tabulation and summation steps! That's far more CPU power than the N=17 steps for the Sharp series to converge to 10 decimal places.

Edit 2 The double precision version for the Sharp PC-1475 to compute up to 20 decimal places:

' driver program
10 DEFDBL: DEFDBL A,S,X-Z: N=10,D=7
20 INPUT "N=";N,"D=";D
30 ERASE A: DIM A(D)
40 GOSUB 100: PRINT S: END

' init differences t[0] to t[d]
100 FOR I=N+1 TO N+D+1: GOSUB 200: A(I-N-1)=Y: NEXT I
' apply Euler transformation to compute differences t[0] to t[d] and store in A(0) to A(D)
110 FOR I=1 TO D: FOR J=D TO I STEP -1: A(J)=A(J)-A(J-1): NEXT J: NEXT I
120 Z=2
130 FOR I=0 TO D: A(I)=A(I)/Z,Z=-Z-Z: NEXT I
' apply Aitken extrapolation as a correction to initialize S2
140 S=0,X=-A(D)*A(D)/(A(D)-A(D-1)),Z=1
' compute S = S1 = sum_0^n (-1)^i*y(i)
150 FOR I=0 TO N: GOSUB 200: S=S+Z*Y,Z=-Z: NEXT I
' compute X = S2 = sum_0^d t[i] in reverse order to retain accuracy
160 FOR I=D TO 0 STEP -1: X=X+A(I): NEXT I
' S = S1+S2 if N is odd, S = S1-S2 otherwise, where Z=1 if N is odd or -1
170 S=S+Z*X
180 RETURN

RUN
N=28
D=1
0.9068996821171089253
0.9068996821171089253#*6/√3
3.14159265358979793285