1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
|
;;;"randinex.scm" Pseudo-Random inexact real numbers for scheme.
;;; Copyright (C) 1991, 1993, 1999 Aubrey Jaffer
;
;Permission to copy this software, to modify it, to redistribute it,
;to distribute modified versions, 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.
;;; Sphere and normal functions corrections from: Harald Hanche-Olsen
;;; Generate an inexact real between 0 and 1.
(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
;;; 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 *2pi (* 8 (atan 1)))
;;@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))
(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 (* *2pi (random:uniform1 state)))
(r (sqrt (* -2 (log (random:uniform1 state))))))
(do! n (* r (cos t)))
(if (positive? n) (do! (- n 1) (* r (sin t)))))))))
;;; 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 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: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 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:solid-sphere! vect . args)
(apply random:hollow-sphere! vect 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) r)
(vector-set! vect n (* r (vector-ref vect n))))))
|
