;;; PI -- Compute PI using bignums.

; See http://mathworld.wolfram.com/Pi.html for the various algorithms.

(library (rnrs-benchmarks pi)
  (export main)
  (import (rnrs) (rnrs r5rs) (rnrs-benchmarks))

  ; Utilities.
  
  (define (width x)
    (let loop ((i 0) (n 1))
      (if (< x n) i (loop (+ i 1) (* n 2)))))
  
  (define (root x y)
    (let loop ((g (expt
                   2
                   (quotient (+ (width x) (- y 1)) y))))
      (let ((a (expt g (- y 1))))
        (let ((b (* a y)))
          (let ((c (* a (- y 1))))
            (let ((d (quotient (+ x (* g c)) b)))
              (if (< d g) (loop d) g)))))))
  
  (define (square-root x)
    (root x 2))
  
  (define (quartic-root x)
    (root x 4))
  
  (define (square x)
    (* x x))
  
  ; Compute pi using the 'brent-salamin' method.
  
  (define (pi-brent-salamin nb-digits)
    (let ((one (expt 10 nb-digits)))
      (let loop ((a one)
                 (b (square-root (quotient (square one) 2)))
                 (t (quotient one 4))
                 (x 1))
        (if (= a b)
            (quotient (square (+ a b)) (* 4 t))
            (let ((new-a (quotient (+ a b) 2)))
              (loop new-a
                    (square-root (* a b))
                    (- t
                              (quotient
                               (* x (square (- new-a a)))
                               one))
                    (* 2 x)))))))
  
  ; Compute pi using the quadratically converging 'borwein' method.
  
  (define (pi-borwein2 nb-digits)
    (let* ((one (expt 10 nb-digits))
           (one^2 (square one))
           (one^4 (square one^2))
           (sqrt2 (square-root (* one^2 2)))
           (qurt2 (quartic-root (* one^4 2))))
      (let loop ((x (quotient
                     (* one (+ sqrt2 one))
                     (* 2 qurt2)))
                 (y qurt2)
                 (p (+ (* 2 one) sqrt2)))
        (let ((new-p (quotient (* p (+ x one))
                                      (+ y one))))
          (if (= x one)
              new-p
              (let ((sqrt-x (square-root (* one x))))
                (loop (quotient
                       (* one (+ x one))
                       (* 2 sqrt-x))
                      (quotient
                       (* one (+ (* x y) one^2))
                       (* (+ y one) sqrt-x))
                      new-p)))))))
  
  ; Compute pi using the quartically converging 'borwein' method.
  
  (define (pi-borwein4 nb-digits)
    (let* ((one (expt 10 nb-digits))
           (one^2 (square one))
           (one^4 (square one^2))
           (sqrt2 (square-root (* one^2 2))))
      (let loop ((y (- sqrt2 one))
                 (a (- (* 6 one) (* 4 sqrt2)))
                 (x 8))
        (if (= y 0)
            (quotient one^2 a)
            (let* ((t1 (quartic-root (- one^4 (square (square y)))))
                   (t2 (quotient
                        (* one (- one t1))
                        (+ one t1)))
                   (t3 (quotient
                        (square (quotient (square (+ one t2)) one))
                        one))
                   (t4 (+ one
                                 (+ t2
                                           (quotient (square t2) one)))))
              (loop t2
                    (quotient
                     (- (* t3 a) (* x (* t2 t4)))
                     one)
                    (* 4 x)))))))
  
  ; Try it.
  
  (define (pies n m s)
    (if (< m n)
        '()
        (let ((bs (pi-brent-salamin n))
              (b2 (pi-borwein2 n))
              (b4 (pi-borwein4 n)))
          (cons (list b2 (- bs b2) (- b4 b2))
                (pies (+ n s) m s)))))
  
  (define expected
    '((314159265358979323846264338327950288419716939937507
       -54
       124)
      (31415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170673
       -51
       -417)
      (3141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408122
       -57
       -819)
      (314159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038195
       -76
       332)
      (31415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019089
       -83
       477)
      (3141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141268
       -72
       -2981)
      (314159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536431
       -70
       -2065)
      (31415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116089
       -79
       1687)
      (3141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118542
       -92
       -2728)
      (314159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194907
       -76
       -3726)))
  
  (define (main . args)
    (run-benchmark
     "pi"
     pi-iters
     (lambda (result) (equal? result expected))
     (lambda (n m s) (lambda () (pies n m s)))
     50
     500
     50)))