source: project/release/4/random-swb/trunk/random-swb.scm @ 14477

Last change on this file since 14477 was 14477, checked in by Ivan Raikov, 10 years ago

random-swb ported to Chicken 4

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-2009 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
30(module random-swb
31
32
33 (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  (import scheme chicken extras srfi-4 matchable)
44
45  (require-extension srfi-4 matchable)
46
47
48(define nbits 30)
49(define max-word (inexact->exact (- (expt 2 nbits) 1)))
50
51(define bit29 #x20000000)
52(define lo29  #x1FFFFFFF)
53
54
55(define N 48)
56(define lag 8)
57(define offset (- N lag))
58
59
60#>
61  unsigned int cminus (unsigned int x, unsigned int y, int b)
62  {
63     if (b > 0) 
64     {
65      return x-y-1;
66      }
67     else
68     {
69      return x-y;
70      }
71  }
72<#
73
74(define cminus (foreign-lambda int "cminus" int int int))
75
76
77
78(define (swb:error x . rest)
79  (let ((port (open-output-string)))
80    (let loop ((objs (cons x rest)))
81      (if (null? objs)
82          (begin
83            (newline port)
84            (error 'random-swb (get-output-string port)))
85          (begin (display (car objs) port)
86                 (display " " port)
87                 (loop (cdr objs)))))))
88
89(define 2^-29 (expt 2 -29))
90
91;;
92;; Random number generator state representation
93;;
94;; * seeds:    seed vector
95;; * borrow:   last borrow
96;; * congseed: congruential seed
97;; * index:    index of next available value in seeds
98;;
99(define-record swb-state seeds borrow congseed index)
100
101(define-record-printer (swb-state x out)
102  (fprintf out "#(SWB seeds=~S borrow=~S congseed=0x~X index=~S)"
103           (swb-state-seeds x)
104           (swb-state-borrow x)
105           (swb-state-congseed x)
106           (swb-state-index x)))
107
108
109
110;;
111;; Linear congruential generator:
112;; multiplication by 48271 mod (2^30 - 1)
113;;
114(define  lcg-a 48271)
115(define  lcg-m max-word)
116(define  lcg-q (quotient lcg-m lcg-a))
117(define  lcg-r (modulo lcg-m lcg-a))
118
119(define (lcg s)
120  (let ((left  (fx* lcg-a (modulo s lcg-q)))
121        (right (fx* lcg-r (quotient s lcg-q))))
122    (if (fx> left right) (fx- left right)
123        (fx+ (fx- lcg-m right) left))))
124
125;;
126;; Fills seed array using subtract-with-borrow generator:
127;;
128;;  x[n] = x[n-lag] - x[n-N] - borrow
129;; Sets index to 1 and returns 0th value.
130;;
131(define (fill! state)
132 
133  (define seeds (swb-state-seeds state))
134  (define index (swb-state-index state))
135  (define borrow (swb-state-borrow state))
136
137  (define (update ix iy b)
138    (let ((x (u32vector-ref seeds ix))
139          (y (u32vector-ref seeds iy)))
140      (let ((z (cminus x y (if b 1 0))))
141        (u32vector-set! seeds iy z)
142        (if b (fx>= y x) (fx> y x)))))
143
144  (define (fillup i b)
145    (if (fx= i lag) b
146        (fillup (fx+ 1 i) (update (fx+ i offset) i b))))
147
148  (define (fillup1 i b)
149    (if (fx= i N) b
150        (fillup1 (fx+ 1 i) (update (fx- i lag) i b))))
151
152  (swb-state-borrow-set! state (fillup1 lag (fillup 0 borrow)))
153  (swb-state-index-set! state 1)
154  (u32vector-ref seeds 0))
155
156;;
157;; Creates an initial seed vector and generator state.  Fills the seed
158;; array one bit at a time by taking the leading bit of the xor of a
159;; shift register and a congruential sequence.  The congruential
160;; generator is (c*48271) mod (2^30 - 1).  The shift register
161;; generator is c(I + L18)(I + R13).  The same congruential generator
162;; continues to be used as a mixing generator with the SWB generator.
163;;
164(define (make-swb-random-state congy shrgx)
165
166  (define (mki v)
167    (match v ((i c s)
168              (let* ((c1 (lcg  c)))
169                (let* ((s1 (bitwise-xor s  (fxshl s  18)))
170                       (s2 (bitwise-xor s1 (fxshr s1 13)))
171                       (i1 (fx+ (bitwise-and lo29 (fxshr i 1))
172                                (bitwise-and bit29 (bitwise-xor c1 s2)))))
173                  (list i1 c1 s2))))
174              (else (swb:error 'mki "invalid argument " v))))
175
176  (define (iterate n v)
177    (if (zero? n) v
178        (iterate (fx- n 1) (mki v))))
179
180  (define (mkseed congx shrgx)
181    (iterate nbits (list 0 congx shrgx)))
182
183  (define (genseed n seeds congx shrgx)
184    (if (zero? n) (values seeds congx)
185        (let ((v (mkseed congx shrgx)))
186          (match v  ((seed congx1 shrgx1)
187                     (genseed (fx- n 1) (cons seed seeds) congx1  shrgx1))
188                 (else (swb:error 'genseed "invalid mkseed return value " v))))))
189
190  (define congx (fx+ 1 (fxshl (bitwise-and congy max-word) 1)))
191
192  (let-values (((seeds congx)  (genseed N (list) congx shrgx)))
193    (make-swb-state (list->u32vector seeds) #f congx 0)))
194
195;;
196;; Computes the next random number. The tweak function xor's the
197;; number from the SWB generator with a number from the linear
198;; congruential generator.
199;;
200(define (swb:random! s)
201  (match s 
202         (($ swb-state seeds borrow congseed index)
203          (let ((tweak (lambda (i) 
204                         (let ((c (lcg congseed)))
205                           (swb-state-congseed-set! s c)
206                           (bitwise-xor i c)))))
207            (if (fx= index N)
208                (tweak (fill! s))
209                (begin
210                  (swb-state-index-set! s (fx+ 1 index))
211                  (tweak (u32vector-ref seeds index))))))
212         (else (swb:error 'swb:random "invalid state " s))))
213
214(define (swb:random-natural! s)
215  (bitwise-and (swb:random! s) lo29))
216
217(define (swb:random-real! s)
218  (let* ((n1 (swb:random-natural! s))
219         (n2 (swb:random-natural! s)))
220    (* 2^-29 (+ n1 (* n2 2^-29)))))
221
222(define (swb:random-range! i j)
223  (if (fx< j i)
224      (swb:error 'swb:random-range "argument j < i")
225      (if (and (fixnum? i) (fixnum? j))
226          (let ((R  (* 2^-29 (fx+ 1 (fx- j i)))))
227            (lambda (s) (fx+ i (inexact->exact (truncate (* R (swb:random-natural! s)))))))
228          (let ((R  (fx+ 1 (fx- j i))))
229            (lambda (s) (+ i (* R (swb:random-real! s))))))))
230         
231         
232;;
233;; Side-effect-free variants of the above functions.
234;;
235(define (swb:random s)
236  (match s 
237         (($ swb-state seeds borrow congseed index)
238          (let ((tweak (lambda (s i) 
239                         (let ((c (lcg congseed)))
240                           (swb-state-congseed-set! s c)
241                           (bitwise-xor i c)))))
242            (if (fx= index N)
243                (let* ((s1  (make-swb-state seeds borrow congseed index))
244                       (i   (fill! s1))
245                       (v   (tweak s1 i)))
246                  (values s1 v))
247                (let* ((s1 (make-swb-state seeds borrow congseed (fx+ 1 index)))
248                       (v  (tweak s1 (u32vector-ref seeds index))))
249                  (values s1 v)))))
250         (else (swb:error 'swb:random "invalid state " s))))
251
252(define (swb:random-natural s)
253  (let-values (((s v) (swb:random s)))
254    (values s (bitwise-and v lo29))))
255
256(define (swb:random-real s)
257  (let-values (((s1 n1) (swb:random-natural s)))
258    (let-values (((s2 n2) (swb:random-natural s1)))
259       (values s2 (* 2^-29 (+ n1 (* n2 2^-29)))))))
260
261(define (swb:random-range i j)
262  (if (fx< j i)
263      (swb:error 'swb:random-range "argument j < i")
264      (if (and (fixnum? i) (fixnum? j))
265          (let ((R  (* 2^-29 (fx+ 1 (fx- j i)))))
266            (lambda (s) 
267              (let-values (((s1 v) (swb:random-natural s)))
268                (values s1 (fx+ i (inexact->exact (truncate (* R v))))))))
269          (let ((R  (fx+ 1 (fx- j i))))
270            (lambda (s) 
271              (let-values (((s1 v)  (swb:random-real s)))
272                  (values s1 (+ i (* R v)))))))))
273         
274         
275)
Note: See TracBrowser for help on using the repository browser.