home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / mitsch75.zip / scheme-7_5_17-src.zip / scheme-7.5.17 / src / runtime / random.scm < prev    next >
Text File  |  2000-04-10  |  8KB  |  222 lines
  1. #| -*-Scheme-*-
  2.  
  3. $Id: random.scm,v 14.23 2000/04/11 03:46:50 cph Exp $
  4.  
  5. Copyright (c) 1993-2000 Massachusetts Institute of Technology
  6.  
  7. This program is free software; you can redistribute it and/or modify
  8. it under the terms of the GNU General Public License as published by
  9. the Free Software Foundation; either version 2 of the License, or (at
  10. your option) any later version.
  11.  
  12. This program is distributed in the hope that it will be useful, but
  13. WITHOUT ANY WARRANTY; without even the implied warranty of
  14. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15. General Public License for more details.
  16.  
  17. You should have received a copy of the GNU General Public License
  18. along with this program; if not, write to the Free Software
  19. Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  20. |#
  21.  
  22. ;;;; Random Number Generator
  23. ;;; package: (runtime random-number)
  24.  
  25. (declare (usual-integrations))
  26.  
  27. ;;; A "subtract-with-borrow" RNG, based on the algorithm from "A New
  28. ;;; Class of Random Number Generators", George Marsaglia and Arif
  29. ;;; Zaman, The Annals of Applied Probability, Vol. 1, No. 3, 1991.
  30.  
  31. ;;; The basic algorithm produces a sequence of values x[n] by
  32. ;;;      let t = (x[n-s] - x[n-r] - borrow[n-1])
  33. ;;;      if (t >= 0)
  34. ;;;        then x[n] = t, borrow[n] = 0
  35. ;;;        else x[n] = t + b, borrow[n] = 1
  36. ;;; where the constants R, S, and B are chosen according to some
  37. ;;; constraints that are explained in the article.  Finding
  38. ;;; appropriate constants is compute-intensive; the constants used
  39. ;;; here are taken from the article and are claimed to represent
  40. ;;; "hundreds of hours" of compute time.  The period of this generator
  41. ;;; is (- (EXPT B R) (EXPT B S)), which is approximately (EXPT 10 414).
  42.  
  43. (define-integrable r 43)
  44. (define-integrable s 22)
  45. (define-integrable b 4294967291 #|(- (expt 2 32) 5)|#)
  46. (define-integrable b. 4294967291. #|(exact->inexact b)|#)
  47.  
  48. (define (random modulus #!optional state)
  49.   (let ((element
  50.      (flo:random-unit
  51.       (guarantee-random-state (if (default-object? state) #f state)
  52.                   'RANDOM))))
  53.     ;; Kludge: an exact integer modulus means that result is an exact
  54.     ;; integer.  Otherwise, the result is a real number.
  55.     (cond ((and (flo:flonum? modulus) (flo:< 0. modulus))
  56.        (flo:* element modulus))
  57.       ((and (int:integer? modulus) (int:< 0 modulus))
  58.        (flo:truncate->exact (flo:* element (int:->flonum modulus))))
  59.       ((and (real? modulus) (< 0 modulus))
  60.        (* (inexact->exact element) modulus))
  61.       (else
  62.        (error:wrong-type-argument modulus "positive real" 'RANDOM)))))
  63.  
  64. (define (flo:random-unit state)
  65.   (let ((mask (set-interrupt-enables! interrupt-mask/gc-ok)))
  66.     (let ((index (random-state-index state))
  67.       (vector (random-state-vector state)))
  68.       (let ((element (flo:vector-ref vector index)))
  69.     (let ((difference
  70.            (flo:- (flo:- (flo:vector-ref vector
  71.                          (if (fix:< (fix:- index s) 0)
  72.                          (fix:+ (fix:- index s) r)
  73.                          (fix:- index s)))
  74.                  element)
  75.               (random-state-borrow state))))
  76.       (if (flo:< difference 0.)
  77.           (begin
  78.         (flo:vector-set! vector index (flo:+ difference b.))
  79.         (set-random-state-borrow! state 1.))
  80.           (begin
  81.         (flo:vector-set! vector index difference)
  82.         (set-random-state-borrow! state 0.)))
  83.       (set-random-state-index! state
  84.                    (if (fix:= (fix:+ index 1) r)
  85.                        0
  86.                        (fix:+ index 1))))
  87.     (set-interrupt-enables! mask)
  88.     (flo:/ element b.)))))
  89.  
  90. (define (random-byte-vector n #!optional state)
  91.   (let ((state (if (default-object? state) #f state))
  92.     (s (make-string n)))
  93.     (do ((i 0 (fix:+ i 1)))
  94.     ((fix:= i n))
  95.       (vector-8b-set! s i (random 256 state)))
  96.     s))
  97.  
  98. (define (make-random-state #!optional state)
  99.   (let ((state (if (default-object? state) #f state)))
  100.     (if (or (eq? #t state) (int:integer? state))
  101.     ;; Use good random source if available
  102.     (if (file-exists? "/dev/urandom")
  103.         (call-with-input-file "/dev/urandom"
  104.           (lambda (port)
  105.         (initial-random-state
  106.          (lambda (b)
  107.            (let outer ()
  108.              (let inner
  109.              ((m #x100)
  110.               (n (char->integer (read-char port))))
  111.                (cond ((< m b)
  112.                   (inner (* m #x100)
  113.                      (+ (* n #x100)
  114.                     (char->integer (read-char port)))))
  115.                  ((< n b) n)
  116.                  (else (outer)))))))))
  117.         (initial-random-state
  118.          (congruential-rng (+ (real-time-clock) 123456789))))
  119.     (copy-random-state
  120.      (guarantee-random-state state 'MAKE-RANDOM-STATE)))))
  121.  
  122. (define (initial-random-state generate-random-seed)
  123.   ;; The numbers returned by GENERATE-RANDOM-SEED are not critical.
  124.   ;; Except for the explicitly disallowed sequences, all other
  125.   ;; sequences produce reasonable results, although some sequences
  126.   ;; might require a small number of initial generation steps to get
  127.   ;; them into the main cycle.  (See the article for details.)
  128.   (let ((seeds (flo:vector-cons r)))
  129.     (let fill ()
  130.       (do ((i 0 (fix:+ i 1)))
  131.       ((fix:= i r))
  132.     (flo:vector-set! seeds i (int:->flonum (generate-random-seed b))))
  133.       ;; Disallow cases with all seeds either 0 or b-1, since they can
  134.       ;; get locked in trivial cycles.
  135.       (if (or (let loop ((i 0))
  136.         (or (fix:= i r)
  137.             (and (flo:= (flo:vector-ref seeds i) 0.)
  138.              (loop (fix:+ i 1)))))
  139.           (let ((b-1 (flo:- b. 1.)))
  140.         (let loop ((i 0))
  141.           (or (fix:= i r)
  142.               (and (flo:= (flo:vector-ref seeds i) b-1)
  143.                (loop (fix:+ i 1)))))))
  144.       (fill)))
  145.     (%make-random-state 0 0. seeds)))
  146.  
  147. (define (congruential-rng seed)
  148.   (let ((a 16807 #|(expt 7 5)|#)
  149.     (m 2147483647 #|(- (expt 2 31) 1)|#))
  150.     (let ((m-1 (- m 1)))
  151.       (let ((seed (+ (int:remainder seed m-1) 1)))
  152.     (lambda (b)
  153.       (let ((n (int:remainder (* a seed) m)))
  154.         (set! seed n)
  155.         (int:quotient (* (- n 1) b) m-1)))))))
  156.  
  157. ;;; The RANDOM-STATE data abstraction must be built by hand because
  158. ;;; the random-number generator is needed in order to build the record
  159. ;;; abstraction.
  160.  
  161. (define-integrable (%make-random-state i b v)
  162.   (vector random-state-tag i b v))
  163.  
  164. (define (random-state? object)
  165.   (and (vector? object)
  166.        (not (fix:= (vector-length object) 0))
  167.        (eq? (vector-ref object 0) random-state-tag)))
  168.  
  169. (define-integrable random-state-tag
  170.   ((ucode-primitive string->symbol) "#[(runtime random-number)random-state]"))
  171.  
  172. (define-integrable (random-state-index s) (vector-ref s 1))
  173. (define-integrable (set-random-state-index! s x) (vector-set! s 1 x))
  174.  
  175. (define-integrable (random-state-borrow s) (vector-ref s 2))
  176. (define-integrable (set-random-state-borrow! s x) (vector-set! s 2 x))
  177.  
  178. (define-integrable (random-state-vector s) (vector-ref s 3))
  179.  
  180. (define (copy-random-state state)
  181.   (%make-random-state (random-state-index state)
  182.               (random-state-borrow state)
  183.               (flo:vector-copy (random-state-vector state))))
  184.  
  185. (define (flo:vector-copy vector)
  186.   (let ((n (flo:vector-length vector)))
  187.     (let ((result (flo:vector-cons n)))
  188.       (do ((i 0 (fix:+ i 1)))
  189.       ((fix:= i n))
  190.     (flo:vector-set! result i (flo:vector-ref vector i)))
  191.       result)))
  192.  
  193. (define (guarantee-random-state state procedure)
  194.   (if state
  195.       (begin
  196.     (if (not (random-state? state))
  197.         (error:wrong-type-argument state "random state" procedure))
  198.     state)
  199.       (let ((state *random-state*))
  200.     (if (not (random-state? state))
  201.         (error "Invalid *random-state*:" state))
  202.     state)))
  203.  
  204. (define *random-state*)
  205.  
  206. (define (initialize-package!)
  207.   (set! *random-state*
  208.     (initial-random-state
  209.      (congruential-rng (+ (real-time-clock) 123456789))))
  210.   unspecific)
  211.  
  212. (define (finalize-random-state-type!)
  213.   (add-event-receiver! event:after-restore
  214.                (lambda ()
  215.              (set! *random-state* (make-random-state #t))
  216.              unspecific))
  217.   (named-structure/set-tag-description! random-state-tag
  218.     (make-define-structure-type 'VECTOR
  219.                 'RANDOM-STATE
  220.                 '(INDEX BORROW VECTOR)
  221.                 '(1 2 3)
  222.                 (standard-unparser-method 'RANDOM-STATE #f))))