Lazy sequence of the digits of pi

piseq.clj

; based on http://bellard.org/pi/pi.c
(defn inv-mod
  "Returns the inverse of (mod x y)"
  ([x y] (inv-mod y x y 1 0))
  ([y u v c a]
     (let [q  (long (Math/floor (/ v u)))
           c' (- a (* q c))
           u' (- v (* q u))]
       (if (zero? u')
         (let [m (mod c y)]
           (if (< m 0) (+ m y) m))
         (recur y u' u c' c)))))

(defn mul-mod
  "Returns (a*b) mod m"
  [a b m]
  (mod (* a b) m))

(defn pow-mod
  "Returns (a^b) mod m"
  ([a b m] (pow-mod a b m 1))
  ([a b m r]
     (let [r' (if (zero? (bit-and b 1)) r (mul-mod r a m))
           b' (bit-shift-right b 1)]
       (if (zero? b')
         r'
         (recur (mul-mod a a m) b' m r')))))

(defn prime? 
  "True if n is prime"
  [n]
  (if (zero? (mod n 2))
    false
    (loop [i 3 r (long (Math/sqrt n))]
      (if (zero? (mod n i))
        false
        (if (> i r)
          true
          (recur (+ i 2) r))))))

(defn next-prime
  "Returns the next prime number immediately after n"
  [n]
  (let [n2 (inc n)]
    (if (prime? n2) n2 (recur n2))))

(defn digits [n]
  (let [N (int (* (+ n 20) (/ (Math/log 10) (Math/log 2))))]
    (loop [sum 0 a 3]
      (let [vmax (int (/ (Math/log (* 2 N)) (Math/log a)))
            av   (reduce * 1 (repeat vmax a))
            s    (loop [k 1, s 0, num 1, den 1, kq 1, kq2 1]
                   (if (<= k N)
                     (let [;; numerator
                           [num-t num-v kq']
                           (if (>= kq a)
                             (loop [t k, v 0]
                               (let [t' (long (/ t a))
                                     v' (dec v)]
                                 (if (zero? (mod t' a))
                                   (recur t' v')
                                   [t' v' 0])))
                             [k 0 (inc kq)])
                           num' (mul-mod num num-t av)
                           ;; denominator
                           [den-t den-v kq2']
                           (let [t (dec (* 2 k))]
                             (if (>= kq2 a)
                               (if (= kq2 a)
                                 (loop [t' t, v 0]
                                   (if (zero? (mod t' a))
                                     (recur (long (/ t' a)) (inc v))
                                     [t' v (+ 2 kq2)]))
                                 [t 0 (+ 2 (- kq2 a))])
                               [t 0 (+ 2 kq2)]))
                           den' (mul-mod den den-t av)]
                       (let [v (+ num-v den-v)]
                         (recur (inc k)
                                (if (> v 0)
                                  (let [s' (+ s
                                              (loop [i v, t (mul-mod (mul-mod (inv-mod den' av) num' av) k av)]
                                                (if (< i vmax)
                                                  (recur (inc i) (mul-mod t a av))
                                                  t)))]
                                    (if (>= s' av) (- s' av) s'))
                                  s)
                                num'
                                den'
                                kq'
                                kq2')))
                     s))
            sum' (mod (+ sum (/ (mul-mod s (pow-mod 10 (dec n) av) av) av)) 1.0)]
        (if (<= a (* 2 N))
          (recur sum' (next-prime a))
          sum)))))