summaryrefslogtreecommitdiffstats
path: root/pi.scm
diff options
context:
space:
mode:
Diffstat (limited to 'pi.scm')
-rw-r--r--pi.scm165
1 files changed, 165 insertions, 0 deletions
diff --git a/pi.scm b/pi.scm
new file mode 100644
index 0000000..bb533bc
--- /dev/null
+++ b/pi.scm
@@ -0,0 +1,165 @@
+;; Copyright (C) 1991, 1993, 1994, 1995 Free Software Foundation, Inc.
+;;
+;; This program is free software; you can redistribute it and/or modify
+;; it under the terms of the GNU General Public License as published by
+;; the Free Software Foundation; either version 2, or (at your option)
+;; any later version.
+;;
+;; This program is distributed in the hope that it will be useful,
+;; but WITHOUT ANY WARRANTY; without even the implied warranty of
+;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+;; GNU General Public License for more details.
+;;
+;; You should have received a copy of the GNU General Public License
+;; along with this software; see the file COPYING. If not, write to
+;; the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+;;
+;; As a special exception, the Free Software Foundation gives permission
+;; for additional uses of the text contained in its release of GUILE.
+;;
+;; The exception is that, if you link the GUILE library with other files
+;; to produce an executable, this does not by itself cause the
+;; resulting executable to be covered by the GNU General Public License.
+;; Your use of that executable is in no way restricted on account of
+;; linking the GUILE library code into it.
+;;
+;; This exception does not however invalidate any other reasons why
+;; the executable file might be covered by the GNU General Public License.
+;;
+;; This exception applies only to the code released by the
+;; Free Software Foundation under the name GUILE. If you copy
+;; code from other Free Software Foundation releases into a copy of
+;; GUILE, as the General Public License permits, the exception does
+;; not apply to the code that you add in this way. To avoid misleading
+;; anyone as to the status of such modified files, you must delete
+;; this exception notice from them.
+;;
+;; If you write modifications of your own for GUILE, it is your choice
+;; whether to permit this exception to apply to your modifications.
+;; If you do not wish that, delete this exception notice.
+
+;;;; "pi.scm", program for computing digits of numerical value of PI.
+;;;; "bigpi.scm", program for computing digits of numerical value of PI.
+;;;; "e.scm", program for computing digits of numerical value of 'e'.
+;;; Authors: Aubrey Jaffer & Jerry D. Hedden
+
+;;; (pi <n> <d>) prints out <n> digits of pi in groups of <d> digits.
+
+;;; 'Spigot' algorithm origionally due to Stanly Rabinowitz.
+;;; This algorithm takes time proportional to the square of <n>/<d>.
+;;; This fact can make comparisons of computational speed between systems
+;;; of vastly differring performances quicker and more accurate.
+
+;;; Try (pi 100 5)
+;;; The digit size <d> will have to be reduced for larger <n> or an
+;;; overflow error will occur (on systems lacking bignums).
+
+;;; It your Scheme has bignums try (pi 1000).
+
+(define (pi n . args)
+ (if (null? args) (bigpi n)
+ (let* ((d (car args))
+ (r (do ((s 1 (* 10 s))
+ (i d (- i 1)))
+ ((zero? i) s)))
+ (n (+ (quotient n d) 1))
+ (m (quotient (* n d 3322) 1000))
+ (a (make-vector (+ 1 m) 2)))
+ (vector-set! a m 4)
+ (do ((j 1 (+ 1 j))
+ (q 0 0)
+ (b 2 (remainder q r)))
+ ((> j n))
+ (do ((k m (- k 1)))
+ ((zero? k))
+ (set! q (+ q (* (vector-ref a k) r)))
+ (let ((t (+ 1 (* 2 k))))
+ (vector-set! a k (remainder q t))
+ (set! q (* k (quotient q t)))))
+ (let ((s (number->string (+ b (quotient q r)))))
+ (do ((l (string-length s) (+ 1 l)))
+ ((>= l d) (display s))
+ (display #\0)))
+ (if (zero? (modulo j 10)) (newline) (display #\ )))
+ (newline))))
+
+;;; (pi <n>) prints out <n> digits of pi.
+
+;;; 'Spigot' algorithm originally due to Stanly Rabinowitz:
+;;;
+;;; PI = 2+(1/3)*(2+(2/5)*(2+(3/7)*(2+ ... *(2+(k/(2k+1))*(4)) ... )))
+;;;
+;;; where 'k' is approximately equal to the desired precision of 'n'
+;;; places times 'log2(10)'.
+;;;
+;;; This version takes advantage of "bignums" in SCM to compute all
+;;; of the requested digits in one pass! Basically, it calculates
+;;; the truncated portion of (PI * 10^n), and then displays it in a
+;;; nice format.
+
+(define (bigpi digits)
+ (let* ((n (* 10 (quotient (+ digits 9) 10))) ; digits in multiples of 10
+ (z (inexact->exact (truncate ; z = number of terms
+ (/ (* n (log 10)) (log 2)))))
+ (q (do ((x 2 (* 10000000000 x)) ; q = 2 * 10^n
+ (i (/ n 10) (- i 1)))
+ ((zero? i) x)))
+ (_pi (number->string ; _pi = PI * 10^n
+ ;; do the calculations in one pass!!!
+ (let pi_calc ((j z) (k (+ z z 1)) (p (+ q q)))
+ (if (zero? j)
+ p
+ (pi_calc (- j 1) (- k 2) (+ q (quotient (* p j) k))))))))
+ ;; print out the result ("3." followed by 5 groups of 10 digits per line)
+ (display (substring _pi 0 1)) (display #\.) (newline)
+ (do ((i 0 (+ i 10)))
+ ((>= i n))
+ (display (substring _pi (+ i 1) (+ i 11)))
+ (display (if (zero? (modulo (+ i 10) 50)) #\newline #\ )))
+ (if (not (zero? (modulo n 50))) (newline))))
+
+;;; (e <n>) prints out <n> digits of 'e'.
+
+;;; Uses the formula:
+;;;
+;;; 1 1 1 1 1
+;;; e = 1 + -- + -- + -- + -- + ... + --
+;;; 1! 2! 3! 4! k!
+;;;
+;;; where 'k' is determined using the desired precision 'n' in:
+;;;
+;;; n < ((k * (ln(k) - 1)) / ln(10))
+;;;
+;;; which uses Stirling's formula for approximating ln(k!)
+;;;
+;;; This program takes advantage of "bignums" in SCM to compute all
+;;; the requested digits at once! Basically, it calculates the
+;;; fractional part of 'e' (i.e., e-2) as a fraction of two bignums
+;;; 'e_n' and 'e_d', determines the integer part of (e_n * 10^n)/e_d,
+;;; and then displays it in a nice format.
+
+(define (e digits)
+ (let* ((n (* 10 (quotient (+ digits 9) 10))) ; digits in multiples of 10
+ (k (do ((i 15 (+ i 1))) ; k = number of terms
+ ((< n (/ (* i (- (log i) 1)) (log 10))) i)))
+ (q (do ((x 1 (* 10000000000 x)) ; q = 10^n
+ (i (/ n 10) (- i 1)))
+ ((zero? i) x)))
+ (_e (let ((ee
+ ; do calculations
+ (let e_calc ((i k) (e_d 1) (e_n 0))
+ (if (= i 1)
+ (cons (* q e_n) e_d)
+ (e_calc (- i 1) (* e_d i) (+ e_n e_d))))))
+ (number->string (+ (quotient (car ee) (cdr ee))
+ ; rounding
+ (if (< (remainder (car ee) (cdr ee))
+ (quotient (cdr ee) 2))
+ 0 1))))))
+ ;; print out the result ("2." followed by 5 groups of 10 digits per line)
+ (display "2.") (newline)
+ (do ((i 0 (+ i 10)))
+ ((>= i n))
+ (display (substring _e i (+ i 10)))
+ (display (if (zero? (modulo (+ i 10) 50)) #\newline #\ )))
+ (if (not (zero? (modulo n 50))) (newline))))