Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Estimating PI with the Gregory-Leibniz series for clojure101
;; 1
;; the nth element of the series
(defn pi-series [n]
(/ (* 4 (Math/pow -1 n)) (inc (* 2 n))))
;; lazy sequence giving the individual terms
(def pi1 (map #(pi-series %) (iterate inc 0)))
;; basically the same, but will use fractions (slower, but nice to look at)
(def pi2 (map #(/ (* 4 (expt -1 %)) (+ (* 2 %) 1)) (iterate inc 0)))
;; tools.pi> (time (double (reduce + (take 1000 pi1))))
;; "Elapsed time: 38.863614 msecs"
;; 3.140592653839794
;; tools.pi> (time (double (reduce + (take 1000 pi2))))
;; "Elapsed time: 5230.369837 msecs"
;; 3.140592653839793
;; (the fraction will be converted by double, but is inspectable)
;; (double (reduce + (take 650 pi)))
;; takes a while to converge, to 2 decimal points
;; this function is not tail recursive and as such will blow the stack for large values
(defn sums [i j]
(cond (zero? i) (reduce + (map #(pi-series %) (range i (+ j 1))))
:else
(/ (+ (sums (- i 1) j) (sums (- i 1) (+ j 1))) 2)))
;; this gets the job done, but it's slow as sin.
;; let's see if we can do something about that
;; 1
(defn pi-formular [n]
(/ (* 4 (Math/pow -1 n)) (inc (* 2 n))))
(defn pi-formular-slow [n]
(/ (* 4 (expt -1 n)) (+ (* 2 n) 1)))
(defn series [f] (map #(f %) (iterate inc 0)))
;; (take 10 (series pi-formular-slow))
;; (4 -4/3 4/5 -4/7 4/9 -4/11 4/13 -4/15 4/17 -4/19)
;; 2
;; first elemet will be whatever the initial element for iterate is
(defn sums []
(map first (iterate (fn [[a b]] [(+ (pi-formular b) a) (inc b)]) [4 1])))
(defn sums-slow []
(map first (iterate (fn [[a b]] [(+ (pi-formular-slow b) a) (inc b)]) [4 1])))
;; by dropping the first element, we do not have to initialize in a special way
(defn sums-slow []
(map first (drop 1 (iterate (fn [[a b]] [(+ (pi-formular-slow b) a) (inc b)]) [0 0]))))
;; by giving the series generator as an argument, we only need one function
(defn sums [f]
(map first (drop 1 (iterate (fn [[a b]] [(+ (f b) a) (inc b)]) [0 0]))))
;(take 10 (sums pi-formular-slow))
;; it would be better if sums accepted a sequence, when we have to do average dampning
(defn sums-seq [sequence]
(map first (drop 1 (iterate (fn [[a b]] [(+ (first b) a) (rest b)]) [0 sequence]))))
;;(take 10 (sums-seq (series pi-formular-slow)))
;; well, that seems to do the job
(deftest test-sums-seq
;; (4 -4/3 4/5 -4/7 4/9 -4/11 4/13 -4/15 4/17 -4/19)
(is (= (list 4
(+ 4 -4/3)
(+ 4 -4/3 4/5)
(+ 4 -4/3 4/5 -4/7)
(+ 4 -4/3 4/5 -4/7 4/9)
(+ 4 -4/3 4/5 -4/7 4/9 -4/11)
(+ 4 -4/3 4/5 -4/7 4/9 -4/11 4/13)
(+ 4 -4/3 4/5 -4/7 4/9 -4/11 4/13 -4/15)
(+ 4 -4/3 4/5 -4/7 4/9 -4/11 4/13 -4/15 4/17)
(+ 4 -4/3 4/5 -4/7 4/9 -4/11 4/13 -4/15 4/17 -4/19))
(take 10 (sums-seq (series pi-formular-slow)))
(take 10 (sums pi-formular-slow)))))
;; 3
(defn damping [seq]
(map #(/ % 2) seq))
;; not the correct way to do this
(defn damping [sequence]
(map first (drop 1 (iterate (fn [[a b]] [(/ (+ (first b) (second b)) 2) (rest b)]) [0 sequence]))))
(deftest test-sums-seq
(is (= (list (/ (+ 4
(+ 4 -4/3)) 2)
(/ (+ (+ 4 -4/3)
(+ 4 -4/3 4/5)) 2)
(/ (+ (+ 4 -4/3 4/5)
(+ 4 -4/3 4/5 -4/7)) 2)
(/ (+ (+ 4 -4/3 4/5 -4/7)
(+ 4 -4/3 4/5 -4/7 4/9)) 2)
(/ (+ (+ 4 -4/3 4/5 -4/7 4/9)
(+ 4 -4/3 4/5 -4/7 4/9 -4/11)) 2)
(/ (+ (+ 4 -4/3 4/5 -4/7 4/9 -4/11)
(+ 4 -4/3 4/5 -4/7 4/9 -4/11 4/13)) 2)
(/ (+ (+ 4 -4/3 4/5 -4/7 4/9 -4/11 4/13)
(+ 4 -4/3 4/5 -4/7 4/9 -4/11 4/13 -4/15)) 2)
(/ (+ (+ 4 -4/3 4/5 -4/7 4/9 -4/11 4/13 -4/15)
(+ 4 -4/3 4/5 -4/7 4/9 -4/11 4/13 -4/15 4/17))2)
(/ (+ (+ 4 -4/3 4/5 -4/7 4/9 -4/11 4/13 -4/15 4/17)
(+ 4 -4/3 4/5 -4/7 4/9 -4/11 4/13 -4/15 4/17 -4/19))2))
'(10/3 46/15 334/105 982/315 10942/3465 140986/45045 28382/9009 2400458/765765 45788882/14549535)
(take 9 (damping (sums-seq (series pi-formular-slow))))
'(10/3 46/15 334/105 982/315 10942/3465 140986/45045 28382/9009 2400458/765765 45788882/14549535)
(take 9 (damping (sums pi-formular-slow)))
'(10/3 46/15 334/105 982/315 10942/3465 140986/45045 28382/9009 2400458/765765 45788882/14549535))))
;; 4
(double (first (damping (damping (damping (damping (damping (sums-seq (series pi-formular-slow)))))))))
(defn sums [i j]
(cond (zero? i) (reduce + (map #(pi-formular %) (range i (+ j 1))))
:else
(/ (+ (sums (- i 1) j) (sums (- i 1) (+ j 1))) 2)))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment