;;;
;;; Computation of the n'th digit of pi in any base 10
;;;
;;; Jens Axel Søgaard, 27. maj 2001, a dull sunday
;;;
;;; Reference: "Computation of the n'th digit of pi in any base in O(n^2)"
;;;            by Fabrice Bellard, 1997
;;;            <http://www.stud.enst.fr/~bellard/pi/pi_n2/pi_n2.html>

;;
;;  MAIN
;;

;; returns the n'th digit of pi in base 10

(define (pi-digit n)
  (let* ((N (inexact->exact (truncate (* (+ n 20) (/ (log 10) (log 2)))))))
    (letrec ((for-p
              (lambda (p sum)
                ;; termination condition
                (if (> p (* 2 N))
                    sum
                    (let* ((vmax    (inexact->exact (truncate (/ (log (* 2 N)) (log p)))))
                           (modulus (expt p vmax)))
                      (letrec ((for-k
                                (lambda (k b A v alpha)
                                  (if (> k N)
                                      alpha
                                      (let* ((v (+ (- v (multiplicity p k))
                                                   (multiplicity p (- (* 2 k) 1))))
                                             (A  (modulo (/ (* (- (* 2 k) 1)
                                                               A)
                                                            (expt p (multiplicity p (- (* 2 k) 1))))
                                                         modulus))
                                             (b (modulo (* b (/ k (expt p (multiplicity p k)))) modulus)))

                                        (for-k (+ k 1) b A v
                                               (if (> v 0)
                                                   (+ alpha (* k b (inv-mod A modulus) (expt p (- vmax v))))
                                                   alpha)))))))
                        (for-p (next-prime p)
                               (mod1 (+ sum
                                        (/ (modulo (* (expt-mod 10 (- n 1) modulus)
                                                      (for-k 1 1 1 0 0))
                                                   modulus)
                                           modulus))))))))))
      (first-digit (for-p 3 0)))))

;;
;; Arithmetical utilities
;;


;; Extended version of Euclids algoritm for finding the gcd
;; returns
;;        (g,a,b)  s.t.  g=gcd(m,n)  and  am+bm=g

(define (gcd-ext m n)
  (letrec ((helper (lambda (m n am bm an bn)
                     (if (> m n)
                         (helper (- m n) n (- am an) (- bm bn) an bn)
                         (if (> n m)
                             (helper m (- n m) am bm (- an am) (- bn bm))
                             (cons m (cons am (cons bm '()))))))))
    (helper m n 1 0 0 1)))


;;  b  s.t.  ab = 1 mod n    ( note: requires gcd(a,n)=1 )
(define inv-mod
  (lambda (a n)
    (cadr (gcd-ext a n))))


;;  base ^ exp  mod    (Abelson et al. p.51)
(define (expt-mod base exp m)
  (cond ((= exp 0) 1)
        ((even? exp)
         (remainder (square (expt-mod base (/ exp 2) m))
                    m))
        (else
         (remainder (* base (expt-mod base (- exp 1) m))
                    m))))

(define (square x) (* x x))


;;;
;;; PRIME RELATED  (Abelson et al. p.50)
;;;

;;; TODO: Replace prime? from O(sqrt(n)) to O(log n) algorithm

(define (smallest-divisor n)
  (find-divisor n 2))

(define (find-divisor n test-divisor)
  (cond ((> (square test-divisor) n) n)
        ((divides? test-divisor n) test-divisor)
        (else (find-divisor n (+ test-divisor 1)))))

(define (divides? a b)
  (= (remainder b a) 0))

(define (prime? n)
  (= n (smallest-divisor n)))

(define (next-prime n)
  (if (= n 2)
      3
      (if (prime? (+ n 2))
          (+ n 2)
          (next-prime (+ n 2)))))

;;
;; FACTORISE
;;

;; Example:  (factor 45) -> ((3 2) (5 1))

(define (factor n)
  (factor-from n 2))

(define (factor-from  n f)
  (if (= n 1)
      '()
      (if (divides? f n)
          (cons (list f (multiplicity f n))
                (factor-from (/ n (expt f (multiplicity f n))) (next-prime f)))
          (factor-from n (next-prime f)))))

(define (multiplicity p n)
  (if (divides? p n)
      (+ 1 (multiplicity p (/ n p)))
      0))


;; Utilies

(define (mod1 x)
  (- x (floor x)))

(define (first-digit x)
  (inexact->exact (floor (* 10 (mod1 x)))))


;;
;; TEST
;;

(define (interval from to)
  (if (> from to)
      '()
      (cons from (interval (+ 1 from) to))))


;; guile> (map pi-digit (interval 1 100))

;; (1 4 1 5 9 2 6 5 3 5
;;  8 9 7 9 3 2 3 8 4 6
;;  2 6 4 3 3 8 3 2 7 9
;;  5 0 2 8 8 4 1 9 7 1
;;  6 9 3 9 9 3 7 5 1 0
;;  5 8 2 0 9 7 4 9 4 4
;;  5 9 2 3 0 7 8 1 6 4
;;  0 6 2 8 6 2 0 8 9 9
;;  8 6 2 8 0 3 4 8 2 5
;;  3 4 2 1 1 7 0 6 7 9)