;;; pi-ratios.scm (ported from ratios.lisp, from CMUCL's cl-bench)
;;;   -- calculate digits of pi using ratios
;;
;; Time-stamp: <2003-12-29 emarsden>
;;
;;
;; This code was posted to comp.lang.lisp on 2001-12-28 by Vladimir
;; Nesterovsky <vnestr@netvision.net.il>, in message
;; <tcio2ucso5eace5dc7cnbevpc2nqukgofh@4ax.com>. 
;;
;; "I played with some code to calculate more then 16 digits of the 
;; number pi. The simplest series for pi/4 is 1 - 1/3 + 1/5 - 1/7 ... 
;; with Euler transformation applied on its partial sums series 
;; several times (as per SICP).
;; 
;; Of course denominators will grow very fast if we'd just
;; perform calculations in a straightforward manner which would
;; make them very slow, so I tried to "adjust" the ratio to having 
;; no more precision than I need (say 1000 digits), which greatly
;; improved the speed and made it possible to calculate much more
;; digits of pi in the same amount of time."
;;
;;
;; sample calls (timings on PIII/550/Win98):
;;   make 20 euler transforms, pre-generate 41 terms 
;;   (2n+1 by default because each transform shortens the list by two)
;; 	(find-pi 20)     ; 0.3 sec CLISP // 5 sec ACL 
;;   make 20 euler transforms, pre-generate 100 terms
;; 	(find-pi 20 100) ; 0.3 sec CLISP // 6 sec ACL
;;   make 20 euler transforms, pre-generate 200 terms
;; 	(find-pi 20 200) ; 0.3 sec // 7 sec
;;   50 transforms, 101 pre-generated terms
;; 	(find-pi 50)     ; 2 sec // 37 sec (35 compiled) 
;;   1000 pre-generated (working only on 101 last) terms for 120 digits
;; 	(find-pi 50 1000 120)     ; 2 sec // 84 sec (80 compiled)
;;   10,000 pre-generated terms for better precision for 150 digits
;; 	(find-pi 50 10000 150)    ; 4.5 sec // ??
;;   find 1040 digits of the number pi
;; 	(find-pi 120 150000 1039) ; 6.5 min, CLISP // ?? ACL

;;;; Comment this out if you want to test in another Scheme
(use numbers)

;; To get the "error" procedure...
(require-extension (srfi 23))

;; Not needed for Chicken, but useful for comparison with other implementations
(define (assert x)
  (if (not x) (error "assertion failed")))

;; Compat
(define-syntax dotimes
  (syntax-rules ()
    ((dotimes (var times) body? ...)
     (do ((var 0 (+ var 1)))
         ((= var times))
       body? ...))))

;; Also not needed for Chicken
(define-syntax do*
  (syntax-rules ()
    ((do* ((var0 init0 step0) ...)
          (test expr0 ...)
          command0 ...)
     (let* ((var0 init0) ...)
       (let loop ((var0 var0) ...)
         (if test
             (begin expr0 ...)
             (begin command0 ...
                    (set! var0 step0) ...
                    (loop var0 ...))))))))

(define (last l) (if (null? (cdr l)) l (last (cdr l))))

(define *report-digits* 55)        ; digits to report
(define *interim-precision* 95)    ; decimal places of interim precision
(define *supress-interim-printout* #t)
(define *supress-interim-dot-printout* #t)

(define (adjust-ratio r n)
  "adjust ratio's precision to n decimal places"
  (let ((base (expt 10 n))
        (rem (- r (round r))))
    (if (< (denominator rem) base)
        r
        (/ (round (* r base)) base))))

;; a series for pi/4 = 
;;   = 1 - 1/3 + 1/5 - 1/7 + 1/9 - 1/11 + ....
(define (s-pi/4 m len-from-end)
  "generate up to m partial sums for pi series"
  (do* ((i0 (- m len-from-end) i0)
        (i  0   (+ i 1))
        (s  0   (adjust-ratio (+ s xi) *interim-precision*))
        (j  1   (- j))
        (n  1   (+ n 2))
        (xi 1   (adjust-ratio (/ j n) *interim-precision*))
        (sums '()  (if (> i i0) (cons s sums) sums)))
       ((>= i m)  (reverse sums))
       (if (and (not *supress-interim-dot-printout*)
                (zero? (remainder i 1000))
                (> i 0))
           (display "."))))

;; the Euler transformation of a partial sums series of alternating series
;; A[n] = S[n+1] - (S[n+1]-S[n])^2/(S[n-1]-2S[n]+S[n+1])
;;   returns (cons 'div-zero seq) on 0/0 (to stop further transformations)
(define (euler-trans seq) 
  "perform a euler transformation on a sequence"
  (let ((a1 (car seq))
        (a2 (cadr seq)))
    (call-with-current-continuation
     (lambda (return)
       (map (lambda(a3)
              (let ((den (+ a1 (* -2 a2) a3)))
                (if (zero? den)        ; sould I use throw/catch here?
                    (return (cons 'div-zero seq)))
                (let ((r (- a3 (/ (sqr (- a3 a2)) den))))
                  (set! a1 a2)
                  (set! a2 a3)
                  (adjust-ratio r *interim-precision*))))
            (cddr seq))))))

;; main function to call
(define (find-pi n . rest)
   "(find-pi n [m d p]) to find d digits of the pi value by performing n euler
   transformations on m terms of pi series with p digits of interim precision."
  ;; For a lack of #!optional args that can refer back to preceding args, we
  ;; have to this ugly hack :(
  (let ((m #f)
        (dig #f)
        (prec #f)
        (sup #f))
    (if (null? rest)
        (set! m (+ (* 2 n) 1))
        (begin (set! m (car rest)) (set! rest (cdr rest))))
    (if (null? rest)
        (set! dig *report-digits*)
        (begin (set! dig (car rest)) (set! rest (cdr rest))))
    (if (null? rest)
        (set! prec (+ dig 40))
        (begin (set! prec (car rest)) (set! rest (cdr rest))))
    (if (null? rest)
        (set! sup *supress-interim-printout*)
        (set! sup (car rest)))
    (let ((*report-digits* dig)
          (*interim-precision* prec)
          (*supress-interim-printout* sup))
      (if (< m (+ (* 2 n) 1))
          (set! m (+ (* 2 n) 1)))
      ;; (format t "~&Generating ~A terms of pi/4 series " m)
      (do ((s   (s-pi/4 m (+ (* 2 n) 1))
                (euler-trans s))
           (i   1   (+ i 1)))
          ((or (> i n) (eq? 'div-zero (car s)))
           ;; (format t "~& After ~A Euler transforms:~%" (- i 1))
           (n-digits (* 4 (car (last s)))
                     *report-digits*))
        (if (not *supress-interim-printout*)
            (begin (newline)
                   (display (last-n-digits
                             (n-digits (* 4 (car (last s)))
                                       *report-digits*) 
                             79))))))))

(define (n-digits v n)
  "return v with n digits after decimal point as integer"
  (round (* v (expt 10 n))))

(define (last-n-digits v n)
  "return last n digits of long integer as integer"
  (floor (- v (remainder v (expt 10 n)))))

(define (sqr x)
  "find the square of number"
  (* x x))

;; entry point
(define (run-pi-ratios)
  (assert (eqv? (find-pi 20 1000 78)
                (find-pi 90 0 78 81))))


;;(find-pi 20 100) ; 55 digits
;;31415926535897932384626433832795028841971693993751058210
;;(find-pi 20 1000 78) ; 78 digits, 1000 terms
;;3141592653589793238462643383279502884197169399375105820974944592307816406286209
;;(find-pi 90 0 78 81) ; 181 terms --
;;       ; like SICP's "tableau" of streams --
;;       ; all in all slower because it produces bigger interim lists to be able
;;       ; to run much more transforms (although it exits in the middle on 0/0).
;;3141592653589793238462643383279502884197169399375105820974944592307816406286209

(display "run-pi-ratios:")
(newline)
(time (dotimes (_ 2) (run-pi-ratios)))
;; EOF