summaryrefslogtreecommitdiffstats
path: root/randinex.scm
diff options
context:
space:
mode:
Diffstat (limited to 'randinex.scm')
-rw-r--r--randinex.scm99
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)