home *** CD-ROM | disk | FTP | other *** search
/ Dream 44 / Amiga_Dream_44.iso / RiscPc / programmation / scm4e2.arc / !Scm / scm / pi < prev    next >
Text File  |  1994-10-28  |  5KB  |  133 lines
  1. ;;;; "pi.scm", program for computing digits of numerical value of PI.
  2. ;;; Copyright (C) 1991 Aubrey Jaffer.
  3. ;;; See the file `COPYING' for terms applying to this program.
  4.  
  5. ;;; (pi <n> <d>) prints out <n> digits of pi in groups of <d> digits.
  6.  
  7. ;;; 'Spigot' algorithm origionally due to Stanly Rabinowitz.
  8. ;;; This algorithm takes time proportional to the square of <n>/<d>.
  9. ;;; This fact can make comparisons of computational speed between systems
  10. ;;; of vastly differring performances quicker and more accurate.
  11.  
  12. ;;; Try (pi 100 5)
  13. ;;; The digit size <d> will have to be reduced for larger <n> or an
  14. ;;; overflow error will occur (on systems lacking bignums).
  15.  
  16. ;;; It your Scheme has bignums try (pi 1000).
  17.  
  18. (define (pi n . args)
  19.   (if (null? args) (bigpi n)
  20.       (let* ((d (car args))
  21.          (r (do ((s 1 (* 10 s))
  22.              (i d (- i 1)))
  23.             ((zero? i) s)))
  24.          (n (+ (quotient n d) 1))
  25.          (m (quotient (* n d 3322) 1000))
  26.          (a (make-vector (+ 1 m) 2)))
  27.     (vector-set! a m 4)
  28.     (do ((j 1 (+ 1 j))
  29.          (q 0 0)
  30.          (b 2 (remainder q r)))
  31.         ((> j n))
  32.       (do ((k m (- k 1)))
  33.           ((zero? k))
  34.         (set! q (+ q (* (vector-ref a k) r)))
  35.         (let ((t (+ 1 (* 2 k))))
  36.           (vector-set! a k (remainder q t))
  37.           (set! q (* k (quotient q t)))))
  38.       (let ((s (number->string (+ b (quotient q r)))))
  39.         (do ((l (string-length s) (+ 1 l)))
  40.         ((>= l d) (display s))
  41.           (display #\0)))
  42.       (if (zero? (modulo j 10)) (newline) (display #\ )))
  43.     (newline))))
  44.  
  45. ;;; "bigpi.scm", program for computing digits of numerical value of PI.
  46. ;;; Copyright (C) 1993, 1994 Jerry D. Hedden
  47. ;;; See the file `COPYING' for terms applying to this program.
  48.  
  49. ;;; (pi <n>) prints out <n> digits of pi.
  50.  
  51. ;;; 'Spigot' algorithm originally due to Stanly Rabinowitz:
  52. ;;;
  53. ;;; PI = 2+(1/3)*(2+(2/5)*(2+(3/7)*(2+ ... *(2+(k/(2k+1))*(4)) ... )))
  54. ;;;
  55. ;;; where 'k' is approximately equal to the desired precision of 'n'
  56. ;;; places times 'log2(10)'.
  57. ;;;
  58. ;;; This version takes advantage of "bignums" in SCM to compute all
  59. ;;; of the requested digits in one pass!  Basically, it calculates
  60. ;;; the truncated portion of (PI * 10^n), and then displays it in a
  61. ;;; nice format.
  62.  
  63. (define (bigpi digits)
  64.   (let* ((n (* 10 (quotient (+ digits 9) 10)))    ; digits in multiples of 10
  65.      (z (inexact->exact (truncate        ; z = number of terms
  66.                  (/ (* n (log 10)) (log 2)))))
  67.      (q (do ((x 2 (* 10000000000 x))    ; q = 2 * 10^n
  68.          (i (/ n 10) (- i 1)))
  69.         ((zero? i)  x)))
  70.      (_pi (number->string            ; _pi = PI * 10^n
  71.            ;; do the calculations in one pass!!!
  72.            (let pi_calc ((j z) (k (+ z z 1)) (p (+ q q)))
  73.          (if (zero? j)
  74.              p
  75.              (pi_calc (- j 1) (- k 2) (+ q (quotient (* p j) k))))))))
  76.     ;; print out the result ("3." followed by 5 groups of 10 digits per line)
  77.     (display (substring _pi 0 1)) (display #\.) (newline)
  78.     (do ((i 0 (+ i 10)))
  79.     ((>= i n))
  80.       (display (substring _pi (+ i 1) (+ i 11)))
  81.       (display (if (zero? (modulo (+ i 10) 50)) #\newline #\ )))
  82.     (if (not (zero? (modulo n 50))) (newline))))
  83.  
  84. ;;; "e.scm", program for computing digits of numerical value of 'e'.
  85. ;;; Copyright (C) 1994 Jerry D. Hedden
  86. ;;; See the file 'COPYING' for terms applying to this program.
  87.  
  88. ;;; (e <n>) prints out <n> digits of 'e'.
  89.  
  90. ;;; Uses the formula:
  91. ;;;
  92. ;;;           1    1    1    1          1
  93. ;;;   e = 1 + -- + -- + -- + -- + ... + --
  94. ;;;           1!   2!   3!   4!         k!
  95. ;;;
  96. ;;; where 'k' is determined using the desired precision 'n' in:
  97. ;;;
  98. ;;;    n  <  ((k * (ln(k) - 1)) / ln(10))
  99. ;;;
  100. ;;; which uses Stirling's formula for approximating ln(k!)
  101. ;;;
  102. ;;; This program takes advantage of "bignums" in SCM to compute all
  103. ;;; the requested digits at once!  Basically, it calculates the
  104. ;;; fractional part of 'e' (i.e., e-2) as a fraction of two bignums
  105. ;;; 'e_n' and 'e_d', determines the integer part of (e_n * 10^n)/e_d,
  106. ;;; and then displays it in a nice format.
  107.  
  108. (define (e digits)
  109.   (let* ((n (* 10 (quotient (+ digits 9) 10)))    ; digits in multiples of 10
  110.      (k (do ((i 15 (+ i 1)))        ; k = number of terms
  111.         ((< n (/ (* i (- (log i) 1)) (log 10)))  i)))
  112.      (q (do ((x 1 (* 10000000000 x))    ; q = 10^n
  113.          (i (/ n 10) (- i 1)))
  114.         ((zero? i)  x)))
  115.      (_e (let ((ee
  116.             ; do calculations
  117.             (let e_calc ((i k) (e_d 1) (e_n 0))
  118.               (if (= i 1)
  119.               (cons (* q e_n) e_d)
  120.               (e_calc (- i 1) (* e_d i) (+ e_n e_d))))))
  121.            (number->string (+ (quotient (car ee) (cdr ee))
  122.                   ; rounding
  123.                   (if (< (remainder (car ee) (cdr ee))
  124.                      (quotient (cdr ee) 2))
  125.                       0 1))))))
  126.     ;; print out the result ("2." followed by 5 groups of 10 digits per line)
  127.     (display "2.") (newline)
  128.     (do ((i 0 (+ i 10)))
  129.     ((>= i n))
  130.       (display (substring _e i (+ i 10)))
  131.       (display (if (zero? (modulo (+ i 10) 50)) #\newline #\ )))
  132.     (if (not (zero? (modulo n 50))) (newline))))
  133.