1 files changed, 99 insertions, 0 deletions
diff --git a/randinex.scm b/randinex.scm
new file mode 100644
index 0000000..1c2b702
--- /dev/null
+++ b/randinex.scm
@@ -0,0 +1,99 @@
+;;;"randinex.scm" Pseudo-Random inexact real numbers for scheme.
+;;; Copyright (C) 1991, 1993 Aubrey Jaffer.
+;
+;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.  I have made no warrantee or representation that the operation of
+;this software will be error-free, and I am under no obligation to
+;provide any services, by way of maintenance, update, or otherwise.
+;
+;3.  In conjunction with products arising from the use of this
+;material, there shall be no use of my name in any advertising,
+;promotional, or sales literature without prior written consent in
+;each case.
+
+;This file is loaded by random.scm if inexact numbers are supported by
+;the implementation.
+
+;;; Fixed sphere and normal functions from: Harald Hanche-Olsen
+
+(define random:float-radix
+  (+ 1 (exact->inexact random:MASK)))
+
+;;; This determines how many chunks will be neccessary to completely
+;;; fill up an inexact real.
+(define (random:size-float l x)
+  (cond ((= 1.0 (+ 1 x)) l)
+	((= 4 l) l)
+	(else (random:size-float (+ l 1) (/ x random:float-radix)))))
+(define random:chunks/float (random:size-float 0 1.0))
+
+(define (random:uniform-chunk n state)
+  (if (= 1 n)
+      (/ (exact->inexact (random:chunk state))
+	 random:float-radix)
+      (/ (+ (random:uniform-chunk (- n 1) state)
+	    (exact->inexact (random:chunk state)))
+	 random:float-radix)))
+
+;;; Generate an inexact real between 0 and 1.
+(define (random:uniform state)
+  (random:uniform-chunk random:chunks/float state))
+
+;;; If x and y are independent standard normal variables, then with
+;;; x=r*cos(t), y=r*sin(t), we find that t is uniformly distributed
+;;; over [0,2*pi] and the cumulative distribution of r is
+;;; 1-exp(-r^2/2).  This latter means that u=exp(-r^2/2) is uniformly
+;;; distributed on [0,1], so r=sqrt(-2 log u) can be used to generate r.
+
+(define (random:normal-vector! vect . args)
+  (let ((state (if (null? args) *random-state* (car args)))
+	(sum2 0))
+    (let ((do! (lambda (k x)
+		 (vector-set! vect k x)
+		 (set! sum2 (+ sum2 (* x x))))))
+      (do ((n (- (vector-length vect) 1) (- n 2)))
+	  ((negative? n) sum2)
+	(let ((t (* 6.28318530717958 (random:uniform state)))
+	      (r (sqrt (* -2 (log (random:uniform state))))))
+	  (do! n (* r (cos t)))
+	  (if (positive? n) (do! (- n 1) (* r (sin t)))))))))
+
+(define random:normal
+  (let ((vect (make-vector 1)))
+    (lambda args 
+      (apply random:normal-vector! vect args)
+      (vector-ref vect 0))))
+
+;;; For the uniform distibution on the hollow sphere, pick a normal
+;;; family and scale.
+
+(define (random:hollow-sphere! vect . args)
+  (let ((ms (sqrt (apply random:normal-vector! vect args))))
+    (do ((n (- (vector-length vect) 1) (- n 1)))
+	((negative? n))
+      (vector-set! vect n (/ (vector-ref vect n) ms)))))
+
+;;; For the uniform distribution on the solid sphere, note that in
+;;; this distribution the length r of the vector has cumulative
+;;; distribution r^n; i.e., u=r^n is uniform [0,1], so r kan be
+;;; generated as r=u^(1/n).
+
+(define (random:solid-sphere! vect . args)
+  (apply random:hollow-sphere! vect args)
+  (let ((r (expt (random:uniform (if (null? args) *random-state* (car args)))
+		 (/ (vector-length vect)))))
+    (do ((n (- (vector-length vect) 1) (- n 1)))
+	((negative? n))
+      (vector-set! vect n (* r (vector-ref vect n))))))
+
+(define (random:exp . args)
+  (let ((state (if (null? args) *random-state* (car args))))
+    (- (log (random:uniform state)))))
+
+(require 'random)