diff options
Diffstat (limited to 'randinex.scm')
| -rw-r--r-- | randinex.scm | 96 |
1 files changed, 63 insertions, 33 deletions
diff --git a/randinex.scm b/randinex.scm index e6dc48b..8a0afd1 100644 --- a/randinex.scm +++ b/randinex.scm @@ -1,5 +1,5 @@ ;;;"randinex.scm" Pseudo-Random inexact real numbers for scheme. -;;; Copyright (C) 1991, 1993 Aubrey Jaffer. +;;; Copyright (C) 1991, 1993, 1999 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 @@ -20,28 +20,49 @@ ;This file is loaded by random.scm if inexact numbers are supported by ;the implementation. -;;; Fixed sphere and normal functions from: Harald Hanche-Olsen +;;; Sphere and normal functions corrections from: Harald Hanche-Olsen ;;; Generate an inexact real between 0 and 1. -(define random:uniform - (letrec ((random:chunks/float ; how many chunks fill an inexact? - (letrec ((random:size-float - (lambda (l x) - (cond ((= 1.0 (+ 1 x)) l) - ((= 4 l) l) - (else (random:size-float (+ l 1) (/ x 256.0))))))) - (random:size-float 0 1.0))) - - (random:uniform-chunk - (lambda (n state) - (if (= 1 n) - (/ (exact->inexact (random:chunk state)) - 256.0) - (/ (+ (random:uniform-chunk (- n 1) state) - (exact->inexact (random:chunk state))) - 256.0))))) - (lambda (state) - (random:uniform-chunk random:chunks/float state)))) +(define random:uniform1 + ; how many chunks fill an inexact? + (do ((random:chunks/float 0 (+ 1 random:chunks/float)) + (smidgen 1.0 (/ smidgen 256.0))) + ((or (= 1.0 (+ 1 smidgen)) (= 4 random:chunks/float)) + (lambda (state) + (do ((cnt random:chunks/float (+ -1 cnt)) + (uni (/ (random:chunk state) 256.0) + (/ (+ uni (random:chunk state)) 256.0))) + ((= 1 cnt) uni)))))) + + +;;@args +;;@args state +;;Returns an uniformly distributed inexact real random number in the +;;range between 0 and 1. +(define (random:uniform . args) + (random:uniform1 (if (null? args) *random-state* (car args)))) + + +;;@args +;;@args state +;;Returns an inexact real in an exponential distribution with mean 1. For +;;an exponential distribution with mean @var{u} use +;;@w{@code{(* @var{u} (random:exp))}}. +(define (random:exp . args) + (- (log (random:uniform1 (if (null? args) *random-state* (car args)))))) + + +;;@args +;;@args state +;;Returns an inexact real in a normal distribution with mean 0 and +;;standard deviation 1. For a normal distribution with mean @var{m} and +;;standard deviation @var{d} use +;;@w{@code{(+ @var{m} (* @var{d} (random:normal)))}}. +(define (random:normal . args) + (let ((vect (make-vector 1))) + (apply random:normal-vector! vect args) + (vector-ref vect 0))) + ;;; 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 @@ -49,6 +70,10 @@ ;;; 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. +;;@args vect +;;@args vect state +;;Fills @1 with inexact real random numbers which are independent +;;and standard normally distributed (i.e., with mean 0 and variance 1). (define (random:normal-vector! vect . args) (let ((state (if (null? args) *random-state* (car args))) (sum2 0)) @@ -57,39 +82,44 @@ (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)))))) + (let ((t (* 6.28318530717958 (random:uniform1 state))) + (r (sqrt (* -2 (log (random:uniform1 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. +;;@args vect +;;@args vect state +;;Fills @1 with inexact real random numbers the sum of whose +;;squares is less than 1.0. Thinking of @1 as coordinates in +;;space of dimension @var{n} = @code{(vector-length @1)}, the +;;coordinates are uniformly distributed within the unit @var{n}-shere. +;;The sum of the squares of the numbers is returned. (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 can be ;;; generated as r=u^(1/n). +;;@args vect +;;@args vect state +;;Fills @1 with inexact real random numbers the sum of whose +;;squares is equal to 1.0. Thinking of @1 as coordinates in space +;;of dimension n = @code{(vector-length @1)}, the coordinates are +;;uniformly distributed over the surface of the unit n-shere. (define (random:solid-sphere! vect . args) (apply random:hollow-sphere! vect args) - (let ((r (expt (random:uniform (if (null? args) *random-state* (car args))) + (let ((r (expt (random:uniform1 (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))))) |