## Two more Pi programs from Valentin Albillo

bbc/electron apps, languages, utils, educational progs, demos + more
BigEd
Posts: 2995
Joined: Sun Jan 24, 2010 10:24 am
Location: West
Contact:

### Two more Pi programs from Valentin Albillo

Here's a simple, fast, and accurate one:

Code: Select all

``````   10 P=1
20 FOR I=1 TO 5
30 P=P+SIN(P)
40 PRINT P
50 NEXT

>RUN
1.84147099
2.80506171
3.13527633
3.14159261
3.14159265
``````
[Edit: Note, modified without permission: Valentin is clear that P should be initialised to 3]

Now, it does use SIN, but as noted in the thread, that's not inherently using pi to calculate pi, and it might even be computationally efficient. Even better, the accuracy [*] doubles [+] at each step, and at each step you only need a suitably (in)accurate calculation of SIN.

Here's a much much slower one, and also very much more obscure:

Code: Select all

``````   10 FOR K=1 TO 7
20 N%=10^K
30 S%=0
40 FOR I=1 TO N%
50 S%=S%+INT(RND(1)/RND(1)+.5)MOD2
60 NEXT
70 PRINT N%,1+4*S%/N%
80 NEXT
``````
(This being a tweaked version - I'll assume a typo in the original. (Edit: I am assured there is no typo there.))

10000 iterations is as much as I'd like to run on a 2MHz 6502, and that only gives us a couple of digits of accuracy. I can run deeper on a faster platform, but it doesn't help accuracy much! Note that the actual results will differ from run to run.

Code: Select all

``````>RUN
10       2.2
100       3.2
1000     3.216
10000     3.164
100000   3.13984
1000000  3.139792
10000000 3.1409552
``````
Previously:
Calculating digits of Pi in Basic
Computing pi in a very unusual way

[*] Edit - oops, as noted below, I should have said number of digits, or bits
[+] Edit - oops, as noted below, it triples rather than doubles
Last edited by BigEd on Wed Mar 18, 2020 9:42 am, edited 4 times in total.

Richard Russell
Posts: 1268
Joined: Sun Feb 27, 2011 10:35 am
Location: Downham Market, Norfolk
Contact:

### Re: Two more Pi programs from Valentin Albillo

BigEd wrote:
Sun Mar 15, 2020 10:32 am
Here's a simple, fast, and accurate one:
That's neat, and what is astonishing is that after only 5 iterations not only is it accurate to the 9 significant figures of 6502 BBC BASIC, but it's also accurate to the 19 significant figures of BBC BASIC for Windows and BBC BASIC for SDL (when running on an x86 platform)!

Code: Select all

``````     1 @%=&1313
10 P=1
20 FOR I=1 TO 5
30 P+=SIN(P)
40 PRINT P
50 NEXT
``````
gives:

Code: Select all

``````1.841470984807896507
2.805061709349729914
3.135276332899716
3.14159261159065315
3.141592653589793238
``````
and PRINT PI gives:

Code: Select all

``````3.141592653589793238
``````
My mind is boggled.

BigEd
Posts: 2995
Joined: Sun Jan 24, 2010 10:24 am
Location: West
Contact:

### Re: Two more Pi programs from Valentin Albillo

Spectacular! On re-rereading the thread, I see the accuracy [*] increases threefold on each iteration - I'd mistakenly said it doubles.

[*] see below - I should have said number of correct digits (or bits)
Last edited by BigEd on Sun Mar 15, 2020 2:05 pm, edited 1 time in total.

Richard Russell
Posts: 1268
Joined: Sun Feb 27, 2011 10:35 am
Location: Downham Market, Norfolk
Contact:

### Re: Two more Pi programs from Valentin Albillo

BigEd wrote:
Sun Mar 15, 2020 12:35 pm
Spectacular! On re-rereading the thread, I see the accuracy increases threefold on each iteration - I'd mistakenly said it doubles.
In fact, AIUI, it's the number of significant figures (or bits) which triples at each iteration, which is a vastly faster convergence than what I originally understood you to be using the word "accuracy" to mean.

BigEd
Posts: 2995
Joined: Sun Jan 24, 2010 10:24 am
Location: West
Contact:

### Re: Two more Pi programs from Valentin Albillo

Oh, good point, I used entirely the wrong word there!

RobC
Posts: 2813
Joined: Sat Sep 01, 2007 10:41 pm
Contact:

### Re: Two more Pi programs from Valentin Albillo

You can visualise the convergence of these dynamical systems, where you iterate x=f(x), by plotting the function (y=x+sin(x) in this case) and the line y=x. You then choose your start point on the y=x line, draw a line vertically to meet the curve f(x) and then draw a horizontal line back to the y=x line. You then repeat this process.

The process will converge at the points where the two lines meet i.e. f(x)=x which happens at every multiple of Pi in this case. Due to the shape of the curve, convergence is stable at odd multiples of Pi (so values around odd multiples get closer to them) and unstable at even multiples (so any slight displacement will cause the process to move away).

Here's a plot illustrating the process with x=0.5 as the starting point:

BigEd
Posts: 2995
Joined: Sun Jan 24, 2010 10:24 am
Location: West
Contact:

### Re: Two more Pi programs from Valentin Albillo

Nice plot! Valentin wasn't best pleased that I started the convergence at 1: for sufficiently small starting values, the convergence will be very slow at first. But I like 1 rather than 3 because it appears to encode less information about the final result.

scruss
Posts: 216
Joined: Sun Jul 01, 2018 4:12 pm
Location: Toronto
Contact:

### Re: Two more Pi programs from Valentin Albillo

It's neat, but it's not doing anything magic. All it's doing is exploiting the fact that sin π = 0. If you start the program with P > π, you'll notice it converges on the nearest multiple of π:

Code: Select all

``````   10 P=RND(100)
15 PRINT "P= ",P
20 FOR I=1 TO 5
30   P=P+SIN(P)
40   PRINT I," ",P," ",P/PI
50 NEXT
``````
with three sample runs as follows:

Code: Select all

``````>RUN
P=                32
1          32.5514267          10.3614409
2          33.4581719          10.6500669
3           34.349083          10.9336527
4          34.5560132          10.9995206
5          34.5575192                  11
>RUN
P=               100
1          99.4936344          31.6698074
2          98.6325845          31.3957268
3          97.6857621          31.0943438
4          97.3936927          31.0013752
5          97.3893723                  31
>RUN
P=                90
1          90.8939967          28.9324577
2          91.1045982          28.9994943
3           91.106187                  29
4           91.106187                  29
5           91.106187                  29
``````
You're effectively dropping marbles into a cosine-shaped trough, since sine is the slope of the cosine graph. I think differential calculus is already known …
Last edited by scruss on Sat Mar 21, 2020 7:02 pm, edited 1 time in total.

BigEd
Posts: 2995
Joined: Sun Jan 24, 2010 10:24 am
Location: West
Contact:

### Re: Two more Pi programs from Valentin Albillo

Yes, indeed, the iteration will also converge to other odd multiples of pi.

I'd say it's pretty good - you don't have to start especially near pi, you don't need pi in your calculation of sin(), and you get three times as many digits each time you iterate. And as Valentin points out, the calculation of sin() doesn't need to have full precision until the final iteration. Oh, and also, sin() isn't necessarily that expensive to compute, if you really roll up your algorithmic sleeves.

RobC
Posts: 2813
Joined: Sat Sep 01, 2007 10:41 pm
Contact:

### Re: Two more Pi programs from Valentin Albillo

scruss wrote:
Sat Mar 21, 2020 6:42 pm
It's neat, but it's not doing anything magic. All it's doing is exploiting the fact that sin π = 0. If you start the program with P > π, you'll notice it converges on the nearest multiple of π.
Not quite - it'll generally converge to the nearest odd multiple of π unless you land bang on an even multiple. This is because x+sin(x) > x in the interval (2nπ, (2n+1)π) and x+sin(x) < x in the interval ((2n+1)π, 2nπ). So, the sequence will move away from even multiples of π and towards odd multiples of π.

scruss
Posts: 216
Joined: Sun Jul 01, 2018 4:12 pm
Location: Toronto
Contact:

### Re: Two more Pi programs from Valentin Albillo

okay, so missing the odd multiple detail blunted my message a bit, but it's really only doing a simple root finding. It's going to appear to converge at times roughly cubically because the sine curve is roughly cubic. SIN is slightly more expensive than ATN (on a BBC Micro, at least) so we don't gain anything by using it.

Even using the demonstrably worst approximation of the sine function y = 2 - 2x/π (for 0 < x ≤ 3π/2) is going to converge onto π. You might say you don't need π in your calculation of sin(), but π is intrinsic to trig functions. You're not using it explicitly in 4*ATN(1) either.

BigEd
Posts: 2995
Joined: Sun Jan 24, 2010 10:24 am
Location: West
Contact:

### Re: Two more Pi programs from Valentin Albillo

If you're not impressed by this sin() approach, that's fine. I am!

If you're interested in a deeper view of pi, I very much recommend this article:
What is 𝜋? (and while we're at it, what's 𝑒?)