Two more Pi programs from Valentin Albillo

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

Two more Pi programs from Valentin Albillo

Post by BigEd » Sun Mar 15, 2020 10:32 am

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.

User avatar
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

Post by Richard Russell » Sun Mar 15, 2020 12:13 pm

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.

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

Re: Two more Pi programs from Valentin Albillo

Post by BigEd » 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.

[*] 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.

User avatar
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

Post by Richard Russell » Sun Mar 15, 2020 2:01 pm

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.

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

Re: Two more Pi programs from Valentin Albillo

Post by BigEd » Sun Mar 15, 2020 2:04 pm

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

Post by RobC » Sat Mar 21, 2020 12:19 pm

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:
convergence.png

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

Re: Two more Pi programs from Valentin Albillo

Post by BigEd » Sat Mar 21, 2020 12:59 pm

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

Post by scruss » 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 π:

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.

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

Re: Two more Pi programs from Valentin Albillo

Post by BigEd » Sat Mar 21, 2020 7:01 pm

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

Post by RobC » Sat Mar 21, 2020 7:10 pm

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

Post by scruss » Sat Mar 21, 2020 10:08 pm

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.

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

Re: Two more Pi programs from Valentin Albillo

Post by BigEd » Sat Mar 21, 2020 10:14 pm

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 𝑒?)

Post Reply

Return to “8-bit acorn software: other”