From: Chris Hanson <org/chris-hanson/cph>
Date: Tue, 9 Feb 1993 00:25:45 +0000 (+0000)
Subject: Completely new generator, interface is upwards compatible with
X-Git-Tag: 20090517-FFI~8547
X-Git-Url: https://birchwood-abbey.net/git?a=commitdiff_plain;h=9fc0ef5dbf38fb83f1db2ec1ca5b88d9acb793a2;p=mit-scheme.git

Completely new generator, interface is upwards compatible with
previous design.  RANDOM accepts any real number as a modulus
argument.  MAKE-RANDOM-STATE additionally accepts an exact integer,
which it uses as a seed to create the state, thus providing a simple
means of producing a specific random state.
---

diff --git a/v7/src/runtime/random.scm b/v7/src/runtime/random.scm
index ccabe1bb2..b9144b38c 100644
--- a/v7/src/runtime/random.scm
+++ b/v7/src/runtime/random.scm
@@ -1,112 +1,170 @@
-;;;; Pseudo-Random number generator for scheme.
-;;; Copyright (C) 1991 Aubrey Jaffer.
-;;; From slib1c4, modified for MIT Scheme by cph.
-;;; $Id: random.scm,v 14.4 1993/01/13 08:48:34 cph Exp $
-;;;
-;;;   (random n)					procedure
-;;;   (random n state)					procedure
-;;; 
-;;; Accepts a positive integer or real n and returns a number of the
-;;; same type between zero (inclusive) and n (exclusive).  The values
-;;; returned have a uniform distribution.
-;;;
-;;; The optional argument state must be of the type produced by
-;;; (make-random-state).  It defaults to the value of the variable
-;;; *random-state*.  This object is used to maintain the state of the
-;;; pseudo-random-number generator and is altered as a side effect of the
-;;; RANDOM operation.
-;;;
-;;;   *random-state*					variable
-;;;
-;;; Holds a data structure that encodes the internal state of the
-;;; random-number generator that RANDOM uses by default.  The nature of
-;;; this data structure is implementation-dependent.  It may be printed
-;;; out and successfully read back in, but may or may not function
-;;; correctly as a random-number state object in another implementation.
-;;;
-;;;   (make-random-state)				procedure
-;;;   (make-random-state state)				procedure
-;;;
-;;; Returns a new object of type suitable for use as the value of the
-;;; variable *random-state* and as second argument to RANDOM.  If argument
-;;; state is given, a copy of it is returned.  Otherwise a copy of
-;;; *random-state* is returned.
-;;;------------------------------------------------------------------
+#| -*-Scheme-*-
+
+$Id: random.scm,v 14.5 1993/02/09 00:25:45 cph Exp $
+
+Copyright (c) 1993 Massachusetts Institute of Technology
+
+This material was developed by the Scheme project at the Massachusetts
+Institute of Technology, Department of Electrical Engineering and
+Computer Science.  Permission to copy this software, to redistribute
+it, and to use it for any purpose is granted, subject to the following
+restrictions and understandings.
+
+1. Any copy made of this software must include this copyright notice
+in full.
+
+2. Users of this software agree to make their best efforts (a) to
+return to the MIT Scheme project any improvements or extensions that
+they make, so that these may be included in future releases; and (b)
+to inform MIT of noteworthy uses of this software.
+
+3. All materials developed as a consequence of the use of this
+software shall duly acknowledge such use, in accordance with the usual
+standards of acknowledging credit in academic research.
+
+4. MIT has made no warrantee or representation that the operation of
+this software will be error-free, and MIT is under no obligation to
+provide any services, by way of maintenance, update, or otherwise.
+
+5. In conjunction with products arising from the use of this material,
+there shall be no use of the name of the Massachusetts Institute of
+Technology nor of any adaptation thereof in any advertising,
+promotional, or sales literature without prior written consent from
+MIT in each case. |#
+
+;;;; Random Number Generator
+;;; package: (runtime random-number)
+
+(declare (usual-integrations))
+
+(define-integrable (int:->flonum n)
+  ((ucode-primitive integer->flonum 2) n #b10))
 
-(define random
-  (let ((state-tap-1 24)
-	(state-size 55)
-	(chunk-size 24)
-	(chunk-sup #x1000000))
-
-    (define (get-bits n state)
-      (let loop ((n n))
-	(let ((p (vector-ref state state-size)))
-	  (let ((i (fix:modulo (fix:- p state-tap-1) state-size))
-		(chunk (vector-ref state p)))
-	    (vector-set! state p (fix:xor (vector-ref state i) chunk))
-	    (vector-set! state state-size (fix:modulo (fix:- p 1) state-size))
-	    (cond ((fix:= n chunk-size)
-		   chunk)
-		  ((fix:< n chunk-size)
-		   (fix:and chunk (fix:- (fix:lsh 1 n) 1)))
-		  (else
-		   (+ chunk (* chunk-sup (loop (fix:- n chunk-size))))))))))
-
-    (define (fix:modulo n d)
-      ;; Specialized for nonnegative D.
-      (let ((r (fix:remainder n d)))
-	(if (or (fix:= r 0)
-		(not (fix:< n 0)))
-	    r
-	    (fix:+ r d))))
-
-    (lambda (modulus #!optional state)
-      (if (not (and (exact-integer? modulus) (> modulus 0)))
-	  (error:wrong-type-argument modulus "exact positive integer" 'RANDOM))
-      (let ((state (if (default-object? state) *random-state* state)))
-	(if (not (random-state? state))
-	    (error:wrong-type-argument state
-				       "random number state"
-				       'MAKE-RANDOM-STATE))
-	(let ((ilen (exact-nonnegative-integer-length modulus))
-	      (state (random-state-bits state)))
-	  (do ((r (get-bits ilen state)
-		  (get-bits ilen state))) ;this could be improved.
-	      ((< r modulus) r)))))))
+;;; A "subtract-with-carry" RNG, based on the algorithm from "A New
+;;; Class of Random Number Generators", George Marsaglia and Arif
+;;; Zaman, The Annals of Applied Probability, Vol. 1, No. 3, 1991.
 
-(define-structure (random-state (constructor %make-random-state))
-  (bits false read-only true))
+;;; The basic algorithm produces a sequence of values x[n] by
+;;;      let t = (x[n-s] - x[n-r] - borrow[n-1])
+;;;      if (t >= 0)
+;;;        then x[n] = t, borrow[n] = 0
+;;;        else x[n] = t + b, borrow[n] = 1
+;;; where the constants R, S, and B are chosen according to some
+;;; constraints that are explained in the article.  Finding
+;;; appropriate constants is compute-intensive; the constants used
+;;; here are taken from the article and are claimed to represent
+;;; "hundreds of hours" of compute time.  The period of this generator
+;;; is (- (EXPT B R) (EXPT B S)), which is approximately (EXPT 10 414).
 
+(define-integrable r 43)
+(define-integrable s 22)
+(define-integrable b 4294967291 #|(- (expt 2 32) 5)|#)
+(define-integrable b. 4294967291. #|(exact->inexact b)|#)
+
+(define (random modulus #!optional state)
+  (if (not (and (real? modulus) (< 0 modulus)))
+      (error:wrong-type-argument modulus "positive real" 'RANDOM))
+  (let ((element
+	 (get-next-element
+	  (guarantee-random-state (if (default-object? state) #f state)
+				  'RANDOM))))
+    ;; Kludge: an exact integer modulus means that result is an exact
+    ;; integer.  Otherwise, the result is a real number.
+    (cond ((flo:flonum? modulus)
+	   (flo:* element modulus))
+	  ((int:integer? modulus)
+	   (flo:truncate->exact (flo:* element (int:->flonum modulus))))
+	  (else
+	   (* (inexact->exact element) modulus)))))
+
+(define (get-next-element state)
+  (let ((mask (set-interrupt-enables! interrupt-mask/gc-ok)))
+    (let ((index (random-state-index state))
+	  (vector (random-state-vector state)))
+      (let ((element (vector-ref vector index)))
+	(let ((difference
+	       (flo:- (flo:- (vector-ref vector
+					 (if (fix:< (fix:- index s) 0)
+					     (fix:+ (fix:- index s) r)
+					     (fix:- index s)))
+			     element)
+		      (random-state-borrow state))))
+	  (if (flo:< difference 0.)
+	      (begin
+		(vector-set! vector index (flo:+ difference b.))
+		(set-random-state-borrow! state 1.))
+	      (begin
+		(vector-set! vector index difference)
+		(set-random-state-borrow! state 0.)))
+	  (set-random-state-index! state
+				   (if (fix:= (fix:+ index 1) r)
+				       0
+				       (fix:+ index 1))))
+	(set-interrupt-enables! mask)
+	(flo:/ element b.)))))
+
 (define (make-random-state #!optional state)
-  (let ((state (if (default-object? state) *random-state* state)))
-    (if (not (random-state? state))
-	(error:wrong-type-argument state
-				   "random number state"
-				   'MAKE-RANDOM-STATE))
-    (%make-random-state (vector-copy (random-state-bits state)))))
-
-(define *random-state*
-  (%make-random-state
-   (vector #xd909ef #xfd330a #xe33f78 #x76783f #xf3675f
-	   #xb54ef8 #x0be455 #xa67946 #x0bcd56 #xfabcde
-	   #x9cbd3e #x3fd3ef #xe064ef #xdddecc #x344442
-	   #x854444 #x4c5192 #xc03662 #x547345 #x70abcd
-	   #x1bbdac #x616c5a #xa982ef #x105996 #x5f0ccc
-	   #x1ea055 #xfe2acd #x1891c1 #xe66902 #x6912bc
-	   #x2678e1 #x612222 #x907abc #x4ad682 #x9cdd14
-	   #x577988 #x5b8924 #x871c9c #xd1e67b #x8b0a32
-	   #x578ef2 #x28274e #x823ef5 #x845678 #xe67890
-	   #x5890ab #x851fa9 #x13efa1 #xb12278 #xdaf805
-	   #xa0befc #x0068a7 #xe024fd #xa7b690 #x27f357
-	   0)))
-
-(define exact-nonnegative-integer-length
-  (let ((powers-of-two
-	 (let loop ((n 1))
-	   (cons n (delay (loop (* 2 n)))))))
-    (lambda (n)
-      (let loop ((powers-of-two powers-of-two) (e 0))
-	(if (< n (car powers-of-two))
-	    e
-	    (loop (force (cdr powers-of-two)) (fix:+ e 1)))))))
\ No newline at end of file
+  (let ((state (if (default-object? state) #f state)))
+    (if (exact-integer? state)
+	(initial-random-state (congruential-rng (+ state 123456789)))
+	(let ((state (guarantee-random-state state 'MAKE-RANDOM-STATE)))
+	  (%make-random-state (random-state-index state)
+			      (random-state-borrow state)
+			      (vector-copy (random-state-vector state)))))))
+
+(define (initial-random-state generate-random-seed)
+  ;; The numbers returned by GENERATE-RANDOM-SEED are not critical.
+  ;; Except for the explicitly disallowed sequences, all other
+  ;; sequences produce reasonable results, although some sequences
+  ;; might require a small number of initial generation steps to get
+  ;; them into the main cycle.  (See the article for details.)
+  (let ((seeds (make-vector r)))
+    (let fill ()
+      (do ((i 0 (fix:+ i 1)))
+	  ((fix:= i r))
+	(vector-set! seeds i (exact->inexact (generate-random-seed b))))
+      ;; Disallow cases with all seeds either 0 or b-1, since they can
+      ;; get locked in trivial cycles.
+      (if (or (let loop ((i 0))
+		(or (fix:= i r)
+		    (and (flo:= (vector-ref seeds i) 0)
+			 (loop (fix:+ i 1)))))
+	      (let ((b-1 (flo:- b. 1.)))
+		(let loop ((i 0))
+		  (or (fix:= i r)
+		      (and (flo:= (vector-ref seeds i) b-1)
+			   (loop (fix:+ i 1)))))))
+	  (fill)))
+    (%make-random-state 0 0. seeds)))
+
+(define (congruential-rng seed)
+  (let ((a 16807 #|(expt 7 5)|#)
+	(m 2147483647 #|(- (expt 2 31) 1)|#))
+    (let ((m-1 (- m 1)))
+      (let ((seed (+ (modulo seed m-1) 1)))
+	(lambda (b)
+	  (let ((n (modulo (* a seed) m)))
+	    (set! seed n)
+	    (quotient (* (- n 1) b) m-1)))))))
+
+(define-structure (random-state (constructor %make-random-state))
+  index
+  borrow
+  vector)
+
+(define (guarantee-random-state state procedure)
+  (if state
+      (begin
+	(if (not (random-state? state))
+	    (error:wrong-type-argument state "random state" procedure))
+	state)
+      (let ((state *random-state*))
+	(if (not (random-state? state))
+	    (error "Invalid *random-state*:" state))
+	state)))
+
+(define *random-state*)
+
+(define (initialize-package!)
+  (set! *random-state* (make-random-state 0))
+  unspecific)
\ No newline at end of file