Embedded Files
Posted by Uncle Bob on 08/05/2009
The SICP book shows a cute algorithm for computing PI using infinite sequences (lazy evaluation). I have adapted it for Clojure and modified it to use Java’s BigDecimal class, so that I can generate arbitrary precision.
First we define a function that generates an “infinite” list of terms of the form 1/1 + 1/3 – 1/5 + 1/7 – ... This is a well known series that very slowly converges on Pi. Notice that I’m using the ‘M’ suffix to denote BigDecimal, and I’m restricting the precision to 300 places. Notice also that I’m using lazy-seq. This function delays the execution of the ‘cons’ term until that part of the list is actually used. This allows the list to appear to be infinite.
(defn pi-summands [n] (lazy-seq (cons (with-precision 300 (/ 1M n)) (map - (pi-summands (+ n 2))))))
Next we have two utility functions. The first multiplies the terms of a list by a scale factor. The other adds two lists together producing a list of sums. If the input lists are ‘infinite’, the output lists will also be ‘infinite’.
(defn scale-stream [stream factor] (map (fn [x] (* x factor)) stream)) (defn add-streams [p q] (map + p q))
partial-sums takes a list of terms produces a list of the sums of those terms. Given [a b c] it produces [a a+b a+b+c]. And, again, the lists can be infinite.
(defn partial-sums [s] (lazy-seq (cons (first s) (add-streams (next s) (partial-sums s)))))
pi-stream mutiplies the partial-sums of the pi-summands by 4. Turning the summing sequence into 4/1 + 4/3 – 4/5…
(def pi-stream (scale-stream (partial-sums (pi-summands 1)) 4))
A simple squaring function.
(defn square [s] (* s s))
Euler’s transform is a fancy way to average three terms of a +/- sequence. It can dramatically increase the precision of an otherwise slowly converging sequence. The euler-transform function produces a sequence of results from the pi-stream. Again, note the constraint I place on BigDecimal precision.
(defn euler-transform [s] (let [s0 (nth s 0) s1 (nth s 1) s2 (nth s 2)] (lazy-seq (cons (with-precision 300 (- s2 (/ (square (- s2 s1)) (+ s0 (* -2 s1) s2)))) (euler-transform (next s))))))
It turns out that you can use Euler’s transform on it’s own results to further increase precision. The make-tableau function arranges the reapplication of a transform so that it can be repeatedly applied to itself.
(defn make-tableau [transform s] (lazy-seq (cons s (make-tableau transform (transform s)))))
And now we grab the first element of each repetition.
(defn accelerated-sequence [transform s] (map first (make-tableau transform s)))
A little utility to print the elements of a list.
(defn print-list [s] (doseq [n s] (println n)))
And now finally we print the first 150 iterations of the accelerated euler-transformation of the pi sequence.
(print-list (take 150 (accelerated-sequence euler-transform pi-stream)))
The 150th result (the first two lines of which are accurate) 3.14159265358979323846264338327950288419716939 93751058209749445923078164062862089986280348253 25256352154500957038098311477352066337027814620 46050836612955894790973622135628850706785335611 34661971615744147647325427374216175843967838084 09640502874355811103562027043480098593797655792 36344321165429010895M
Comments
Michael Feathers 31 minutes later:
I love laziness. A while back, I wrote a Mandelbrot Set generator in Haskell. The computation I formed was over an infinite grid ([0..] in x and y). Afterward, I used this function to “pluck” out interesting areas. With laziness, these areas were computed on demand:
window :: (Int, Int) -> (Int, Int) -> [[a]] -> [[a]] window (x0, x1) (y0, y1) = range y0 y1 . map (range x0 x1) where range m n = take n . drop m
This is a form of abstraction which is harder to achieve in non-lazy languages. You form a general computation and impose the edges later.
Tyler Jennings about 8 hours later:
You inspired me to try this in Scala. I had some trouble reading the prefix notation, but I think I figured it out. My result for the 150th value in the sequence is close to yours, but not identical. May have to do with the type of rounding I’m doing on the division.
Doesn’t feel as nice as your clojure version, but I’m a Scala newbie. There is probably a much better way.
Jeff about 22 hours later:
Instead of using lazy-seq “manually” you could use some of the build in Clojure functions. For example, you can write pi-summands using the standard iterate function (in a slightly evil feeling way and I’m sure there’s a better way to do this!)
(def pi-summands (map first (iterate (fn [[x y]] (let [n (+ y 2)] [(/ 1 n) n])) [1 1])))
Clojure also has some nifty short-hand syntax
(defn scale-stream [stream factor] (map #(* % factor) stream))
Where # indicates an anonymous function of one argument (represented by %).
Tyler Jennings 1 day later:
I thought I’d attempt a clone of your algorithm in Scala, since I’m learning that. Here’s my best result so far:
Still some rough edges in 2.7.5 I felt compelled to work around (BigDecimal support, no cons(::) operator for streams. The Clojure version feels cleaner to me – even though I’m not comfortable with prefix notation.
All the type definitions are required, of course. However, I feel like they’re more clutter than useful to the developer.
piyo 1 day later:
Great stuff. I learned something new.
@Jeff:You missed the part where the summands are negative on every other term. Uncle Bob’s explicit lazy-seq is really novel with the “map -”!
@Uncle Bob: You can use the deconstructing bind (“Vector binding-exprs”) for your euler-transform instead of let:
(defn euler-transform [[s0 s1 s2 :as s]] ...)
Also I find I get a “Divide by zero error” and/or “Stack overflow error” if I try to calculate a lot of terms, like say 350, with your implementation. i.e.
(print-list (take 50 (drop 300 pi-terms)))
-> ... java.lang.ArithmeticException: Divide by zero
BTW, 300 terms yields 126 correct digits, while my implementation of theCalculating Pi using the van Wijngaarden transformation yields 296 digits after 300 terms, which comes closer to the limits set by the with-precision. (Actually the comparison is more like apples to oranges, because the latter is calculating more. Oh well.) Still, 126 digits ought to be enough for anybody (in this universe).
Heikki 1 day later:
The Leibniz series actually converges on Pi/4, but the second paragraph says it converges on Pi. It took me a while to understand what the multiplication was needed for.
However, it’s very nice to read and try to grasp lisp after 8 or so years from taking the SICP course.
Page updated
Google Sites
Report abuse