# [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 *)
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)

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.

```