141 lines
3.9 KiB
Scheme
141 lines
3.9 KiB
Scheme
;;; PI -- Compute PI using bignums.
|
|
|
|
; See http://mathworld.wolfram.com/Pi.html for the various algorithms.
|
|
|
|
(import (scheme base)
|
|
(scheme write)
|
|
(scheme read))
|
|
|
|
(define div quotient)
|
|
|
|
; 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
|
|
(div (+ (width x) (- y 1)) y))))
|
|
(let ((a (expt g (- y 1))))
|
|
(let ((b (* a y)))
|
|
(let ((c (* a (- y 1))))
|
|
(let ((d (div (+ 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 (div (square one) 2)))
|
|
(t (div one 4))
|
|
(x 1))
|
|
(if (= a b)
|
|
(div (square (+ a b)) (* 4 t))
|
|
(let ((new-a (div (+ a b) 2)))
|
|
(loop new-a
|
|
(square-root (* a b))
|
|
(- t
|
|
(div
|
|
(* 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 (div
|
|
(* one (+ sqrt2 one))
|
|
(* 2 qurt2)))
|
|
(y qurt2)
|
|
(p (+ (* 2 one) sqrt2)))
|
|
(let ((new-p (div (* p (+ x one))
|
|
(+ y one))))
|
|
(if (= x one)
|
|
new-p
|
|
(let ((sqrt-x (square-root (* one x))))
|
|
(loop (div
|
|
(* one (+ x one))
|
|
(* 2 sqrt-x))
|
|
(div
|
|
(* 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)
|
|
(div one^2 a)
|
|
(let* ((t1 (quartic-root (- one^4 (square (square y)))))
|
|
(t2 (div
|
|
(* one (- one t1))
|
|
(+ one t1)))
|
|
(t3 (div
|
|
(square (div (square (+ one t2)) one))
|
|
one))
|
|
(t4 (+ one
|
|
(+ t2
|
|
(div (square t2) one)))))
|
|
(loop t2
|
|
(div
|
|
(- (* 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 (main)
|
|
(let* ((count (read))
|
|
(input1 (read))
|
|
(input2 (read))
|
|
(input3 (read))
|
|
(output (read))
|
|
(s4 (number->string count))
|
|
(s3 (number->string input3))
|
|
(s2 (number->string input2))
|
|
(s1 (number->string input1))
|
|
(name "pi"))
|
|
(run-r7rs-benchmark
|
|
(string-append name ":" s1 ":" s2 ":" s3 ":" s4)
|
|
count
|
|
(lambda ()
|
|
(pies (hide count input1)
|
|
(hide count input2)
|
|
(hide count input3)))
|
|
(lambda (result) (equal? result output)))))
|
|
|
|
(include "src/common.sch")
|