Another short ALGOL 68 program that I wrote last year while trying to get to grips with the language again, and as before this was based on a existing PASCAL program that I’d written previously.

The original program used Gregory’s Series to find Pi but this converges very slowly and while it is useful for using up spare clock cycles it isn’t much good for finding Pi.

After doing a bit of research I found an alternative related series that converges much faster and modified the original program to use that instead.

```
PROGRAM series
{ Program to calculate Pi - 05 Aug 12
Copyright (C) 2012 - MEJT
This program is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the Free
Software Foundation, either version 3 of the License, or (at your option)
any later version.
This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of ERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License along
with this program. If not, see http://www.gnu.org/licenses/
Uses a Gregory's series expansion to find Pi, an example of which is the
Leibnitz Formula.
Pi 1 1 1 1
---- = 1 - - + - - - + - - ...
4 3 5 7 9
This series is very very slow to converge so the following series
expansion has been used.
1 Pi 1 1 1
atan(-------) = -- = ------- - --------------- + ----------------- - ...
sqrt(3) 6 sqrt(3) 3 * 3 * sqrt(3) 5 * 3^2 * sqrt(3)
By expanding this and simplifying we get a series that converges much
faster and appears to yeld aproximatly 16 decimal places in just 31
iterations.
Pi 1 1 1
------------- = 1 - ----- + --------- - ------------- + ...
2 * sqrt(3) 3 * 3 5 * 3 * 3 7 * 3 * 3 * 3
When executed this program prints successive approximations for Pi.
1 1 3.4641016151377546
2 9 3.0792014356780041
3 45 3.1561814715699542
4 189 3.1378528915956804
5 729 3.1426047456630847
6 2673 3.1413087854628836
7 9477 3.1416743126988377
8 32805 3.1415687159417843
9 111537 3.1415997738115058
10 373977 3.1415905109380801
11 1240029 3.1415933045030816
12 4074381 3.1415924542876464
13 13286025 3.1415927150203798
14 43046721 3.1415926345473139
15 138706101 3.1415926595217137
16 444816117 3.1415926517339976
17 1420541793 3.1415926541725754
18 4519905705 3.1415926534061652
19 14334558093 3.1415926536478261
20 45328197213 3.1415926535714034
21 142958160441 3.1415926535956350
22 449795187729 3.1415926535879335
23 1412147682405 3.1415926535903866
24 4424729404869 3.1415926535896037
25 13839047287569 3.1415926535898540
26 43211719081593 3.1415926535897738
27 134718888901437 3.1415926535897996
28 419407861674285 3.1415926535897913
29 1303977169932777 3.1415926535897939
30 4049192264528097 3.1415926535897931
31 12559359057773589 3.1415926535897933
32 38913423965888661 3.1415926535897932
}
BEGIN
REAL term := 1.0, factor;
REAL number := 1.0, sum := 0.0;
INT count := 0, limit := 31;
factor := 2 * sqrt (3.0);
WHILE (count <= limit) DO
term := (2 * count + 1) * 3.0 ^ (count);
IF ODD count THEN
sum := sum - (1.0 / term)
ELSE
sum := sum + (1.0 / term)
FI;
count +:= 1;
write ((fixed (count, -3, 0), " ", fixed (term, -18, 0), " "));
write ((fixed (sum * factor, -18, 16), " ", newline))
OD
END
FINISH
```

Advertisements