For some reason I started thinking about summing up some reciprocals, specifically computing pi using the famously slow (but historically important) series by Leibniz
So I did some programming, and here's what happened next. I don't include any code, so you can enjoy writing your own - nothing here takes more than a dozen lines.
So, if it weren't so terribly slow (wikipedia says "Calculating π to 10 correct decimal places using direct summation of the series requires precisely five billion terms"), we'd just be alternately adding and subtracting the reciprocals of the odd numbers, and multiplying by four. If we do just that, we take ages to get only very modest results. But it is very simple to program...
It turns out that we can add up ten thousand terms in a minute, in Basic, on a Beeb, we get just a few digits of accuracy, and we don't get much more if we double the work:
I reckoned that the alternate partial sums were always overshooting or undershooting but by just a little less each time... should be possible to get a better result by adding the final term halved. After some experimentation and a little thought, I realised that was just the same as averaging the last two partial sums.
Using that idea just 100 terms, and just 1 second, gets us a bit of accuracy:
But these results too seem to undershoot then overshoot, so it looks like we could average these partial sums too... in fact I think you can keep doing that, and it turns out to be, I think, something a little like Euler's acceleration idea. Except that was proved whereas this is just playing.
So here we see we can get plenty of accuracy in just a screen's worth of results, in a couple of seconds:
It turns out that if you do enough maths, you can get all the averaging done in the construction of the series, and compute a different series instead, which is equivalent. But that's then a program of a different form - simpler and faster, but different. For one thing, it computes half pi rather than a quarter pi.
I did wonder a bit about perhaps trying to use Kahan's method to improve the accuracy of summation, but once I was adding rather few terms that went out the window. Likewise thoughts about using a double precision idea.
(For more on the history of the series, previous work, and similar things, perhaps see here.)
Computing pi with Leibniz
Re: Computing pi with Leibniz
Some time has passed, so I thought I'd share some of my basic Basic efforts.
Before I do, I'd like to mention this amazing pi approximation, found by Gerson and posted over on hpmuseum. In Basic: x=16*LN878:y=LN(x/LNx):PRINT y;'ABS(PI-y).
So, onto my meanderings. The trajectory, broadly, was either to increase accuracy or decrease the runtime or increase the usefulness of the output.
Super simple starting code (owlet):
Next for 1000 iterations in 6.6s showing how averaging two results is a big improvement in accuracy:
Another level of averaging is even better - see here.
If we run the computation the other way, adding up the smallest terms first, we avoid a lot of rounding errors and it looks like we need rather fewer terms - 1400 terms, 3 seconds, here):
Then I did this which seemed even better but possibly it's done an extra level of averaging.
I've a note-to-self saying
My next claim is
Here's the code for the "if you do enough maths" reformulation - runs in a second:
And back to the Leibniz series, here's code which runs in 1.6 seconds and tabulates the partial sums and their averages and the difference to pi:
Hope this is of interest to someone! Please post any observations, improvements, or related experiments.
Before I do, I'd like to mention this amazing pi approximation, found by Gerson and posted over on hpmuseum. In Basic: x=16*LN878:y=LN(x/LNx):PRINT y;'ABS(PI-y).
So, onto my meanderings. The trajectory, broadly, was either to increase accuracy or decrease the runtime or increase the usefulness of the output.
Super simple starting code (owlet):
Code: Select all
REM MODE 4
@%=&1A0A
D=1
S=1
M=1
REPEAT
D=D+2
M=-M
S1=S
S=S+M/D
PRINT 4*S,2*S+2*S1
UNTILD>20
Code: Select all
REM MODE 4
T=TIME
@%=&1A0A
D%=1
S=1
M%=1
FORI%=1TO1000
D%=D%+2
M%=-M%
S1=S
S=S+M%/D%
IF I%<990 NEXT ELSE PRINT 4*S,2*S+2*S1
NEXT
PRINT (TIME-T)/100;" seconds"
If we run the computation the other way, adding up the smallest terms first, we avoid a lot of rounding errors and it looks like we need rather fewer terms - 1400 terms, 3 seconds, here):
Code: Select all
REM MODE 4
T=TIME
@%=&2090C
S=0:P=0:Q=0
M%=-1
FORD%=1401TO1 STEP -2
M%=-M%
S=S+M%/D%
NEXT
P=4*S
Q=4*(S-1/1403)
R=4*(S-1/1403+1/1405)
PRINTP,(P+Q)/2,((P+Q)/2+(Q+R)/2)/2
@%=&202
PRINT (TIME-T)/100;" seconds"
I've a note-to-self saying
but I can't vouch for that.turns outis a tad more accurate, summing the head in rolling up orderCode: Select all
P1=P:P=4*(1+(-1/3+(1/5+(-1/7+(1/9+S+M%/(D%+2)/2)))))
My next claim is
Code here:the triple averaging gets an exact pi at D=367 which is 184 terms or so
and just under 1.5 seconds
Code: Select all
T=TIME
@%=&2090C
S=0:P=0:R=0
M%=-1
FORD%=11TO399 STEP 2
S=S+M%/D%
M%=-M%
IF D%<346 NEXT
P1=P:P=4*(1+(-1/3+(1/5+(-1/7+(1/9+S+M%/(D%+2)/2)))))
R1=R:R=(P+P1)/2
IFR1>3 PRINT (R+R1)/2,(R+R1)/2-PI,P-PI
IF (R+R1)/2<>PI NEXT
@%=&A0A
PRINT D%'(TIME-T)/100;" seconds"
Code: Select all
MODE 4:T%=TIME
@%=&2090D
N=1:D=3:S=1:T=1
REPEAT
PRINT 2*S,2*S-PI
T=T*N/D
S=S+T
N=N+1:D=D+2
UNTIL N>29
@%=&A0A
PRINT (TIME-T%)/100;" seconds"
Code: Select all
MODE 4:T=TIME
@%=&20609
S=0:P=0:R=0:Q=0
M%=1
FORD%=1TO55 STEP 2
S=S+M%/D%
M%=-M%
P1=P:P=4*S
Q1=Q:Q=(P+P1)/2
R1=R:R=(Q+Q1)/2
PRINT (P+P1)/2,(Q+Q1)/2,(R+R1)/2,ABS(PI-(R+R1)/2)
NEXT
@%=&A0A
PRINT (TIME-T)/100;" seconds"
- TobyLobster
- Posts: 618
- Joined: Sat Aug 31, 2019 7:58 am
- Contact:
Re: Computing pi with Leibniz
For the 'runs in a second' code, most of the time appears to be printing the partial results. If you move line 50 after the REPEAT-UNTIL loop, the runtime decreases to 0.3 seconds.
Re: Computing pi with Leibniz
Nice experiments, it's good to see the various evolutions. In the "enough maths" version, I think you could initialise S and T to 2 and then you wouldn't need to double S before printing it.
Re: Computing pi with Leibniz
Thanks both - there's another little thing which is to make N and D integer variables. And a FOR..NEXT is usually quicker than REPEAT..UNTIL but in this case we don't loop many times. And then squashing onto fewer lines can be a small win...
(On the sidetrack of that spectacular pi approximation, good to 10 digits, Gerson has elaborated in the linked thread as to how he went about finding it. He also linked to an interesting paper about how best to measure the efficiency of an approximation formula in terms of size and accuracy - not just counting the digits of the constants, which is what I was inclined to do.)
(On the sidetrack of that spectacular pi approximation, good to 10 digits, Gerson has elaborated in the linked thread as to how he went about finding it. He also linked to an interesting paper about how best to measure the efficiency of an approximation formula in terms of size and accuracy - not just counting the digits of the constants, which is what I was inclined to do.)