source: project/random-swb/trunk/random-swb.scm @ 6984

Last change on this file since 6984 was 6984, checked in by Ivan Raikov, 12 years ago

removed dependence on easyffi

File size: 7.6 KB
Line 
1
2;;
3;; A port of the SML/NJ implementation of a random number generator
4;; using the subtract-with-borrow (SWB) method described by Marsaglia
5;; and Zaman in _A New Class of Random Number Generators_,
6;; Ann. Applied Prob. 1(3), 1991, pp. 462-480.
7;;
8;; The SWB generator is a 30-bit generator with lags 48 and 8. It has
9;; period (2^1487 - 2^247)/105 or about 10^445. The SWB generator acts
10;; locally like a lagged Fibonacci generator and thus it is combined
11;; with a linear congruential generator (48271*a)mod(2^30-1).
12;;
13;;
14;; Copyright 2007 Ivan Raikov
15;;
16;; This program is free software: you can redistribute it and/or
17;; modify it under the terms of the GNU General Public License as
18;; published by the Free Software Foundation, either version 3 of the
19;; License, or (at your option) any later version.
20;;
21;; This program is distributed in the hope that it will be useful, but
22;; WITHOUT ANY WARRANTY; without even the implied warranty of
23;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
24;; General Public License for more details.
25;;
26;; A full copy of the GPL license can be found at
27;; <http://www.gnu.org/licenses/>."))))
28
29(require-extension srfi-4)
30
31(define-extension random-swb)
32
33(declare (export make-swb-random-state
34                 swb:random!
35                 swb:random-natural!
36                 swb:random-real!
37                 swb:random-range!
38                 swb:random
39                 swb:random-natural
40                 swb:random-real
41                 swb:random-range))
42
43
44(define nbits 30)
45(define max-word (inexact->exact (- (expt 2 nbits) 1)))
46
47(define bit29 #x20000000)
48(define lo29  #x1FFFFFFF)
49
50
51(define N 48)
52(define lag 8)
53(define offset (- N lag))
54
55
56#>
57  unsigned int cminus (unsigned int x, unsigned int y, int b)
58  {
59     if (b > 0) 
60     {
61      return x-y-1;
62      }
63     else
64     {
65      return x-y;
66      }
67  }
68<#
69
70(define cminus (foreign-lambda int "cminus" int int int))
71
72
73
74(define (swb:error x . rest)
75  (let ((port (open-output-string)))
76    (let loop ((objs (cons x rest)))
77      (if (null? objs)
78          (begin
79            (newline port)
80            (error 'random-swb (get-output-string port)))
81          (begin (display (car objs) port)
82                 (display " " port)
83                 (loop (cdr objs)))))))
84
85(define 2^-29 (expt 2 -29))
86
87;;
88;; Random number generator state representation
89;;
90;; * seeds:    seed vector
91;; * borrow:   last borrow
92;; * congseed: congruential seed
93;; * index:    index of next available value in seeds
94;;
95(define-record swb-state seeds borrow congseed index)
96
97(define-record-printer (swb-state x out)
98  (fprintf out "#(SWB seeds=~S borrow=~S congseed=0x~X index=~S)"
99           (swb-state-seeds x)
100           (swb-state-borrow x)
101           (swb-state-congseed x)
102           (swb-state-index x)))
103
104
105
106;;
107;; Linear congruential generator:
108;; multiplication by 48271 mod (2^30 - 1)
109;;
110(define  lcg-a 48271)
111(define  lcg-m max-word)
112(define  lcg-q (quotient lcg-m lcg-a))
113(define  lcg-r (modulo lcg-m lcg-a))
114
115(define (lcg s)
116  (let ((left  (fx* lcg-a (modulo s lcg-q)))
117        (right (fx* lcg-r (quotient s lcg-q))))
118    (if (fx> left right) (fx- left right)
119        (fx+ (fx- lcg-m right) left))))
120
121;;
122;; Fills seed array using subtract-with-borrow generator:
123;;
124;;  x[n] = x[n-lag] - x[n-N] - borrow
125;; Sets index to 1 and returns 0th value.
126;;
127(define (fill! state)
128 
129  (define seeds (swb-state-seeds state))
130  (define index (swb-state-index state))
131  (define borrow (swb-state-borrow state))
132
133  (define (update ix iy b)
134    (let ((x (u32vector-ref seeds ix))
135          (y (u32vector-ref seeds iy)))
136      (let ((z (cminus x y (if b 1 0))))
137        (u32vector-set! seeds iy z)
138        (if b (fx>= y x) (fx> y x)))))
139
140  (define (fillup i b)
141    (if (fx= i lag) b
142        (fillup (fx+ 1 i) (update (fx+ i offset) i b))))
143
144  (define (fillup1 i b)
145    (if (fx= i N) b
146        (fillup1 (fx+ 1 i) (update (fx- i lag) i b))))
147
148  (swb-state-borrow-set! state (fillup1 lag (fillup 0 borrow)))
149  (swb-state-index-set! state 1)
150  (u32vector-ref seeds 0))
151
152;;
153;; Creates an initial seed vector and generator state.  Fills the seed
154;; array one bit at a time by taking the leading bit of the xor of a
155;; shift register and a congruential sequence.  The congruential
156;; generator is (c*48271) mod (2^30 - 1).  The shift register
157;; generator is c(I + L18)(I + R13).  The same congruential generator
158;; continues to be used as a mixing generator with the SWB generator.
159;;
160(define (make-swb-random-state congy shrgx)
161
162  (define (mki v)
163    (match v ((i c s)
164              (let* ((c1 (lcg  c)))
165                (let* ((s1 (bitwise-xor s  (fxshl s  18)))
166                       (s2 (bitwise-xor s1 (fxshr s1 13)))
167                       (i1 (fx+ (bitwise-and lo29 (fxshr i 1))
168                                (bitwise-and bit29 (bitwise-xor c1 s2)))))
169                  (list i1 c1 s2))))
170              (else (swb:error 'mki "invalid argument " v))))
171
172  (define (iterate n v)
173    (if (zero? n) v
174        (iterate (fx- n 1) (mki v))))
175
176  (define (mkseed congx shrgx)
177    (iterate nbits (list 0 congx shrgx)))
178
179  (define (genseed n seeds congx shrgx)
180    (if (zero? n) (values seeds congx)
181        (let ((v (mkseed congx shrgx)))
182          (match v  ((seed congx1 shrgx1)
183                     (genseed (fx- n 1) (cons seed seeds) congx1  shrgx1))
184                 (else (swb:error 'genseed "invalid mkseed return value " v))))))
185
186  (define congx (fx+ 1 (fxshl (bitwise-and congy max-word) 1)))
187
188  (let-values (((seeds congx)  (genseed N (list) congx shrgx)))
189    (make-swb-state (list->u32vector seeds) #f congx 0)))
190
191;;
192;; Computes the next random number. The tweak function xor's the
193;; number from the SWB generator with a number from the linear
194;; congruential generator.
195;;
196(define (swb:random! s)
197  (match s 
198         (($ swb-state seeds borrow congseed index)
199          (let ((tweak (lambda (i) 
200                         (let ((c (lcg congseed)))
201                           (swb-state-congseed-set! s c)
202                           (bitwise-xor i c)))))
203            (if (fx= index N)
204                (tweak (fill! s))
205                (begin
206                  (swb-state-index-set! s (fx+ 1 index))
207                  (tweak (u32vector-ref seeds index))))))
208         (else (swb:error 'swb:random "invalid state " s))))
209
210(define (swb:random-natural! s)
211  (bitwise-and (swb:random! s) lo29))
212
213(define (swb:random-real! s)
214  (let* ((n1 (swb:random-natural! s))
215         (n2 (swb:random-natural! s)))
216    (* 2^-29 (+ n1 (* n2 2^-29)))))
217
218(define (swb:random-range! i j)
219  (if (fx< j i)
220      (swb:error 'swb:random-range "argument j < i")
221      (if (and (fixnum? i) (fixnum? j))
222          (let ((R  (* 2^-29 (fx+ 1 (fx- j i)))))
223            (lambda (s) (fx+ i (inexact->exact (truncate (* R (swb:random-natural! s)))))))
224          (let ((R  (fx+ 1 (fx- j i))))
225            (lambda (s) (+ i (* R (swb:random-real! s))))))))
226         
227         
228;;
229;; Side-effect-free variants of the above functions.
230;;
231(define (swb:random s)
232  (match s 
233         (($ swb-state seeds borrow congseed index)
234          (let ((tweak (lambda (s i) 
235                         (let ((c (lcg congseed)))
236                           (swb-state-congseed-set! s c)
237                           (bitwise-xor i c)))))
238            (if (fx= index N)
239                (let* ((s1  (make-swb-state seeds borrow congseed index))
240                       (i   (fill! s1))
241                       (v   (tweak s1 i)))
242                  (values s1 v))
243                (let* ((s1 (make-swb-state seeds borrow congseed (fx+ 1 index)))
244                       (v  (tweak s1 (u32vector-ref seeds index))))
245                  (values s1 v)))))
246         (else (swb:error 'swb:random "invalid state " s))))
247
248(define (swb:random-natural s)
249  (let-values (((s v) (swb:random s)))
250    (values s (bitwise-and v lo29))))
251
252(define (swb:random-real s)
253  (let-values (((s1 n1) (swb:random-natural s)))
254    (let-values (((s2 n2) (swb:random-natural s1)))
255       (values s2 (* 2^-29 (+ n1 (* n2 2^-29)))))))
256
257(define (swb:random-range i j)
258  (if (fx< j i)
259      (swb:error 'swb:random-range "argument j < i")
260      (if (and (fixnum? i) (fixnum? j))
261          (let ((R  (* 2^-29 (fx+ 1 (fx- j i)))))
262            (lambda (s) 
263              (let-values (((s1 v) (swb:random-natural s)))
264                (values s1 (fx+ i (inexact->exact (truncate (* R v))))))))
265          (let ((R  (fx+ 1 (fx- j i))))
266            (lambda (s) 
267              (let-values (((s1 v)  (swb:random-real s)))
268                  (values s1 (+ i (* R v)))))))))
269         
270         
Note: See TracBrowser for help on using the repository browser.