[Oberon] Fun with π
Hans Klaver
hklaver at dds.nl
Tue Mar 15 00:44:50 CET 2022
Hi all,
Here are a few programs to have fun with π.
The first is just a joke, and shows how even Oberon code can be obfuscated ;-)
It calculates π by measuring the surface of a round disc. Feeding it larger discs
increases the number of decimals of π, but only very gradually: you need a ten times
larger disc to gain one decimal precision.
The second one is more interesting because the algorithm used produces quite a lot
of digits of pi using only ordinary integer operations; no need for floating point or
some arbitrary-precision datatype.
In Oberon System V5 it can calculate about 10,000 digits, and several hundreds of
thousands of digits can be produced by a slight variation of the algorithm run from
the command line (see the third program).
--
Hans Klaver
MODULE Westley; (* hk 3.14.2022 *)
(*
After an obfuscated C program by Brian Westley (IOCCC, 1988).
Select (the start of) a circle and middle click Westley.FOO @
Note that the underscore is not rendered well in Oberon10 font, so use Courier10 font to draw circles.
Separate two circles with "~~~".
*)
IMPORT Foo := In, FoO := Out;
PROCEDURE FOO*; (* just foo'ing around *)
CONST BAR = 9X; Bar = 0DX; baR = 20X; bar = 2DX; BaR = 5FX;
VAR foo: CHAR;
F,OO: REAL;
f,O0: INTEGER;
BEGIN f:=0;O0:=0;
Foo.Open; Foo.Char(foo);
WHILE Foo.Done & (foo # "~") DO
WHILE (foo = BAR) OR (foo = Bar) OR (foo = baR) DO
Foo.Char(foo)
ELSIF foo = BaR DO
DEC(f); DEC(O0);
Foo.Char(foo)
ELSIF foo = bar DO
DEC(f);
Foo.Char(foo); Foo.Char(foo)
END;
Foo.Char(foo); Foo.Char(foo)
END;
F := FLT(f);
OO := FLT(O0);
FoO.Real(4.0*(-F)/OO/OO, 12); FoO.Ln
END FOO;
END Westley.FOO @
Edit.ChangeFont Courier10.Fnt
_-_-_-_
_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_
_-_-_-_
~~~
_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_ _-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_
~~~
MODULE PiSpigot4; (* hk 3.14.2022 *) (* PiSpigot4.Compute *)
(*
Computes pi to (a few more than) N decimal digits.
After an obfuscated C program by Dik T. Winter at CWI Amsterdam,
(that computes 800 digits of pi) and its rewrite in C by Ben Lynn.
int a=10000,b,c=2800,d,e,f[2801],g;main(){for(;b-c;)f[b++]=a/5;
for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)for(b=c;d+=f[b]*a,
f[b]=d%--g,d/=g--,--b;d*=b);}
It uses a 'spigot' algorithm (dripping out digits as if from a tap,
in this special case four by four) and does not use the digits after
they are computed.
Based on evaluation of a mathematical series that is a mixed-radix
representation of pi with variable base = ( 1/3, 2/5, 3/7, 4/9, ...);
in this base pi is simply a list of 2's.
Lit.:
S. Rabinowitz & S. Wagon, A Spigot Algorithm for the Digits of π.
The American Mathematical Monthly, Vol. 102, No. 3 (Mar., 1995), pp. 195-203
Jörg Arndt and Christoph Haenel. 'π Unleashed', 2nd ed., 2001.
https://www.maa.org/sites/default/files/pdf/pubs/amm_supplements/Monthly_Reference_12.pdf
https://crypto.stanford.edu/pbc/notes/pi/code.html
*)
IMPORT Out;
CONST
N = 10000; (* number of digits of pi to compute; max. 10015 on Oberon V5 *)
len = (N DIV 4 + 1) * 14;
VAR f: ARRAY len OF INTEGER; (* array of 4-digit decimals *)
PROCEDURE Compute*;
CONST
a = 10000; (* new base, 4 decimal digits *)
aDIV5 = a DIV 5;
VAR
b, (* numerator previous base *)
c, (* counter *)
d, (* accumulator & carry *)
e, (* saves previous 4 digits *)
g: INTEGER; (* denominator previous base *)
PROCEDURE WritePaddedInt (x, n: INTEGER; ch: CHAR);
(* Write integer x in a field of at least n characters, padded
to the left with enough characters ch to fill the whole field *)
VAR i, x0: INTEGER;
BEGIN i := 0; x0 := x;
REPEAT x0 := x0 DIV 10; INC(i) UNTIL x0 = 0; (* x has i digits *)
WHILE n > i DO Out.Char(ch); DEC(n) END; Out.Int(x, 0)
END WritePaddedInt;
BEGIN ASSERT(N < 54933); (* at digit 54936 pi has 4 zero's at an index divisible by 4
which this version of the algorithm cannot handle properly *)
FOR b := 0 TO len - 1 DO f[b] := aDIV5 END; (* initialize the array with 2000's *)
e := 0;
FOR c := len - 1 TO 1 BY -14 DO (* outer loop: 4 digits / loop *)
d := 0; (* initialize accumulator & carry *)
b := c;
REPEAT (* inner loop: radix conversion *)
INC(d, f[b] * a);
g := b + b - 1; (* denominator previous base *)
f[b] := d MOD g;
d := d DIV g; (* save carry *)
DEC(b);
IF b # 0 THEN d := d * b END (* accum. := accum. * num. prev. base *)
UNTIL b = 0;
WritePaddedInt(e + d DIV a, 4, "0"); (* print previous 4 digits;
the leading zeros are significant! *)
IF c MOD 196 < 14 THEN Out.Ln END; (* to fit into the System.Log *)
e := d MOD a (* saves current 4 digits *)
END; Out.Ln
END Compute;
END PiSpigot4.Compute
MODULE PiSpigot;
(*
Oberon-07 program to compute N digits of π.
Compiled with OBNC for the POSIX console (Windows, macOS, Linux).
After a C program in the book 'π Unleashed' (p. 81), 2nd ed., 2001,
by Jörg Arndt and Christoph Haenel.
It uses a 'spigot' algorithm by Rabinowitz & Wagon, as modified by Haenel.
The algorithm drips out digits of π, one per loop, as if from a leaky tap,
and does not use the digits after they are computed. It only uses INTEGER,
no need for floating point or 'bignum'.
Based on evaluation of a mathematical series:
π = 2 + 1/3(2 + 2/5(2 + 3/7(2 + 4/9(2 + ...))))
in which π can be seen as a number in a number system with a variable base.
The 'mixed-radix' base b = (1/3, 2/5, 3/7, 4/9, ...) and now π itself
simply is written as an infinite list of 2's. The rest of the algorithm
is a radix conversion to the decimal number system.
The original algorithm (and a Pascal program) is described in the paper:
A Spigot Algorithm for the Digits of π, by Stanley Rabinowitz & Stan Wagon,
Amer. Math. Monthly, March 1995, 195-203.
https://www.maa.org/sites/default/files/pdf/pubs/amm_supplements/Monthly_Reference_12.pdf
Bug fixes for the algorithm and an earlier C program by Christoph Haenel:
https://www.jjj.de/hfloat/spigot.haenel.txt
For an explanation of how the algorithm works read the above R&W paper,
or see the book 'π-Unleashed' by Arndt & Haenel,
or see: http://www.pi314.net/eng/goutte.php
*)
IMPORT Out, Input;
CONST
N = 20000; (* max. N = 209602 before 'Segmentation fault' *)
len = 10*N DIV 3 + 1; (* ?? 10*N + 3 + 1 *)
VAR
a: ARRAY len OF INTEGER;
t: INTEGER; (* time *)
PROCEDURE GetClock (): INTEGER;
(* Time in seconds *)
BEGIN RETURN Input.Time() DIV Input.TimeUnit
END GetClock;
PROCEDURE Calculate (n: INTEGER);
(* Calculate the first n decimal digits of π *)
VAR
i, p, q,
nines: INTEGER; (* number of 9's next to each other *)
BEGIN (* ASSERT(n <= N); *)
i := 0; p := -1; q := -1; nines := 0;
FOR i := 0 TO len - 1 DO a[i] := 2 END; (* fill a[] with 2's *)
WHILE n >= 0 DO
(* Normalise and compensate for the very first digit *)
q := 0;
FOR i := len - 1 TO 1 BY -1 DO (* work from right to left *)
q := q + 10*a[i]; (* q := carry + 10*a[i] *)
a[i] := q MOD (i+i+1); (* a[i] := remainder *)
q := q DIV (i+i+1); (* quotient := FLOOR(q/(2*i+1))*i *)
q := q * i (* now q = carry *)
END;
(* The first digit is post-processed *)
q := q + 10*a[0]; (* q := carry + 10*a[0] *)
a[0] := q MOD 10;
q := q DIV 10; (* q: next provisional digit *)
IF q = 9 THEN INC(nines) (* q = 9: increment no. of 9's *)
ELSE
IF p >= 0 THEN Out.Int(p + q DIV 10, 0) END; (* p: previous provis. digit *)
IF n < nines THEN nines := n END; (* adjust digits to print *)
DEC(n, nines + 1);
DEC(nines);
WHILE nines >= 0 DO (* print 9's or 0's *)
IF q = 10 THEN Out.Int(0, 0) ELSE Out.Int(9, 0) END;
DEC(nines)
END;
nines := 0;
IF q = 10 THEN p := 0 ELSE p := q END (* set previous provisional digit *)
END
END; Out.Ln
END Calculate;
BEGIN
t := GetClock();
Calculate(N);
t := GetClock() - t;
Out.String("run-time = "); Out.Int(t, 0); Out.String(" s"); Out.Ln
END PiSpigot.
More information about the Oberon
mailing list