Computing Pi and e: Fourteen Classical Algorithms in Exact Arithmetic
The Oldest Calculation in Mathematics
Few numbers have occupied mathematicians as persistently as 
, the ratio of a
circle’s circumference to its diameter, and 
, the base of the natural
logarithm. Archimedes trapped between two polygons around 250 BCE,
proving 
with nothing but geometry and
patience. The Madhava school in Kerala and later Gregory and Leibniz in Europe
discovered that falls out of infinitely long sums of simple fractions.
Viete wrote the first infinite product in 1593, Wallis another in 1656, Machin
found a formula in 1706 that let him compute 100 digits by hand, and Euler gave
us the continued fraction for in 1737. The story did not end with the
computer age: the Chudnovsky brothers published the series behind every modern
record computation of in 1988, and Bailey, Borwein, and Plouffe stunned
everyone in 1995 with a formula that produces a single hexadecimal digit of
without computing any of the digits before it.
Why should a programmer care about this history? Because these algorithms are a perfect laboratory for three ideas that matter far beyond recreational mathematics.
First, convergence rates. Each algorithm is a process that produces better and better approximations, and the processes differ by enormous factors. The Leibniz series needs roughly ten times more terms for each additional correct digit, while the Chudnovsky series buys about 14 digits per term. Running both side by side makes the abstract idea of a convergence rate visceral.
Second, exact arithmetic. A 64 bit floating point number carries about 16
decimal digits of precision, so it is useless for computing the 1000th digit of
. The classical workaround, and the one used throughout this chapter, is
to represent a real number as an exact integer 
paired with a scale 
,
interpreted as the value 
. All arithmetic is done on integers, so
there is no rounding error to accumulate; the only error is the mathematical
error of truncating an infinite process, which we can estimate on paper.
Third, verification. If you compute 500 digits of , how do you know
they are right? The program in this chapter answers that question the way
numerical analysts do: compute the same constant twice with two mathematically
independent algorithms and count how many leading digits agree.
The example code for this chapter implements nine algorithms for and
five for , wraps them in a small registry with metadata about each method,
verifies every run against a trusted reference value, and presents the whole
thing through both an interactive console menu and a command line interface. A
RackUnit test suite pins the fast algorithms against published digit tables and
checks the structural guarantees of the slow ones. The code is in the files
math.rkt and math-tests.rkt.
A Tour of the Algorithms
Before looking at the code, it helps to have the formulas in front of you and to know what kind of behavior to expect from each one.
Infinite series for pi
The Gregory-Leibniz series is the simplest of all:
It is also nearly useless in practice. The error after 
terms is about
, so a million terms buys only about six digits. Nilakantha’s 14th
century improvement,
converges much faster because the denominators grow like 
; each tenfold
increase in terms yields about three extra digits.
Machin’s formula from 1706,
converts the problem into evaluating the arctangent series
at two small arguments,
where it converges rapidly. This was the standard hand computation method for
two centuries, and it is still good for about 1.4 digits per term.
Products and iterations for pi
The Wallis product and Viete’s nested radicals both express as an
infinite product:
Archimedes’ method is the most geometric. Inscribe and circumscribe hexagons in
a circle, then repeatedly double the number of sides. The perimeters of the
inscribed polygons climb toward from below while the circumscribed
perimeters descend toward it from above, so at every step is trapped
inside a known bracket. In the code this appears in its modern algebraic form as
a pair of interleaved harmonic and geometric means, 
and 
, starting from 
and 
.
The bracket guarantee is a property we can test directly.
The Gauss-Legendre AGM iteration replaces the geometric picture with pure
arithmetic: repeatedly replace 
and 
by their arithmetic and geometric
means, accumulate a correction term 
, and compute
. Its convergence is quadratic, meaning the number
of correct digits roughly doubles at every step. Four iterations already give
more than 40 digits.
The Chudnovsky series is the industrial strength method:
Each term contributes about 14.18 digits. The code evaluates it with binary splitting, a divide and conquer technique that computes the product tree of the factorials recursively instead of term by term, which is what makes the method fast enough to hold every modern record.
The BBP formula is the odd one out:
Because the weights are powers of 
, the fractional part of the partial
sum up to position 
determines hexadecimal digit of directly.
You can ask for the millionth hex digit without knowing any of the earlier ones.
The algorithm is genuinely defined in terms of floating point fractional parts,
so it is the single place in the chapter where exact integer arithmetic is set
aside, by design.
Three views of e
The constant has three classical definitions, and the program implements
all of them. The exponential series converges at a factorial rate,
the limit definition 
converges only
linearly, and Euler’s continued fraction
produces a sequence of exact rational convergents, each buying roughly one more
digit. The program also keeps the partial sums of the series as exact rationals,
so you can see as a single enormous fraction.
The Program at a Glance
The code lives in math.rkt. Rather than a collection of loose scripts, the
program is built around a method registry: each algorithm is described by a
method structure carrying its id, title, formula, convergence notes, input
specifications (with defaults and caps), and a function that runs it. The
interactive menu, the command line interface, the comparison report, and the
test suite all work off this single table, so adding a fourteenth or fifteenth
algorithm means adding one entry, not touching any user interface code.
Every method takes small integer inputs given as key=value pairs, for example
digits=1000 or terms=500000. Here is the method table the program prints,
which is also the vocabulary for everything that follows:
1 # id method input
2 -- ---------------- -------------------------------------- -------------
3 1 leibniz Gregory-Leibniz series terms
4 2 nilakantha Nilakantha series terms
5 3 wallis Wallis product terms
6 4 viete Viete's nested radicals terms
7 5 archimedes Archimedes polygon doubling iterations
8 6 machin Machin-like arctangent formula digits
9 7 agm Gauss-Legendre / Brent-Salamin AGM digits
10 8 chudnovsky Chudnovsky binary splitting digits
11 9 bbp BBP hex-digit spigot start,count
12
13 10 e-series Exponential series to N digits digits
14 11 e-series-terms Exponential series, first n terms terms
15 12 e-series-exact Exponential series as an exact fraction terms
16 13 e-limit Limit definition (1 + 1/n)^n n
17 14 e-cf Continued fraction convergents terms
The rest of this section walks through the complete listing of math.rkt in
seven parts, following the section banners already in the file.
Part 1: File Header and Exports
The file opens with a usage comment that doubles as the manual for the command
line interface, followed by the provide list. Notice that the exports include
both the raw algorithm functions (pi-leibniz, pi-chudnovsky, and friends)
and the registry machinery (method, run-method, method-lookup), which is
what allows the test file to exercise everything through the same interface the
user sees.
1 #lang racket
2
3 ;; math.rkt — classical algorithms for pi and e, with a console UI.
4 ;;
5 ;; racket math.rkt interactive menu
6 ;; racket math.rkt --list show the method table
7 ;; racket math.rkt --compare run every method once
8 ;; racket math.rkt chudnovsky digits=1000 run one method
9 ;; racket math.rkt --show 40 leibniz terms=500000 (flags precede the id)
10 ;;
11 ;; Every algorithm works in exact integer arithmetic. A real number x is
12 ;; represented as the pair (v, s) meaning v / 10^s, so results are exact to
13 ;; the requested number of digits with no floating point drift. The one
14 ;; exception is the BBP spigot, which is defined in terms of flonum
15 ;; fractional parts.
16
17 (require racket/cmdline
18 racket/format)
19
20 (provide (struct-out fp)
21 fp->decimal
22 fp->digits
23 fp-display-str
24 fp->hex-frac-digits
25 count-matching-digits
26 reference-value
27 pi-leibniz pi-nilakantha pi-wallis pi-viete pi-archimedes
28 pi-machin pi-agm pi-chudnovsky pi-bbp-hex
29 e-series e-series-terms e-series-exact e-limit e-continued-fraction
30 (struct-out input-spec)
31 (struct-out method)
32 (struct-out outcome)
33 format-number-lines
34 show-limit
35 methods
36 method-lookup
37 run-method)
Part 2: The Fixed-Point Number Helpers
This is the foundation everything else stands on. An fp value packages an
exact integer (or, for the continued fraction method, an exact rational) with a
decimal scale. The central function is fp->decimal, which rescales a value to
a requested number of fractional digits and splits it into integer and fraction
parts. Two details are worth studying. First, digits are truncated, not
rounded, because the goal is to match published digit tables exactly; a
rounded last digit might differ from the table even when the computation is
correct. Second, the rescaling happens before the numerator and denominator are
taken, which matters when v is an exact rational: normalizing a rational
changes its denominator, so splitting first and scaling second would give wrong
answers. The comment in the code records this lesson.
count-matching-digits compares two digit strings and reports how many leading
characters agree; it is the workhorse of the verification reports and the test
suite. fp->hex-frac-digits converts the fractional part of a fixed-point value
to hexadecimal digits, which the test suite uses to cross-check the BBP spigot
against the decimal expansion.
1 ;; =====================================================================
2 ;; Section 1: numeric helpers
3 ;; =====================================================================
4
5 ;; A fixed-point value: the exact number `v` scaled by 10^-`s`.
6 ;; `v` is normally an exact integer; the continued fraction method returns
7 ;; an exact rational with s = 0.
8 (struct fp (v s) #:transparent)
9
10 (define (pow10 n) (expt 10 n))
11
12 (define (log10r x)
13 (if (<= x 1) 0.0 (/ (log x) (log 10))))
14
15 (define (log2r x)
16 (if (<= x 1) 0.0 (/ (log x) (log 2))))
17
18 (define (left-pad str n ch)
19 (if (>= (string-length str) n)
20 str
21 (string-append (make-string (- n (string-length str)) ch) str)))
22
23 ;; Round an exact rational n/d (both positive) to the nearest integer.
24 (define (round-quotient n d)
25 (quotient (+ (* 2 n) d) (* 2 d)))
26
27 ;; Split x into integer part and `digits` fractional digits, returning
28 ;; (values int-str frac-str digits-str). The digits are truncated, not
29 ;; rounded, so the output matches published digit tables of pi and e exactly.
30 (define (fp->decimal x digits)
31 (define e (- (fp-s x) digits))
32 ;; Rescale first, then split: normalising the rational changes the
33 ;; denominator, so numerator/denominator must be taken after scaling.
34 (define scaled (if (> e 0)
35 (/ (fp-v x) (pow10 e))
36 (* (fp-v x) (pow10 (- e)))))
37 (define total (quotient (numerator scaled) (denominator scaled)))
38 (define divisor (pow10 digits))
39 (define int-part (quotient total divisor))
40 (define frac-part (remainder total divisor))
41 (define int-str (number->string int-part))
42 (define frac-str (left-pad (number->string frac-part) digits #\0))
43 (values int-str frac-str (string-append int-str frac-str)))
44
45 ;; Just the digit string, without the integer/fraction split.
46 (define (fp->digits x digits)
47 (define-values (_i _f all) (fp->decimal x digits))
48 all)
49
50 (define (fp-display-str x digits)
51 (define-values (int-str frac-str _a) (fp->decimal x digits))
52 (if (zero? digits) int-str (string-append int-str "." frac-str)))
53
54 ;; Number of leading digits shared by two digit strings.
55 (define (count-matching-digits a b)
56 (define n (min (string-length a) (string-length b)))
57 (let loop ([i 0])
58 (cond [(= i n) n]
59 [(char=? (string-ref a i) (string-ref b i)) (loop (add1 i))]
60 [else i])))
61
62 ;; Hexadecimal digits of the fractional part of a fixed-point value.
63 (define (fp->hex-frac-digits x count)
64 (define den (pow10 (fp-s x)))
65 (define-values (_q rem0) (quotient/remainder (fp-v x) den))
66 (for/fold ([acc '()] [rem rem0] #:result (list->string (reverse acc)))
67 ([_i (in-range count)])
68 (define r (* rem 16))
69 (define-values (d next) (quotient/remainder r den))
70 (values (cons (string-ref "0123456789abcdef" d) acc) next)))
Part 3: The Pi Algorithms
Here is the first remarkable thing about the code: every algorithm is a direct
transliteration of its formula into integer arithmetic. Look at pi-leibniz.
The accumulator is a plain integer, each term is S / (2k + 1) where
, and the sign alternates with (even? k). There is no floating
point anywhere, so the only error is the truncation of the series itself, the
error we predicted on paper.
Two idioms deserve attention. pi-viete and pi-archimedes need square roots
at full precision, and Racket’s integer-sqrt gives exactly that: the floor of
the true square root of an exact integer. And pi-archimedes returns two
values, the lower and upper bounds, using Racket’s multiple return values; the
caller in the method registry uses both to report how many digits the bracket
pins down.
The historical progression is visible in the function names themselves: series from the 14th and 17th centuries, Viete’s product from 1593, and Archimedes’ doubling that is over two thousand years old, all expressed in the same twenty lines of style.
1 ;; =====================================================================
2 ;; Section 2: pi
3 ;; =====================================================================
4
5 ;; pi/4 = 1 - 1/3 + 1/5 - 1/7 + ... (Gregory-Leibniz, 1350s/1670s)
6 ;; Error after n terms is about 1/(2n): painfully linear convergence.
7 (define (pi-leibniz terms s)
8 (define S (pow10 s))
9 (define acc
10 (for/fold ([acc 0]) ([k (in-range terms)])
11 (define term (quotient S (+ (* 2 k) 1)))
12 (if (even? k) (+ acc term) (- acc term))))
13 (* 4 acc))
14
15 ;; pi = 3 + 4/(2*3*4) - 4/(4*5*6) + 4/(6*7*8) - ... (Nilakantha, 14th c.)
16 ;; Error after n terms is about 1/(4n^3): three digits per tenfold of terms.
17 (define (pi-nilakantha terms s)
18 (define S (pow10 s))
19 (for/fold ([acc (* 3 S)]) ([k (in-range 1 (add1 terms))])
20 (define d (* (* 2 k) (+ (* 2 k) 1) (+ (* 2 k) 2)))
21 (define term (quotient (* 4 S) d))
22 (if (odd? k) (+ acc term) (- acc term))))
23
24 ;; pi/2 = (2*2)/(1*3) * (4*4)/(3*5) * (6*6)/(5*7) * ... (Wallis, 1656)
25 (define (pi-wallis terms s)
26 (define S (pow10 s))
27 (define prod
28 (for/fold ([p S]) ([n (in-range 1 (add1 terms))])
29 (define two-n (* 2 n))
30 (define factor (round-quotient (* two-n two-n S)
31 (* (sub1 two-n) (add1 two-n))))
32 (quotient (* p factor) S)))
33 (* 2 prod))
34
35 ;; 2/pi = sqrt(1/2) * sqrt(1/2 + 1/2*sqrt(1/2)) * ... (Viete, 1593)
36 ;; Each c_k = sqrt((1 + c_(k-1))/2), starting from c_0 = 0.
37 (define (pi-viete terms s)
38 (define S (pow10 s))
39 (define prod
40 (for/fold ([c 0] [p S] #:result p) ([_k (in-range terms)])
41 (define c2 (integer-sqrt (* S (quotient (+ S c) 2))))
42 (values c2 (quotient (* p c2) S))))
43 (quotient (* 2 S S) prod))
44
45 ;; Archimedes' polygon doubling, written as the Pfaff/Borchardt pair of
46 ;; means: a <- 2ab/(a+b) (harmonic, circumscribed), b <- sqrt(a*b)
47 ;; (geometric, inscribed). Both converge to pi and always bracket it,
48 ;; gaining about 0.6 digits per doubling. Returns (values lower upper).
49 (define (pi-archimedes iterations s)
50 (define S (pow10 s))
51 (define a0 (integer-sqrt (* 12 S S))) ; 2*sqrt(3), circumscribed hexagon
52 (define b0 (* 3 S)) ; 3, inscribed hexagon
53 (let loop ([a a0] [b b0] [k iterations])
54 (if (zero? k)
55 (values b a)
56 (let* ([a2 (quotient (* 2 a b) (+ a b))]
57 [b2 (integer-sqrt (* a2 b))])
58 (loop a2 b2 (sub1 k))))))
59
60 ;; arctan(1/x) = sum (-1)^k / ((2k+1) x^(2k+1)), in fixed point.
61 (define (arctan-recip x s)
62 (define S (pow10 s))
63 (define x2 (* x x))
64 (let loop ([k 0] [num (quotient S x)] [acc 0])
65 (if (zero? num)
66 acc
67 (let ([term (quotient num (+ (* 2 k) 1))])
68 (loop (add1 k)
69 (quotient num x2)
70 (if (even? k) (+ acc term) (- acc term)))))))
71
72 ;; Machin's formula (1706): pi = 16*arctan(1/5) - 4*arctan(1/239).
73 ;; About 1.4 digits per term of the arctan series.
74 (define (pi-machin s)
75 (- (* 16 (arctan-recip 5 s))
76 (* 4 (arctan-recip 239 s))))
77
78 ;; Gauss-Legendre / Brent-Salamin AGM iteration (1975-76).
79 ;; Digits roughly double at every step.
80 (define (pi-agm s)
81 (define S (pow10 s))
82 (define iterations (+ 1 (exact-ceiling (log2r s))))
83 (let loop ([a S]
84 [b (integer-sqrt (quotient (* S S) 2))] ; 1/sqrt(2)
85 [t (quotient S 4)]
86 [p 1]
87 [k iterations])
88 (if (zero? k)
89 (quotient (* (+ a b) (+ a b)) (* 4 t))
90 (let* ([a2 (quotient (+ a b) 2)]
91 [b2 (integer-sqrt (* a b))]
92 [d (- a a2)]
93 [t2 (- t (quotient (* p (* d d)) S))]
94 [p2 (* 2 p)])
95 (loop a2 b2 t2 p2 (sub1 k))))))
96
97 (define CHUD-Q 10939058860032000) ; 640320^3 / 24
98
99 ;; Binary splitting of the Chudnovsky series (1988):
100 ;; 1/pi = 12 * sum (-1)^k (6k)!(13591409 + 545140134k) / ((3k)!(k!)^3 640320^(3k+3/2))
101 ;; Each term buys about 14.18 digits, which is why it holds the records.
102 (define (chudnovsky-split a b)
103 (if (= 1 (- b a))
104 (if (zero? a)
105 (values 1 1 13591409)
106 (let ([P (* (- (* 6 a) 5) (- (* 2 a) 1) (- (* 6 a) 1))])
107 (values P
108 (* (expt a 3) CHUD-Q)
109 (* P (+ 13591409 (* 545140134 a)) (if (even? a) 1 -1)))))
110 (let ([m (quotient (+ a b) 2)])
111 (define-values (Pl Ql Tl) (chudnovsky-split a m))
112 (define-values (Pr Qr Tr) (chudnovsky-split m b))
113 (values (* Pl Pr)
114 (* Ql Qr)
115 (+ (* Tl Qr) (* Pl Tr))))))
116
117 (define (pi-chudnovsky s)
118 (define terms (+ 2 (quotient s 14)))
119 (define-values (_P Q T) (chudnovsky-split 0 terms))
120 (define root (integer-sqrt (* 10005 (pow10 s) (pow10 s)))) ; sqrt(10005)*10^s
121 (quotient (* Q 426880 root) T))
122
123 ;; Bailey-Borwein-Plouffe (1995): hex digit `d` of pi without any of the
124 ;; preceding digits.
125 ;; pi = sum 16^-k (4/(8k+1) - 2/(8k+4) - 1/(8k+5) - 1/(8k+6))
126 ;; This one genuinely needs flonum fractional parts.
127 (define (modpow16 e m)
128 (let loop ([e e] [b (modulo 16 m)] [acc 1])
129 (cond [(zero? e) acc]
130 [(even? e) (loop (quotient e 2) (modulo (* b b) m) acc)]
131 [else (loop (sub1 e) b (modulo (* acc b) m))])))
132
133 (define (bbp-series j d)
134 (define head
135 (for/fold ([acc 0.0]) ([k (in-range d)])
136 (define m (+ (* 8 k) j))
137 (define term (/ (exact->inexact (modpow16 (- d k 1) m))
138 (exact->inexact m)))
139 (define x (+ acc term))
140 (- x (floor x))))
141 (define tail
142 (for/fold ([acc 0.0]) ([k (in-range d (+ (* 2 d) 10))])
143 (+ acc (/ (expt 16.0 (- d k 1)) (exact->inexact (+ (* 8 k) j))))))
144 (define x (+ head tail))
145 (- x (floor x)))
146
147 (define (pi-bbp-hex-digit d)
148 (define x (- (* 4.0 (bbp-series 1 d))
149 (* 2.0 (bbp-series 4 d))
150 (bbp-series 5 d)
151 (bbp-series 6 d)))
152 (define frac (- x (floor x)))
153 (string (string-ref "0123456789ABCDEF"
154 (inexact->exact (floor (* 16.0 frac))))))
155
156 ;; Hex digits of pi starting at 1-based position `start`.
157 (define (pi-bbp-hex start count)
158 (apply string-append
159 (for/list ([d (in-range start (+ start count))])
160 (pi-bbp-hex-digit d))))
The modern algorithms show why the registry’s metadata about convergence
matters. pi-machin loops until the current numerator is zero at the working
scale, so it stops exactly when further terms can no longer change the result.
pi-agm computes its iteration count from 
of the target scale,
because quadratic convergence means the iteration count grows only
logarithmically with the number of digits wanted. chudnovsky-split is the
most sophisticated recursion in the file: it returns three values 
, 
,
and 
for a range of terms and combines subranges with the identities
, 
, 
, turning an
-term sum into an 
deep product tree.
The BBP section is the one place flonums appear. modpow16 computes
by repeated squaring, keeping numbers small, and bbp-series
accumulates only the fractional part of each partial sum. The trade-off is
explicit: we gave up exactness, and in return we can teleport to digit .
Part 4: The Algorithms for e
The three classical views of translate just as directly. e-series
accumulates 
by dividing the previous term by 
and stops when the
next term is zero at the working scale. e-series-exact is more interesting:
it keeps the partial sum as one exact rational with denominator 
,
built with pure integer arithmetic so that Racket never has to normalize the
fraction by a greatest common divisor computation. e-limit evaluates
with binary exponentiation on fixed-point values, and
e-continued-fraction runs the standard convergent recurrence
over Euler’s pattern of partial quotients
.
1 ;; =====================================================================
2 ;; Section 3: e
3 ;; =====================================================================
4
5 ;; e = sum 1/k!, run until the next term no longer affects the scale.
6 (define (e-series s)
7 (define S (pow10 s))
8 (let loop ([k 1] [term S] [acc S])
9 (define next (quotient term k))
10 (if (zero? next) acc (loop (add1 k) next (+ acc next)))))
11
12 ;; The same series truncated after exactly `terms` terms (k = 0 .. terms-1).
13 (define (e-series-terms terms s)
14 (define S (pow10 s))
15 (let loop ([k 1] [term S] [acc S])
16 (if (>= k terms)
17 acc
18 (let ([next (quotient term k)])
19 (loop (add1 k) next (+ acc next))))))
20
21 ;; The exact rational partial sum sum_{k=0}^{terms-1} 1/k!, built with
22 ;; integer arithmetic only so no gcd normalisation blowup occurs.
23 (define (e-series-exact terms)
24 (define n (max 1 terms))
25 (define Q (for/fold ([f 1]) ([i (in-range 2 n)]) (* f i))) ; (n-1)!
26 (define T
27 (let loop ([k 0] [m Q] [acc 0])
28 (if (= k n)
29 acc
30 (loop (add1 k) (quotient m (add1 k)) (+ acc m)))))
31 (/ T Q))
32
33 ;; (1 + 1/n)^n by binary exponentiation on fixed-point values.
34 (define (fp-pow base exponent s)
35 (define S (pow10 s))
36 (let loop ([e exponent] [b base] [acc S])
37 (cond [(zero? e) acc]
38 [(even? e) (loop (quotient e 2) (quotient (* b b) S) acc)]
39 [else (loop (sub1 e) b (quotient (* acc b) S))])))
40
41 (define (e-limit n s)
42 (define S (pow10 s))
43 (define base (+ S (quotient S (max 1 n))))
44 (fp-pow base (max 1 n) s))
45
46 ;; e = [2; 1,2,1, 1,4,1, 1,6,1, 1,8,1, ...] (Euler, 1737)
47 ;; Returns the exact rational of the `terms`-th convergent.
48 (define (e-cf-partial i)
49 (cond [(= i 0) 2]
50 [(= 2 (modulo i 3)) (* 2 (quotient (add1 i) 3))]
51 [else 1]))
52
53 (define (e-continued-fraction terms)
54 (define n (max 1 terms))
55 (define-values (p q)
56 (for/fold ([p-prev2 0] [p-prev1 1]
57 [q-prev2 1] [q-prev1 0]
58 #:result (values p-prev1 q-prev1))
59 ([i (in-range n)])
60 (define a (e-cf-partial i))
61 (values p-prev1 (+ (* a p-prev1) p-prev2)
62 q-prev1 (+ (* a q-prev1) q-prev2))))
63 (/ p q))
Note the asymmetry between the two families: the e functions take the scale
s as an argument, exactly like the functions, except
e-series-exact and e-continued-fraction, which return exact rationals
instead of scaled integers. That small difference is why the fp struct allows
v to be a rational with 
: the formatting machinery can treat both
representations uniformly.
Part 5: Reference Values and the Registry Scaffolding
The verification story starts here. reference-value computes a trusted high
precision value of either constant, using the Chudnovsky series for and
the exponential series for , and caches it in a hash table keyed by
constant and digit count. Every report and every test that asks “how many
digits are correct” is really asking “how many leading digits agree with this
independently computed reference.”
Below it are the three structures that define the little framework: input-spec
describes one named integer input with a label, default, and cap; method
bundles everything the user interface needs to know about an algorithm;
outcome is what running a method produces, either a fixed-point value plus
formatted lines, or, for the BBP spigot, a string of hex digits. The
estimate-digits helper implements the rule of thumb mentioned in the theory
section: if a method gains rate digits per tenfold increase in its input, the
suggested display precision is 
.
1 ;; =====================================================================
2 ;; Section 4: reference values and verification
3 ;; =====================================================================
4
5 (define reference-cache (make-hash))
6
7 ;; A trusted high-precision value, used to say how many digits a method
8 ;; actually got right.
9 (define (reference-value constant digits)
10 (hash-ref! reference-cache
11 (list constant digits)
12 (lambda ()
13 (define s (+ digits 10))
14 (fp (case constant
15 [(pi) (pi-chudnovsky s)]
16 [(e) (e-series s)])
17 s))))
18
19 ;; =====================================================================
20 ;; Section 5: the method table
21 ;; =====================================================================
22
23 (struct input-spec (key label default cap) #:transparent)
24 (struct method (id
25 title
26 constant
27 inputs ; (listof input-spec)
28 auto-digits ; (-> hash? natural?) suggested display digits
29 run ; (-> hash? natural? outcome?)
30 formula
31 convergence
32 compare-inputs)
33 #:transparent)
34 (struct outcome (value lines hex) #:transparent)
35
36 (define GUARD 25)
37
38 ;; Suggest how many digits to display, from the expected truncation error:
39 ;; `rate` digits gained per tenfold increase in the input, plus `offset`.
40 (define (estimate-digits input rate offset)
41 (max 1 (min 200 (inexact->exact
42 (floor (+ offset (* rate (log10r (max 2 input)))))))))
43
44 ;; Build the outcome for an ordinary fixed-point result.
45 (define (fp-outcome v digits show)
46 (define-values (int-str frac-str _all) (fp->decimal v digits))
47 (outcome v (format-number-lines int-str frac-str show) #f))
48
49 (define (group10 str)
50 (define n (string-length str))
51 (for/list ([i (in-range 0 n 10)])
52 (substring str i (min n (+ i 10)))))
53
54 (define (format-number-lines int-str frac-str show)
55 (define total (string-length frac-str))
56 (define shown (min total (max 0 show)))
57 (define groups (group10 (substring frac-str 0 shown)))
58 (define per-line 5)
59 (define chunks
60 (let loop ([g groups])
61 (cond [(null? g) '()]
62 [(<= (length g) per-line) (list g)]
63 [else (cons (take g per-line) (loop (drop g per-line)))])))
64 (define body
65 (if (null? chunks)
66 (list (string-append " " int-str (if (zero? total) "" ".")))
67 (cons (string-append " " int-str "." (string-join (car chunks) " "))
68 (for/list ([c (in-list (cdr chunks))])
69 (string-append " " (string-join c " "))))))
70 (append body
71 (if (> total shown)
72 (list (format " ... showing ~a of ~a digits; last 10: ~a"
73 shown total
74 (string-join (take-right (group10 frac-str) 1) " ")))
75 '())))
76
77 (define (short-num str limit)
78 (define n (string-length str))
79 (if (<= n limit) str (string-append (substring str 0 limit) "...")))
80
81 (define (exact-rational-lines r)
82 (list (format " exact rational: ~a / ~a"
83 (short-num (number->string (numerator r)) 40)
84 (short-num (number->string (denominator r)) 40))
85 (format " (~a-digit numerator, ~a-digit denominator)"
86 (string-length (number->string (numerator r)))
87 (string-length (number->string (denominator r))))))
The GUARD constant of 25 extra digits appears in every method’s run
function. Working at higher precision than you display is the standard defense
against truncation artifacts in the intermediate arithmetic: the guard digits
absorb the accumulated integer divisions so that the displayed digits are clean.
Part 6: The Method Registry
Each entry pairs the raw algorithm with everything the user interface needs: the
formula as a display string, a plain language convergence note, input defaults
and caps, a suggested digit count, and the inputs used by the --compare
report. Read the leibniz entry closely and the pattern for all fourteen is
clear: its scale request adds guard digits plus 
of the term count,
because dividing 
by numbers up to 
loses about that many digits of
headroom.
1 ;; ---- pi methods -----------------------------------------------------
2
3 (define m-leibniz
4 (method 'leibniz
5 "Gregory-Leibniz series"
6 'pi
7 (list (input-spec 'terms "series terms" 100000 10000000))
8 (lambda (in) (estimate-digits (hash-ref in 'terms) 1 -1))
9 (lambda (in digits)
10 (define s (+ digits GUARD (exact-ceiling (log10r (max 2 (hash-ref in 'terms))))))
11 (fp-outcome (fp (pi-leibniz (hash-ref in 'terms) s) s) digits (show-limit)))
12 "pi/4 = 1 - 1/3 + 1/5 - 1/7 + 1/9 - ..."
13 "linear: error ~ 1/(2n), so one extra digit costs ten times the terms"
14 (hash 'terms 200000)))
15
16 (define m-nilakantha
17 (method 'nilakantha
18 "Nilakantha series"
19 'pi
20 (list (input-spec 'terms "series terms" 1000 1000000))
21 (lambda (in) (estimate-digits (hash-ref in 'terms) 3 0))
22 (lambda (in digits)
23 (define s (+ digits GUARD 10))
24 (fp-outcome (fp (pi-nilakantha (hash-ref in 'terms) s) s) digits (show-limit)))
25 "pi = 3 + 4/(2*3*4) - 4/(4*5*6) + 4/(6*7*8) - ..."
26 "cubic: error ~ 1/(4n^3), so three extra digits per tenfold of terms"
27 (hash 'terms 5000)))
28
29 (define m-wallis
30 (method 'wallis
31 "Wallis product"
32 'pi
33 (list (input-spec 'terms "product factors" 10000 1000000))
34 (lambda (in) (estimate-digits (hash-ref in 'terms) 1 -1))
35 (lambda (in digits)
36 (define s (+ digits GUARD (exact-ceiling (log10r (max 2 (hash-ref in 'terms))))))
37 (fp-outcome (fp (pi-wallis (hash-ref in 'terms) s) s) digits (show-limit)))
38 "pi/2 = (2*2)/(1*3) * (4*4)/(3*5) * (6*6)/(5*7) * ..."
39 "linear: error ~ pi/(4n), so one extra digit costs ten times the factors"
40 (hash 'terms 50000)))
41
42 (define m-viete
43 (method 'viete
44 "Viete's nested radicals"
45 'pi
46 (list (input-spec 'terms "radical factors" 30 2000))
47 (lambda (in) (min 200 (max 1 (inexact->exact (floor (* 0.6 (hash-ref in 'terms)))))))
48 (lambda (in digits)
49 (define s (+ digits GUARD))
50 (fp-outcome (fp (pi-viete (hash-ref in 'terms) s) s) digits (show-limit)))
51 "2/pi = sqrt(1/2) * sqrt(1/2 + 1/2*sqrt(1/2)) * ..."
52 "linear: about 0.6 digits per factor (the first infinite product for pi, 1593)"
53 (hash 'terms 40)))
54
55 (define m-archimedes
56 (method 'archimedes
57 "Archimedes polygon doubling"
58 'pi
59 (list (input-spec 'iterations "polygon doublings" 20 5000))
60 (lambda (in) (min 200 (max 1 (inexact->exact (floor (* 0.6 (hash-ref in 'iterations)))))))
61 (lambda (in digits)
62 (define s (+ digits GUARD))
63 (define iters (hash-ref in 'iterations))
64 (define sides-desc
65 (if (<= iters 40)
66 (format "~a-gon" (* 6 (expt 2 iters)))
67 (format "6*2^~a-gon" iters)))
68 (define-values (lo hi) (pi-archimedes iters s))
69 (define lower (fp lo s))
70 (define upper (fp hi s))
71 (define-values (int-str frac-str _a) (fp->decimal lower digits))
72 (define-values (_i2 frac-str2 _a2) (fp->decimal upper digits))
73 (define pinned
74 (max 0 (sub1 (count-matching-digits (string-append int-str frac-str)
75 (string-append int-str frac-str2)))))
76 (outcome lower
77 (append (list (format " inscribed ~a (lower bound):" sides-desc))
78 (format-number-lines int-str frac-str 60)
79 (list (format " circumscribed ~a (upper bound): ~a"
80 sides-desc
81 (fp-display-str upper digits))
82 (format " digits pinned by the bracket: ~a" pinned)))
83 #f))
84 "a <- 2ab/(a+b); b <- sqrt(a*b), from the hexagon a=2*sqrt(3), b=3"
85 "linear: about 0.6 digits per doubling, but every step brackets pi exactly"
86 (hash 'iterations 30)))
87
88 (define m-machin
89 (method 'machin
90 "Machin-like arctangent formula"
91 'pi
92 (list (input-spec 'digits "decimal digits" 100 1000000))
93 (lambda (in) (hash-ref in 'digits))
94 (lambda (in digits)
95 (define s (+ digits GUARD))
96 (fp-outcome (fp (pi-machin s) s) digits (show-limit)))
97 "pi = 16*arctan(1/5) - 4*arctan(1/239)"
98 "about 1.4 digits per arctan term; the classic hand-computation formula"
99 (hash 'digits 50)))
100
101 (define m-agm
102 (method 'agm
103 "Gauss-Legendre / Brent-Salamin AGM"
104 'pi
105 (list (input-spec 'digits "decimal digits" 100 1000000))
106 (lambda (in) (hash-ref in 'digits))
107 (lambda (in digits)
108 (define s (+ digits GUARD))
109 (fp-outcome (fp (pi-agm s) s) digits (show-limit)))
110 "a,b <- (a+b)/2, sqrt(ab); t <- t - p(a-a')^2; pi ~ (a+b)^2/(4t)"
111 "quadratic: the number of correct digits roughly doubles each iteration"
112 (hash 'digits 50)))
113
114 (define m-chudnovsky
115 (method 'chudnovsky
116 "Chudnovsky binary splitting"
117 'pi
118 (list (input-spec 'digits "decimal digits" 1000 2000000))
119 (lambda (in) (hash-ref in 'digits))
120 (lambda (in digits)
121 (define s (+ digits GUARD))
122 (fp-outcome (fp (pi-chudnovsky s) s) digits (show-limit)))
123 "1/pi = 12 * sum (-1)^k (6k)!(13591409+545140134k) / ((3k)!(k!)^3 640320^(3k+3/2))"
124 "about 14.18 digits per term; used for every modern pi record"
125 (hash 'digits 50)))
126
127 (define m-bbp
128 (method 'bbp
129 "BBP hex-digit spigot"
130 'pi
131 (list (input-spec 'start "starting hex digit position (1-based)" 1 200000)
132 (input-spec 'count "how many hex digits" 40 1000))
133 (lambda (in) 0)
134 (lambda (in digits)
135 (define start (hash-ref in 'start))
136 (define count (hash-ref in 'count))
137 (define hx (pi-bbp-hex start count))
138 (outcome #f
139 (list (format " hex digits ~a..~a of pi, after the decimal point:"
140 start (+ start count -1))
141 (if (= 1 start)
142 (format " 3.~a" hx)
143 (format " ...~a" hx)))
144 hx))
145 "pi = sum 16^-k (4/(8k+1) - 2/(8k+4) - 1/(8k+5) - 1/(8k+6))"
146 "a spigot: digit d in base 16 with none of the earlier digits computed"
147 (hash 'start 1 'count 24)))
148
149 ;; ---- e methods ------------------------------------------------------
150
151 (define m-e-series
152 (method 'e-series
153 "Exponential series to N digits"
154 'e
155 (list (input-spec 'digits "decimal digits" 100 1000000))
156 (lambda (in) (hash-ref in 'digits))
157 (lambda (in digits)
158 (define s (+ digits GUARD))
159 (fp-outcome (fp (e-series s) s) digits (show-limit)))
160 "e = sum_{k>=0} 1/k! = 1 + 1 + 1/2 + 1/6 + 1/24 + ..."
161 "factorial convergence: term k is worth about k*log10(k/e) digits"
162 (hash 'digits 50)))
163
164 (define m-e-series-terms
165 (method 'e-series-terms
166 "Exponential series, first n terms"
167 'e
168 (list (input-spec 'terms "series terms" 40 200000))
169 (lambda (in)
170 (define n (max 2 (hash-ref in 'terms)))
171 (min 200 (max 1 (inexact->exact
172 (floor (+ 0.5 (- (* n (log10r n)) (* n (log10r (exp 1))))))))))
173 (lambda (in digits)
174 (define n (hash-ref in 'terms))
175 (define s (+ digits GUARD 10))
176 (fp-outcome (fp (e-series-terms n s) s) digits (show-limit)))
177 "e ~ sum_{k=0}^{n-1} 1/k!"
178 "error is the tail 1/n! ~ about n*log10(n/e) digits from n terms"
179 (hash 'terms 40)))
180
181 (define m-e-series-exact
182 (method 'e-series-exact
183 "Exponential series as an exact fraction"
184 'e
185 (list (input-spec 'terms "series terms" 25 20000))
186 (lambda (in)
187 (define n (max 2 (hash-ref in 'terms)))
188 (min 200 (max 1 (inexact->exact
189 (floor (+ 0.5 (- (* n (log10r n)) (* n (log10r (exp 1))))))))))
190 (lambda (in digits)
191 (define r (e-series-exact (hash-ref in 'terms)))
192 (define v (fp r 0))
193 (define-values (int-str frac-str _a) (fp->decimal v digits))
194 (outcome v
195 (append (format-number-lines int-str frac-str (show-limit))
196 (exact-rational-lines r))
197 #f))
198 "e ~ sum_{k=0}^{n-1} 1/k!, kept as one exact rational"
199 "same series as above, but the answer is an exact fraction rather than a decimal"
200 (hash 'terms 25)))
201
202 (define m-e-limit
203 (method 'e-limit
204 "Limit definition (1 + 1/n)^n"
205 'e
206 (list (input-spec 'n "n" 100000 10000000))
207 (lambda (in) (estimate-digits (hash-ref in 'n) 1 -1))
208 (lambda (in digits)
209 (define n (hash-ref in 'n))
210 (define s (+ digits GUARD (exact-ceiling (log2r (max 2 n)))))
211 (fp-outcome (fp (e-limit n s) s) digits (show-limit)))
212 "e = lim (1 + 1/n)^n"
213 "linear: error ~ e/(2n), one extra digit per tenfold of n"
214 (hash 'n 100000)))
215
216 (define m-e-cf
217 (method 'e-cf
218 "Continued fraction convergents"
219 'e
220 (list (input-spec 'terms "partial quotients" 40 100000))
221 (lambda (in) (min 200 (max 1 (quotient (* 9 (hash-ref in 'terms)) 10))))
222 (lambda (in digits)
223 (define r (e-continued-fraction (hash-ref in 'terms)))
224 (define v (fp r 0))
225 (define-values (int-str frac-str _a) (fp->decimal v digits))
226 (define partials
227 (for/list ([i (in-range (min 24 (hash-ref in 'terms)))])
228 (number->string (e-cf-partial i))))
229 (define tail (string-join (cdr partials) ", "))
230 (outcome v
231 (append (format-number-lines int-str frac-str (show-limit))
232 (list (format " convergent: [~a~a~a~a]"
233 (car partials)
234 (if (string=? tail "") "" "; ")
235 tail
236 (if (>= (hash-ref in 'terms) 24) ", ..." ""))
237 (format " from ~a partial quotients"
238 (hash-ref in 'terms)))
239 (exact-rational-lines r))
240 #f))
241 "e = [2; 1,2,1, 1,4,1, 1,6,1, 1,8,1, ...]"
242 "the partial quotients repeat 1,2k,1; each convergent buys about one digit"
243 (hash 'terms 40)))
244
245 (define methods
246 (list m-leibniz m-nilakantha m-wallis m-viete m-archimedes
247 m-machin m-agm m-chudnovsky m-bbp
248 m-e-series m-e-series-terms m-e-series-exact m-e-limit m-e-cf))
249
250 (define (method-lookup id-or-index)
251 (define n (and (string? id-or-index) (string->number id-or-index)))
252 (cond [(and (exact-integer? n) (<= 1 n (length methods)))
253 (list-ref methods (sub1 n))]
254 [(string? id-or-index)
255 (define sym (string->symbol id-or-index))
256 (findf (lambda (m) (eq? (method-id m) sym)) methods)]
257 [else #f]))
258
259 (define (run-method m inputs digits)
260 ((method-run m) inputs digits))
The archimedes entry shows why the registry stores functions rather than
precomputed data: its run lambda prints the polygon size (a 30 doubling run is
a polygon with 
sides), both bounds, and
the count of digits the two bounds agree on, which is information only
Archimedes’ method can certify. The bbp entry is the only one whose outcome
has #f for the value and a string in the hex field, a difference the
reporting code checks for explicitly.
Part 7: Reporting, the Console UI, and the Command Line
The final stretch of the file turns outcomes into printed reports and wires up
the two front ends. verification-lines is the interesting piece: for ordinary
methods it compares the result against the cached reference digit by digit and
reports either full agreement or the exact position where the digits diverge,
showing both strings around the divergence point. For the BBP outcome it
converts the reference value to hexadecimal and compares in base 16. Inputs to
the program are small and uniform: each method declares its input-spec list,
the command line accepts them as key=value pairs after the method id (or with
repeated --set key value flags), and the interactive loop prompts for them one
at a time with the defaults shown in brackets.
1 ;; =====================================================================
2 ;; Section 6: reporting and verification
3 ;; =====================================================================
4
5 (define show-limit (make-parameter 200))
6
7 (define (verification-lines m out digits inputs)
8 (define hex-str (outcome-hex out))
9 (cond
10 [hex-str
11 (define start (hash-ref inputs 'start 1))
12 (define count (string-length hex-str))
13 (cond
14 [(> (+ start count) 20000)
15 (list " check: skipped (position too deep for a full-precision reference)")]
16 [else
17 (define ref (reference-value 'pi (exact-ceiling (+ 20 (* 1.21 (+ start count))))))
18 (define ref-hex (fp->hex-frac-digits ref count))
19 (define match (count-matching-digits (string-downcase hex-str) ref-hex))
20 (if (= match count)
21 (list (format " check: all ~a hex digits match the reference pi" match))
22 (list (format " check: ~a of ~a hex digits match" match count)
23 (format " got: ~a" hex-str)
24 (format " ref: ~a" (string-upcase ref-hex))))])]
25 [(outcome-value out)
26 => (lambda (v)
27 (define ref (reference-value (method-constant m) (+ digits 10)))
28 (define got (fp->digits v digits))
29 (define exp (fp->digits ref digits))
30 (define match (count-matching-digits got exp))
31 (define correct (max 0 (sub1 match)))
32 (cond
33 [(>= match (string-length exp))
34 (list (format " check: all ~a displayed digits match the reference (~a at ~a digits)"
35 correct (method-id (reference-method (method-constant m)))
36 (+ digits 10)))]
37 [else
38 (define at (max 0 (sub1 match)))
39 (define lo (max 0 (- at 4)))
40 (define hi (min (string-length got) (+ at 12)))
41 (list (format " check: ~a correct decimal digits, diverges at digit ~a"
42 correct (add1 correct))
43 (format " got: ...~a|~a"
44 (substring got lo at)
45 (substring got at hi))
46 (format " ref: ...~a|~a"
47 (substring exp lo at)
48 (substring exp at hi))
49 (format " (reference: ~a at ~a digits)"
50 (method-id (reference-method (method-constant m)))
51 (+ digits 10)))]))]
52 [else '()]))
53
54 (define (reference-method constant)
55 (case constant
56 [(pi) m-chudnovsky]
57 [(e) m-e-series]))
58
59 (define (format-elapsed ms)
60 (cond [(< ms 1) (format "~a us" (exact-round (* ms 1000)))]
61 [(< ms 1000) (format "~a ms" (~r ms #:precision 2))]
62 [else (format "~a s" (~r (/ ms 1000) #:precision 2))]))
63
64 (define (print-report m out digits elapsed-ms inputs)
65 (define input-desc
66 (string-join (for/list ([spec (method-inputs m)])
67 (format "~a=~a" (input-spec-key spec)
68 (hash-ref inputs (input-spec-key spec))))
69 ", "))
70 (displayln (format " ~a" (method-title m)))
71 (displayln (format " ~a" (method-formula m)))
72 (displayln (format " ~a" (method-convergence m)))
73 (displayln (format " inputs: ~a" input-desc))
74 (for-each displayln (outcome-lines out))
75 (for-each displayln (verification-lines m out digits inputs))
76 (displayln (format " time: ~a" (format-elapsed elapsed-ms))))
77
78 ;; =====================================================================
79 ;; Section 7: console UI
80 ;; =====================================================================
81
82 (define (print-banner)
83 (displayln "=======================================================================")
84 (displayln " math.rkt - classical computations of pi and e")
85 (displayln "=======================================================================")
86 (displayln " Exact integer / fixed-point arithmetic throughout (the BBP spigot is")
87 (displayln " the one flonum algorithm, by design).")
88 (displayln ""))
89
90 (define (print-method-table)
91 (displayln " # id method input")
92 (displayln " -- ---------------- -------------------------------------- -------------")
93 (for ([(m i) (in-indexed methods)])
94 (when (and (> i 0)
95 (not (eq? (method-constant m) (method-constant (list-ref methods (sub1 i))))))
96 (displayln ""))
97 (printf " ~a ~a ~a ~a~n"
98 (left-pad (number->string (add1 i)) 2 #\space)
99 (~a (method-id m) #:min-width 16)
100 (~a (method-title m) #:min-width 38)
101 (string-join (for/list ([spec (method-inputs m)])
102 (format "~a" (input-spec-key spec)))
103 ","))))
104
105 (define (print-help m)
106 (displayln (format " ~a [~a]" (method-title m) (method-id m)))
107 (displayln (format " computes: ~a" (method-constant m)))
108 (displayln (format " formula: ~a" (method-formula m)))
109 (displayln (format " convergence: ~a" (method-convergence m)))
110 (displayln " inputs:")
111 (for ([spec (method-inputs m)])
112 (displayln (format " ~a: ~a (default ~a, max ~a)"
113 (input-spec-key spec)
114 (input-spec-label spec)
115 (input-spec-default spec)
116 (input-spec-cap spec)))))
117
118 ;; Read a line, or return eof.
119 (define (read-choice prompt-str)
120 (display prompt-str)
121 (flush-output)
122 (define line (read-line))
123 (if (eof-object? line) eof (string-downcase (string-trim line))))
124
125 (define (prompt-int label default cap)
126 (let loop ()
127 (define raw (read-choice (format " ~a [~a]: " label default)))
128 (cond
129 [(eof-object? raw) eof]
130 [(string=? raw "") default]
131 [else
132 (define n (string->number raw))
133 (cond
134 [(not (and (number? n) (exact-integer? n) (positive? n)))
135 (displayln " please enter a positive whole number (or enter for the default)")
136 (loop)]
137 [(> n cap)
138 (displayln (format " ~a is above the limit ~a - using ~a" n cap cap))
139 cap]
140 [else n])])))
141
142 (define (prompt-inputs m)
143 (let loop ([specs (method-inputs m)] [acc (hash)])
144 (cond [(null? specs) acc]
145 [else
146 (define spec (car specs))
147 (define v (prompt-int (input-spec-label spec)
148 (input-spec-default spec)
149 (input-spec-cap spec)))
150 (if (eof-object? v)
151 eof
152 (loop (cdr specs) (hash-set acc (input-spec-key spec) v)))])))
153
154 (define (display-digits-for m inputs)
155 (define from-inputs (hash-ref inputs 'digits #f))
156 (if from-inputs
157 from-inputs
158 (prompt-int "decimal digits to display"
159 ((method-auto-digits m) inputs)
160 1000000)))
161
162 (define (execute m inputs digits)
163 (define t0 (current-inexact-milliseconds))
164 (define out
165 (with-handlers ([exn:fail?
166 (lambda (e)
167 (printf " error: ~a~n" (exn-message e))
168 #f)])
169 (run-method m inputs digits)))
170 (define elapsed (- (current-inexact-milliseconds) t0))
171 (when out
172 (print-report m out digits elapsed inputs)))
173
174 (define (run-interactive m)
175 (displayln "")
176 (displayln (format "-- ~a --" (method-id m)))
177 (define inputs (prompt-inputs m))
178 (unless (eof-object? inputs)
179 (define digits (display-digits-for m inputs))
180 (unless (eof-object? digits)
181 (execute m inputs digits))))
182
183 (define (compare-all)
184 (displayln "")
185 (displayln " Every method at its default input:")
186 (displayln "")
187 (printf " ~a ~a ~a ~a ~a~n"
188 (~a "method" #:min-width 16)
189 (~a "input" #:min-width 18)
190 (~a "time" #:min-width 12)
191 (~a "correct digits" #:min-width 15)
192 "value")
193 (for ([m (in-list methods)])
194 (define inputs (method-compare-inputs m))
195 (define digits (min 20 (max 1 ((method-auto-digits m) inputs))))
196 (define t0 (current-inexact-milliseconds))
197 (define out
198 (with-handlers ([exn:fail? (lambda (e) #f)])
199 (run-method m inputs digits)))
200 (define elapsed (- (current-inexact-milliseconds) t0))
201 (define input-desc
202 (string-join (for/list ([k (in-list (hash-keys inputs))])
203 (format "~a=~a" k (hash-ref inputs k)))
204 ","))
205 (cond
206 [(not out)
207 (printf " ~a ~a ~a ~a ~a~n"
208 (~a (method-id m) #:min-width 16)
209 (~a input-desc #:min-width 18)
210 (~a (format-elapsed elapsed) #:min-width 12)
211 (~a "FAILED" #:min-width 15)
212 "-")]
213 [(outcome-hex out)
214 (printf " ~a ~a ~a ~a ~a~n"
215 (~a (method-id m) #:min-width 16)
216 (~a input-desc #:min-width 18)
217 (~a (format-elapsed elapsed) #:min-width 12)
218 (~a (format "~a hex" (count-hex-matches out)) #:min-width 15)
219 (format "3.~a" (outcome-hex out)))]
220 [else
221 (define v (outcome-value out))
222 (define ref (reference-value (method-constant m) (+ digits 10)))
223 (define match (count-matching-digits (fp->digits v digits) (fp->digits ref digits)))
224 (printf " ~a ~a ~a ~a ~a~n"
225 (~a (method-id m) #:min-width 16)
226 (~a input-desc #:min-width 18)
227 (~a (format-elapsed elapsed) #:min-width 12)
228 (~a (max 0 (sub1 match)) #:min-width 15)
229 (fp-display-str v (min digits 18)))])))
230
231 (define (count-hex-matches out)
232 (define hx (outcome-hex out))
233 (define ref (reference-value 'pi (exact-ceiling (+ 10 (* 1.21 (string-length hx))))))
234 (count-matching-digits (string-downcase hx)
235 (fp->hex-frac-digits ref (string-length hx))))
236
237 (define (start-repl)
238 (print-banner)
239 (print-method-table)
240 (displayln "")
241 (displayln " Enter a number or id to run a method. Other commands:")
242 (displayln " c compare every method with its defaults")
243 (displayln " h <id> help for one method")
244 (displayln " s <n> digits of output to display (now 200)")
245 (displayln " l list methods again")
246 (displayln " q quit")
247 (let loop ()
248 (define cmd (read-choice "\nchoice> "))
249 (cond
250 [(eof-object? cmd) (displayln "\nbye")]
251 [(string=? cmd "") (loop)]
252 [(member cmd '("q" "quit" "exit" "bye")) (displayln "bye")]
253 [(member cmd '("l" "list" "?")) (print-method-table) (loop)]
254 [(member cmd '("c" "compare")) (compare-all) (loop)]
255 [(string-prefix? cmd "h")
256 (handle-help (string-trim (substring cmd 1)))
257 (loop)]
258 [(string-prefix? cmd "s")
259 (handle-show (string-trim (substring cmd 1)))
260 (loop)]
261 [else
262 (define m (method-lookup cmd))
263 (if m
264 (begin (run-interactive m) (loop))
265 (begin (printf " unknown choice: ~a (try 'l')~n" cmd) (loop)))])))
266
267 (define (handle-help arg)
268 (cond [(string=? arg "")
269 (displayln " usage: h <method-id or number>")
270 (print-method-table)]
271 [else
272 (define m (method-lookup arg))
273 (if m (print-help m) (printf " unknown method: ~a~n" arg))]))
274
275 (define (handle-show arg)
276 (define n (string->number arg))
277 (cond [(not (and (number? n) (exact-integer? n) (positive? n)))
278 (displayln (format " digits to display is currently ~a; usage: s <n>" (show-limit)))]
279 [(> n 1000000)
280 (displayln " that is a lot; capped at 1000000")
281 (show-limit 1000000)]
282 [else
283 (show-limit n)
284 (displayln (format " will display up to ~a digits" n))]))
285
286 ;; =====================================================================
287 ;; Section 8: command line entry point
288 ;; =====================================================================
289
290 (module+ main
291 (define mode 'interactive)
292 (define id-arg #f)
293 (define sets '())
294 (define show 200)
295
296 (command-line
297 #:program "math.rkt"
298 #:once-each
299 [("--list") "Print the method table and exit"
300 (set! mode 'list)]
301 [("--compare") "Run every method with its default inputs"
302 (set! mode 'compare)]
303 [("--show") n "Digits of the result to display (default 200)"
304 (set! show (string->number n))]
305 #:multi
306 [("--set") key value "Set an input, e.g. --set digits 1000"
307 (set! sets (cons (cons (string->symbol key) value) sets))]
308 #:args arg-list
309 ;; Flags must precede the method id; after it, inputs are key=value.
310 (when (pair? arg-list)
311 (set! mode 'run)
312 (set! id-arg (car arg-list))
313 (for ([extra (in-list (cdr arg-list))])
314 (match (string-split extra "=")
315 [(list "show" v) (set! show (string->number v))]
316 [(list key v) (set! sets (cons (cons (string->symbol key) v) sets))]
317 [_ (raise-user-error 'math.rkt "expected key=value, got: ~a" extra)]))))
318
319 (show-limit (if (and (number? show) (positive? show)) show 200))
320
321 ;; Inputs come from --set, falling back to each method's default.
322 (define (inputs-from-sets m)
323 (for/fold ([h (hash)]) ([spec (in-list (method-inputs m))])
324 (define key (input-spec-key spec))
325 (define raw (assv key sets))
326 (define n (and raw (string->number (cdr raw))))
327 (hash-set h key
328 (if (and (exact-integer? n) (positive? n))
329 (min n (input-spec-cap spec))
330 (input-spec-default spec)))))
331
332 (case mode
333 [(list)
334 (print-method-table)]
335 [(compare)
336 (compare-all)]
337 [(run)
338 (define m (and (string? id-arg) (method-lookup id-arg)))
339 (cond
340 [(not m)
341 (printf "unknown method: ~a~n" id-arg)
342 (print-method-table)
343 (exit 1)]
344 [else
345 (define inputs (inputs-from-sets m))
346 (define digits (max 1 ((method-auto-digits m) inputs)))
347 (execute m inputs digits)])]
348 [else
349 (start-repl)]))
Two Racket idioms in this section are worth adding to your toolbox. The cond
clause with => in verification-lines binds the non-false test value to a
name, avoiding a nested if. And module+ main is Racket’s way of separating
“running the file” from “requiring the file”: the test suite can require
math.rkt and get all the exports without triggering the command line
interface, because the main submodule only runs when the file is executed
directly.
The Test Suite
Tests for numerical code need trusted data, and the most trustworthy data
available for this problem is the published digit expansion of the constants
themselves. The test file opens by embedding 100 digits of , 100 digits
of , and 32 hexadecimal digits of as string constants. Every
subsequent check compares computed digits against these tables, or against
values produced by an independent algorithm that has itself been pinned against
the tables.
1 ;; Published expansions of pi and e (fractional digits).
2 (define PI-50 "14159265358979323846264338327950288419716939937510")
3 (define PI-99
4 (string-append PI-50 "5820974944592307816406286208998628034825342117067"))
5 (define E-50 "71828182845904523536028747135266249775724709369995")
6 (define E-99
7 (string-append E-50 "9574966967627724076630353547594571382178525166427"))
8
9 ;; Hex expansion of pi after the point.
10 (define PI-HEX-32 "243F6A8885A308D313198A2E03707344")
One subtlety shapes all the string comparisons. Integer division truncates, so a
computed value at scale can sit one unit below the true value in the final
digit. The helper frac asks for five more digits than it returns and
discards them, which guarantees the last kept digit is a true digit of the
constant and not an artifact of truncation. With that helper in place, the test
strategy has three layers:
- Pin the fast methods (Chudnovsky, Machin, AGM, and the exponential series) against the published tables, then cross-check them against each other at 1000 digits, far beyond anything floating point could reach.
- Verify the slow historical methods against the fast references, at the accuracy their theory predicts, with a digit of slack.
- Check the structural properties each algorithm guarantees: Archimedes’ bounds bracket at every iteration, Viete’s product increases monotonically,
consecutive exact partial sums of differ by exactly , and the
continued fraction convergents match hand computed fractions like
, 
, and 
.
Here is the complete test file:
1 #lang racket
2
3 ;; Tests for math.rkt. Run with: racket math-tests.rkt (or: raco test .)
4 ;;
5 ;; Strategy:
6 ;; 1. pin the fast methods against the published digits of pi and e,
7 ;; 2. cross-check the slower historical methods against those,
8 ;; 3. check the structural properties each algorithm guarantees
9 ;; (Archimedes brackets pi, Viete climbs monotonically, the exact
10 ;; rational partial sums of e match their hand-computed fractions).
11 ;;
12 ;; String comparisons are made on digits that rounding cannot touch: the
13 ;; helper `frac` asks for one more digit than it returns, so the final
14 ;; digit of every expectation is a true digit of the constant.
15
16 (require rackunit
17 rackunit/text-ui
18 "math.rkt")
19
20 ;; Published expansions of pi and e (fractional digits).
21 (define PI-50 "14159265358979323846264338327950288419716939937510")
22 (define PI-99
23 (string-append PI-50 "5820974944592307816406286208998628034825342117067"))
24 (define E-50 "71828182845904523536028747135266249775724709369995")
25 (define E-99
26 (string-append E-50 "9574966967627724076630353547594571382178525166427"))
27
28 ;; Hex expansion of pi after the point.
29 (define PI-HEX-32 "243F6A8885A308D313198A2E03707344")
30
31 (define S 130) ; working scale: comfortably more precision than any check needs
32
33 ;; The first `n` fractional digits of a fixed-point value. A few extra
34 ;; digits are computed and discarded so that a value which sits one unit
35 ;; below the truth at its working scale cannot shift the last digit shown.
36 (define (frac x n)
37 (substring (fp->digits x (+ n 5)) 1 (add1 n)))
38
39 (define (correct-digits x ref n)
40 (max 0 (sub1 (count-matching-digits (fp->digits x n) (fp->digits ref n)))))
41
42 (define PI-REF (fp (pi-chudnovsky S) S))
43 (define E-REF (fp (e-series S) S))
44
45 (define decimal-tests
46 (test-suite
47 "fixed-point formatting"
48 (check-equal? (fp-display-str (fp 314159 5) 4) "3.1415") ; truncated, not rounded
49 (check-equal? (fp-display-str (fp 314159 5) 5) "3.14159")
50 (check-equal? (fp-display-str (fp 314159 5) 7) "3.1415900")
51 (check-equal? (fp-display-str (fp 5 1) 0) "0") ; 0.5 with no digits left
52 (check-equal? (fp-display-str (fp 35 1) 0) "3")
53 (check-equal? (fp-display-str (fp 1/2 0) 3) "0.500") ; exact rational at scale 0
54 (check-equal? (fp-display-str (fp 8/3 0) 6) "2.666666")
55 ;; the denominator must be taken after rescaling, not before
56 (check-equal? (fp-display-str (fp (e-series-exact 40) 0) 20)
57 "2.71828182845904523536")
58 (check-equal? (count-matching-digits "314159" "314159") 6)
59 (check-equal? (count-matching-digits "314159" "314259") 3)
60 (check-equal? (count-matching-digits "31" "314159") 2)
61 (check-equal? (fp->hex-frac-digits (fp 5 1) 4) "8000") ; 0.5 = 0.8 in hex
62 (check-equal? (fp->hex-frac-digits (fp 3141592653589793238 18) 4) "243f")))
63
64 (define fast-pi-tests
65 (test-suite
66 "fast pi algorithms"
67 (check-equal? (frac PI-REF 99) PI-99)
68 (check-equal? (frac (fp (pi-machin S) S) 99) PI-99)
69 (check-equal? (frac (fp (pi-agm S) S) 99) PI-99)
70 (check-equal? (frac (reference-value 'pi 99) 99) PI-99)
71 ;; far beyond flonum range the three independent algorithms must agree
72 (check-equal? (fp->digits (fp (pi-chudnovsky 1025) 1025) 1000)
73 (fp->digits (fp (pi-machin 1025) 1025) 1000))
74 (check-equal? (fp->digits (fp (pi-chudnovsky 1025) 1025) 1000)
75 (fp->digits (fp (pi-agm 1025) 1025) 1000))))
76
77 (define archimedes-tests
78 (test-suite
79 "Archimedes polygon iteration"
80 (let-values ([(lo hi) (pi-archimedes 0 S)])
81 (check-equal? (fp-display-str (fp lo S) 6) "3.000000")
82 (check-equal? (fp-display-str (fp hi S) 6) "3.464101")) ; hexagon: 3 and 2*sqrt(3)
83 (for ([iters (in-list '(1 5 20 60))])
84 (define-values (lo hi) (pi-archimedes iters S))
85 (check-true (< lo (fp-v PI-REF) hi)
86 (format "~a doublings must bracket pi" iters)))
87 ;; one doubling cuts the bracket by roughly a factor of four
88 (let-values ([(lo1 hi1) (pi-archimedes 20 S)]
89 [(lo2 hi2) (pi-archimedes 21 S)])
90 (check-true (< (* 3 (- hi2 lo2)) (- hi1 lo1))))
91 ;; and the bracket keeps tightening all the way down
92 (check-true (>= (correct-digits (let-values ([(lo _h) (pi-archimedes 200 S)])
93 (fp lo S))
94 PI-REF 130)
95 119))))
96
97 (define viete-tests
98 (test-suite
99 "Viete's nested radicals"
100 (check-true (< (pi-viete 1 S) (fp-v PI-REF)))
101 (for/fold ([prev (pi-viete 1 S)]) ([n (in-range 2 30)])
102 (define cur (pi-viete n S))
103 (check-true (> cur prev) "Viete's product increases toward pi")
104 (check-true (< cur (fp-v PI-REF)) "Viete's product stays below pi")
105 cur)
106 (check-equal? (frac (fp (pi-viete 250 S) S) 99) PI-99)))
107
108 (define slow-pi-tests
109 (test-suite
110 "historical pi series"
111 (check-equal? (fp-display-str (fp (pi-leibniz 1 S) S) 3) "4.000")
112 (check-equal? (fp-display-str (fp (pi-leibniz 2 S) S) 3) "2.666") ; 8/3
113 (check-equal? (fp-display-str (fp (pi-nilakantha 1 S) S) 3) "3.166") ; 3 + 1/6
114 ;; an alternating series errs by less than its first omitted term: 1/(2n)
115 (for ([n (in-list '(10 1000 100000))])
116 (define err (abs (- (pi-leibniz n S) (fp-v PI-REF))))
117 (check-true (< err (quotient (* 4 (expt 10 S)) (* 2 n)))
118 (format "Leibniz error at ~a terms exceeds 1/(2n)" n)))
119 ;; measured accuracy at these term counts, with a digit of slack
120 (check-true (>= (correct-digits (fp (pi-leibniz 200000 S) S) PI-REF 20) 4))
121 (check-true (>= (correct-digits (fp (pi-nilakantha 20000 S) S) PI-REF 20) 12))
122 (check-true (>= (correct-digits (fp (pi-wallis 100000 S) S) PI-REF 20) 4))
123 ;; Nilakantha beats Leibniz badly at equal term counts
124 (check-true (> (correct-digits (fp (pi-nilakantha 1000 S) S) PI-REF 30)
125 (correct-digits (fp (pi-leibniz 1000 S) S) PI-REF 30)))))
126
127 (define bbp-tests
128 (test-suite
129 "BBP spigot"
130 (check-equal? (pi-bbp-hex 1 32) PI-HEX-32)
131 ;; deep digits must agree with the decimal expansion converted to hex
132 (let ([ref (reference-value 'pi 400)])
133 (check-equal? (string-downcase (pi-bbp-hex 100 16))
134 (substring (fp->hex-frac-digits ref 300) 99 115))
135 (check-equal? (string-downcase (pi-bbp-hex 250 8))
136 (substring (fp->hex-frac-digits ref 300) 249 257)))))
137
138 (define e-tests
139 (test-suite
140 "e algorithms"
141 (check-equal? (frac E-REF 99) E-99)
142 (check-equal? (frac (reference-value 'e 99) 99) E-99)
143 ;; 500 terms of the series settle 1000 digits
144 (check-equal? (fp->digits (fp (e-series 1025) 1025) 1000)
145 (fp->digits (fp (e-series-terms 500 1025) 1025) 1000))
146 ;; exact rational partial sums, hand-checked
147 (check-equal? (e-series-exact 1) 1)
148 (check-equal? (e-series-exact 2) 2)
149 (check-equal? (e-series-exact 3) 5/2)
150 (check-equal? (e-series-exact 4) 8/3)
151 (check-equal? (e-series-exact 5) 65/24)
152 (check-equal? (e-series-exact 6) 163/60)
153 ;; consecutive partial sums differ by exactly the next term, 1/k!
154 (check-equal? (- (e-series-exact 41) (e-series-exact 40))
155 (/ 1 (for/product ([i (in-range 1 41)]) i)))
156 (check-true (< (e-series-exact 40) (e-series-exact 60)))
157 ;; the decimal view of the exact sum agrees with the fixed-point one
158 (check-equal? (fp->digits (fp (e-series-exact 60) 0) 50)
159 (fp->digits (fp (e-series-terms 60 80) 80) 50))
160 ;; (1 + 1/n)^n climbs toward e from below and never passes it
161 (for/fold ([prev 0]) ([n (in-list '(1 10 1000 100000))])
162 (define cur (e-limit n S))
163 (check-true (> cur prev) "the limit definition increases with n")
164 (check-true (< cur (fp-v E-REF)) "the limit definition stays below e")
165 cur)
166 (check-true (>= (correct-digits (fp (e-limit 1000000 S) S) E-REF 20) 5))
167 (check-true (> (correct-digits (fp (e-limit 10000000 S) S) E-REF 20)
168 (correct-digits (fp (e-limit 1000 S) S) E-REF 20)))
169 ;; continued fraction convergents, hand-checked
170 (check-equal? (e-continued-fraction 1) 2)
171 (check-equal? (e-continued-fraction 2) 3)
172 (check-equal? (e-continued-fraction 3) 8/3)
173 (check-equal? (e-continued-fraction 4) 11/4)
174 (check-equal? (e-continued-fraction 5) 19/7)
175 (check-equal? (e-continued-fraction 6) 87/32)
176 (check-equal? (e-continued-fraction 7) 106/39)
177 (check-equal? (e-continued-fraction 8) 193/71)
178 (check-equal? (fp->digits (fp (e-continued-fraction 80) 0) 50)
179 (fp->digits (reference-value 'e 50) 50))))
180
181 (define method-table-tests
182 (test-suite
183 "method table"
184 (check-equal? (length methods) 14)
185 (check-equal? (length (remove-duplicates (map method-id methods))) 14)
186 (check-equal? (length (filter (lambda (m) (eq? 'pi (method-constant m))) methods)) 9)
187 (check-equal? (length (filter (lambda (m) (eq? 'e (method-constant m))) methods)) 5)
188 ;; lookup by index and by id agree
189 (for ([i (in-range 1 (add1 (length methods)))])
190 (define m (method-lookup (number->string i)))
191 (check-not-false m)
192 (check-eq? m (method-lookup (symbol->string (method-id m)))))
193 (check-false (method-lookup "no-such-method"))
194 (check-false (method-lookup "99"))
195 ;; every method runs at its default inputs and reports something usable
196 (for ([m (in-list methods)])
197 (define inputs (method-compare-inputs m))
198 (for ([spec (in-list (method-inputs m))])
199 (check-true (hash-has-key? inputs (input-spec-key spec))
200 (format "~a compare-inputs is missing ~a"
201 (method-id m) (input-spec-key spec)))
202 (check-true (<= (hash-ref inputs (input-spec-key spec)) (input-spec-cap spec))
203 (format "~a compare-inputs exceeds the cap on ~a"
204 (method-id m) (input-spec-key spec))))
205 (define digits (max 1 ((method-auto-digits m) inputs)))
206 (define out (run-method m inputs digits))
207 (check-pred outcome? out (format "~a produced no outcome" (method-id m)))
208 (check-pred pair? (outcome-lines out) (format "~a printed nothing" (method-id m)))
209 (check-not-false (or (outcome-value out) (outcome-hex out))
210 (format "~a produced no value" (method-id m))))))
211
212 (define outcome-format-tests
213 (test-suite
214 "outcome formatting"
215 (let*-values ([(one-seventh)
216 (values (fp (quotient (expt 10 65) 7) 65))] ; 1/7 = 0.142857...
217 [(int-str frac-str digits-str)
218 (fp->decimal one-seventh 60)]
219 [(bbp-out)
220 (values (run-method (method-lookup "bbp") (hash 'start 1 'count 8) 1))]
221 [(arch-lines)
222 (values (outcome-lines
223 (run-method (method-lookup "archimedes")
224 (hash 'iterations 40) 24)))])
225 (check-equal? int-str "0")
226 (check-equal? frac-str (string-append "1428571428" "5714285714" "2857142857"
227 "1428571428" "5714285714" "2857142857"))
228 (check-equal? digits-str (string-append "0" frac-str))
229 ;; 60 digits = 6 groups of ten, five groups per line
230 (check-equal? (format-number-lines int-str frac-str 200)
231 (list " 0.1428571428 5714285714 2857142857 1428571428 5714285714"
232 " 2857142857"))
233 ;; a short display limit truncates and says so, showing the tail
234 (check-equal? (format-number-lines int-str frac-str 25)
235 (list " 0.1428571428 5714285714 28571"
236 " ... showing 25 of 60 digits; last 10: 2857142857"))
237 ;; zero digits prints no decimal point
238 (check-equal? (format-number-lines "3" "" 200) (list " 3"))
239 ;; the BBP outcome carries hex digits instead of a decimal value
240 (check-false (outcome-value bbp-out))
241 (check-equal? (outcome-hex bbp-out) (substring PI-HEX-32 0 8))
242 ;; Archimedes reports both bounds plus how many digits they pin down
243 (check-true (string-contains? (car arch-lines) "inscribed 6597069766656-gon"))
244 (check-true (ormap (lambda (l) (string-contains? l "circumscribed")) arch-lines))
245 (check-true (ormap (lambda (l) (string-contains? l "digits pinned by the bracket: 24"))
246 arch-lines)))))
247
248 (run-tests
249 (test-suite
250 "math.rkt"
251 decimal-tests
252 fast-pi-tests
253 archimedes-tests
254 viete-tests
255 slow-pi-tests
256 bbp-tests
257 e-tests
258 method-table-tests
259 outcome-format-tests))
The method-table-tests suite deserves a special mention because it tests the
framework rather than the mathematics: it asserts that the registry holds 14
methods with unique ids, that lookup by number and by name agree, and that every
single method runs to a usable outcome at its default inputs. If you add a
fifteenth method and forget a compare input or exceed a cap, this suite fails
immediately and tells you which method is misconfigured.
Running the Code
The program needs only Racket itself, with no external packages. Running the
file with no arguments starts the interactive menu; the --list, --compare,
and --show flags and the key=value input syntax cover the non-interactive
uses. The following outputs were captured from real runs of the program.
Listing the method table:
1 $ racket math.rkt --list
2 # id method input
3 -- ---------------- -------------------------------------- -------------
4 1 leibniz Gregory-Leibniz series terms
5 2 nilakantha Nilakantha series terms
6 3 wallis Wallis product terms
7 4 viete Viete's nested radicals terms
8 5 archimedes Archimedes polygon doubling iterations
9 6 machin Machin-like arctangent formula digits
10 7 agm Gauss-Legendre / Brent-Salamin AGM digits
11 8 chudnovsky Chudnovsky binary splitting digits
12 9 bbp BBP hex-digit spigot start,count
13
14 10 e-series Exponential series to N digits digits
15 11 e-series-terms Exponential series, first n terms terms
16 12 e-series-exact Exponential series as an exact fraction terms
17 13 e-limit Limit definition (1 + 1/n)^n n
18 14 e-cf Continued fraction convergents terms
Running one method, Machin’s formula, at 20 digits:
1 $ racket math.rkt machin digits=20
2 Machin-like arctangent formula
3 pi = 16*arctan(1/5) - 4*arctan(1/239)
4 about 1.4 digits per arctan term; the classic hand-computation formula
5 inputs: digits=20
6 3.1415926535 8979323846
7 check: all 20 displayed digits match the reference (chudnovsky at 30 digits)
8 time: 20 us
Archimedes’ method shows both sides of its bracket, along with the polygon size and the number of digits the bracket certifies:
1 $ racket math.rkt --show 60 archimedes iterations=25
2 Archimedes polygon doubling
3 a <- 2ab/(a+b); b <- sqrt(a*b), from the hexagon a=2*sqrt(3), b=3
4 linear: about 0.6 digits per doubling, but every step brackets pi exactly
5 inputs: iterations=25
6 inscribed 201326592-gon (lower bound):
7 3.1415926535 89793
8 circumscribed 201326592-gon (upper bound): 3.141592653589793
9 digits pinned by the bracket: 15
10 check: all 15 displayed digits match the reference (chudnovsky at 25 digits)
11 time: 49 us
The BBP spigot pulling 16 hexadecimal digits out of nowhere:
1 $ racket math.rkt bbp start=1 count=16
2 BBP hex-digit spigot
3 pi = sum 16^-k (4/(8k+1) - 2/(8k+4) - 1/(8k+5) - 1/(8k+6))
4 a spigot: digit d in base 16 with none of the earlier digits computed
5 inputs: start=1, count=16
6 hex digits 1..16 of pi, after the decimal point:
7 3.243F6A8885A308D3
8 check: all 16 hex digits match the reference pi
9 time: 59 us
The Chudnovsky method at 500 digits, with the display limited to the first 100:
1 $ racket math.rkt --show 100 chudnovsky digits=500
2 Chudnovsky binary splitting
3 1/pi = 12 * sum (-1)^k (6k)!(13591409+545140134k) / ((3k)!(k!)^3 640320^(3k+3/2))
4 about 14.18 digits per term; used for every modern pi record
5 inputs: digits=500
6 3.1415926535 8979323846 2643383279 5028841971 6939937510
7 5820974944 5923078164 0628620899 8628034825 3421170679
8 ... showing 100 of 500 digits; last 10: 8301194912
9 check: all 500 displayed digits match the reference (chudnovsky at 510 digits)
10 time: 76 us
The comparison report runs all fourteen methods at their default inputs and tabulates time, correct digits, and value:
1 $ racket math.rkt --compare
2
3 Every method at its default input:
4
5 method input time correct digits value
6 leibniz terms=200000 7.23 ms 4 3.1415
7 nilakantha terms=5000 458 us 11 3.14159265358
8 wallis terms=50000 13 ms 3 3.141
9 viete terms=40 46 us 20 3.141592653589793238
10 archimedes iterations=30 40 us 18 3.141592653589793238
11 machin digits=50 6 us 20 3.141592653589793238
12 agm digits=50 12 us 20 3.141592653589793238
13 chudnovsky digits=50 5 us 20 3.141592653589793238
14 bbp start=1,count=24 116 us 24 hex 3.243F6A8885A308D313198A2E
15 e-series digits=50 6 us 20 2.718281828459045235
16 e-series-terms terms=40 6 us 20 2.718281828459045235
17 e-series-exact terms=25 9 us 20 2.718281828459045235
18 e-limit n=100000 8 us 4 2.7182
19 e-cf terms=40 11 us 20 2.718281828459045235
Finally, the test suite:
1 $ racket math-tests.rkt
2 246 success(es) 0 failure(s) 0 error(s) 246 test(s) run
Interpreting the Results
The comparison table is the payoff of the entire chapter, because it makes the theoretical convergence rates from the opening section visible as measured facts.
Look at the two extremes. Leibniz burns 200,000 terms and 7.23 milliseconds to
produce 4 correct digits. Chudnovsky computes 50 correct digits in 5
microseconds, using about 5 terms of its series. That is roughly a thousandfold
difference in work for more than ten times the accuracy, and it is exactly what
the error bounds predicted: error for Leibniz versus 14.18 digits per
term for Chudnovsky. The Nilakantha row sits between them, 11 digits from 5,000
terms, matching its cubic convergence. When theory and measurement line up this
cleanly, you know both the formulas and the implementations are right.
The e-limit row teaches the complementary lesson. The limit definition
is how is usually introduced, yet 100,000 exponentiations
worth of input yields only 4 digits. The continued fraction, which almost never
appears in a first calculus course, produces 20 digits from 40 partial
quotients. The choice of representation matters more than the choice of
programming language or the speed of the machine.
The Archimedes run deserves a second look. After 25 doublings the program is
tracking polygons with 201,326,592 sides, and it reports “digits pinned by the
bracket: 15.” That number is different in kind from every other digit count in
the chapter. It is not a statistical statement about agreement with another
algorithm; it is a proof. The true value of is guaranteed to lie between
the inscribed and circumscribed values, so those 15 digits are certain.
Archimedes had a theorem, not just an approximation, and the program preserves
that property.
The verification lines answer the question posed at the start: how do you know
the digits are right? Every run is checked against the Chudnovsky series
computed at higher precision, and every run against the exponential
series. The test suite goes further and makes the fast algorithms check each
other at 1000 digits, where the probability that two independent buggy
implementations agree by chance is effectively zero. The BBP check is the most
satisfying: the spigot works in floating point in base 16, the reference works
in exact integers in base 10, and the hex digits match anyway.
Timing is the one place to be careful with interpretation. The microsecond timings for Machin, AGM, and Chudnovsky at 50 digits reflect that Racket’s arithmetic on small bignums is extremely fast; at 100,000 digits the picture separates, with binary splitting pulling ahead of the AGM’s repeated full precision square roots. The compare table is a convergence demonstration, not a benchmark.
Wrap Up
This chapter used two of the oldest computational problems in mathematics as a
vehicle for three transferable ideas. The fixed-point representation
replaces floating point with exact integer arithmetic whenever
results are judged digit by digit, and the same pattern serves in financial
calculation and anywhere else rounding is the enemy. The method registry shows
how a table of structures with function fields can replace a pile of special
case user interface code; the menu, the command line, the comparison report,
and the test suite all iterate the same fourteen entries. And the layered test
strategy, pin against published data, cross-check independent algorithms, then
verify structural guarantees, is a template for testing any numerical code
where a single “expected value” cannot be typed into the test file by hand.
Along the way the chapter walked through twenty five centuries of mathematical history, from Archimedes’ bracketed polygons to a spigot formula that produces isolated hexadecimal digits, and saw each era’s algorithm run, verified, in a few microseconds.
Optional Practice Problems
The following exercises extend the example code. They are ordered roughly from gentle to ambitious.
Euler’s arctangent identity. Machin’s formula is not the only identity of its kind. Euler showed that
.
Implement (pi-euler-arctan s)using the existingarctan-reciphelper, register it as a fifteenth method in themethodstable, and confirm thatracket math-tests.rktstill passes (remember thatmethod-table-testscounts the methods). How many correct digits does it produce at the same scale wherepi-machinproduces 50? Explain the difference in terms of the convergence of
for larger 
.Measure the convergence rates yourself. Write a small script that requires
math.rktand, for the Leibniz, Nilakantha, and Wallis methods, computes the number of correct digits (usingcorrect-digits-style comparison againstreference-value) at inputs of 100, 1000, 10000, and 100000. Tabulate digits gained per tenfold increase and compare against the theoretical rates of 1 and 3 digits per tenfold stated in the registry. Why does the measured Nilakantha rate sometimes look better than the theory?Ramanujan’s series. In 1914 Ramanujan published the astonishing formula
Implement it in the same fixed-point style as pi-chudnovsky (you do not
need binary splitting; a direct term recurrence is fine at moderate scales),
register it as a method, and verify it against the reference at 200 digits.
About how many digits does each term buy, and how does that compare with
Chudnovsky’s 14.18?
A spigot for log 2. The BBP idea works for any constant with a suitable series. The natural logarithm of 2 satisfies
, a BBP-type series in base 2.
Adapt modpow16andbbp-seriesto extract binary (or hexadecimal) digits of
at a given position, and test them against exact rational
partial sums of the series converted to hex with fp->hex-frac-digits.Continued fractions are best approximations. A convergent
of the
continued fraction for satisfies
. Add a test suite that checks this inequality exactly
for the first 40 convergents: compute the rational error using
e-series-exactat high term count as the reference value of, and
verify 
with pure rational
arithmetic. Then check the stronger property that each convergent is closer
to than any fraction with a smaller denominator, for denominators up to
1000.