summaryrefslogtreecommitdiffstats
path: root/randinex.scm
diff options
context:
space:
mode:
authorBryan Newbold <bnewbold@robocracy.org>2017-02-20 00:05:28 -0800
committerBryan Newbold <bnewbold@robocracy.org>2017-02-20 00:05:28 -0800
commitbd9733926076885e3417b74de76e4c9c7bc56254 (patch)
tree2c99dced547d48407ad44cb0e45e31bb4d02ce43 /randinex.scm
parentfa3f23105ddcf07c5900de47f19af43d1db1b597 (diff)
downloadslib-bd9733926076885e3417b74de76e4c9c7bc56254.tar.gz
slib-bd9733926076885e3417b74de76e4c9c7bc56254.zip
Import Upstream version 2c7upstream/2c7
Diffstat (limited to 'randinex.scm')
-rw-r--r--randinex.scm96
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)))))