After reading Carl Sagan’s Contact I’m crushing on π, so I looked it up on Wikipedia, where, of course, several algorithms for computing it are described in detail. One of them, the Brent-Salamin algorithm, has the really neat property of doubling the number of digits it gets right on each iteration.
The algorithm starts with four terms (images ripped from Wikipedia):
And then iterates on them like:
Whenever you get tired of iterating, you can just plug the latest versions of the terms into:
In Clojure, we’ll start by defining a function to do the heavy lifting:
(defn brent-salamin-term
"compute terms of the the Brent-Salamin algorithm"
[])
We’ll want two versions of this function, one that will be called with no arguments, to get the initial values of the parameters:
([] [1 (/ 1 (sqrt 2)) (/ 1 4) 1])
And a second which will implement a single iteration of the algorithm:
([a b t p]
(let [A (/ (+ a b) 2)
B (sqrt (* a b))
T (- t (* p (sqr (- a A))))
P (* 2 p)]
[A B T P])))
Finally, we’ll define a short function to calculate π given Brent-Salamin terms:
(defn brent-salamin-pi
[a b t _]
(/ (sqr (+ a b)) (* 4 t)))
For idiomatic access to the values generated by this algorithm, we’ll define a lazy sequence:
(def brent-salamin-terms
(iterate #(apply brent-salamin-term %) (brent-salamin-term)))
And another for the values of π:
(def brent-salamin-pis (map #(apply brent-salamin-pi %) brent-salamin-terms))
Things look pretty good:
user> (take 5 brent-salamin-pis)
(2.914213562373095 3.1405792505221686 3.141592646213543
3.141592653589794 3.141592653589794)
But wait - we should probably be careful with precision. Ideally, I’d like to be able to get an arbitrarily precise version of π. Let’s take a look at what kind of numeric types we’re dealing with by printing the class of the inputs to our terms algorithm:
([a b t p]
(prn (class a))
(prn (class b))
(prn (class t))
(prn (class p))
(let [A (/ (+ a b) 2)
B (sqrt (* a b))
T (- t (* p (sqr (- a A))))
P (* 2 p)]
[A B T P]))
user> (apply brent-salamin-term
(apply brent-salamin-term (brent-salamin-term)))
java.lang.Integer
java.lang.Double
clojure.lang.Ratio
java.lang.Integer
java.lang.Double
java.lang.Double
java.lang.Double
java.lang.Integer
[0.8472249029234942 0.8472012667468914 0.22847329108090064 4]
Ugh - what a mess. We start off with integers in the first iteration, and move quickly to doubles. That’s fine for the real world, where according to Wikipedia “a physicist needs only 39 digits of Pi to make a circle the size of the observable universe accurate to one atom of hydrogen,” but we’re trying for gold here. Let’s start by using BigDecimals wherever we can (see sidebar).
This is better thanks to Clojure’s smart arithmetic operators - if they receive a BigDecimal they’ll do BigDecimal operations - but there’s still a problem in our sqrt function. Unfortunately, Java doesn’t give us a BigDecimal square root function out of the box, so we’ll have to roll our own. Michael Gilleland at Merriam Park Software explains Heron’s method in Big Square Roots, and while the Java implementation leaves something to be desired, a basic Clojure implementation is actually quite nice:
(defn sqrt
[n]
(loop [n n g 1]
(let [next-guess (/ (+ (/ n g) g) 2)]
(if (= g next-guess) g (recur n next-guess)))))
That is, until we compile the file, at which point everything explodes!
error: java.lang.ArithmeticException: Non-terminating decimal expansion;
no exact representable decimal result.
The problem here is that we haven’t specified a precision for the BigInteger operations that happen when we create our lazy sequences at compile time, and indeed, we don’t want to. Instead, we’ll make them into functions so that callers can choose whatever precision they want when they are created. This gets everything compiling, and we jump over to the REPL to give it a whirl:
user> (with-precision 10 (take 10 (etc/brent-salamin-pis)))
java.lang.ArithmeticException: Non-terminating decimal expansion;
no exact representable decimal result.
The problem (and it’s the last one, I promise!) is that with-precision creates a thread local binding which isn’t available when the lazy sequence is realized later on. Fortunately, Clojure 1.1 includes a handy function for just this purpose, bound-fn. bound-fn creates a function that carries thread local bindings along into whatever context it is evaluated.
(defn brent-salamin-terms
[]
(iterate (bound-fn [t] (apply brent-salamin-term t)) (brent-salamin-term)))
(defn brent-salamin-pis
[]
(map (bound-fn [t] (apply brent-salamin-pi t)) (brent-salamin-terms)))
Finally - success!
user> (with-precision 10 (take 5 (etc/brent-salamin-pis)))
(2.914213562M 3.140579251M 3.141592648M 3.141592656M 3.141592656M)
This turned into a bit of a slog at the end, but only when we started caring quite a bit about precision. The default behavior of the language leads to very elegant code organized around lazy sequences that will work in almost all real world applications, and that seems like a major win.