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
|
;;;"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
;;; 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))))
;;; 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 can 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)))))
|
