[plt-scheme] Re: Generating Random bignums

From: Jean-Paul.ROY at unice.fr (Jean-Paul.ROY at unice.fr)
Date: Mon Jan 8 05:54:50 EST 2007

I use for several years the following code (due to Dorai Sitaram I  
think) to generate random bignums and floats :

(module srandom mzscheme
   (provide srandom)

   (define (current-time)
     (expt (current-milliseconds) 2))

   (define srandom:initialize-state #f)
   (define srandom:seed-state #f)
   (define srandom #f)

   (let ((ij #f)
         (kl #f)
         (u #f)
         (c #f)
         (cd #f)
         (cm #f)
         (i97 #f)
         (j97 #f))

     (set! srandom:initialize-state
           (lambda ()
             (set! ij 1802)
             (set! kl 9373)
             (set! u (make-vector 97 0))
             (set! c (/ 362436.0 16777216.0))
             (set! cd (/ 7654321.0 16777216.0))
             (set! cm (/ 16777213.0 16777216.0))
             (set! i97 96)
             (set! j97 32)))

     (srandom:initialize-state)

     (set! srandom:seed-state
           (lambda z
             (let ((z-len (length z))
                   (tt (current-time)))
               (set! ij
                     (if (>= z-len 1) (car z)
                         (modulo tt 31329)))
               (set! kl
                     (if (>= z-len 2) (cadr z)
                         (modulo tt 30082)))
               (let ((i (+ (modulo (inexact->exact (floor (/ ij 177)))
                                   177) 2))
                     (j (+ (modulo ij 177) 2))
                     (k (+ (modulo (inexact->exact (floor (/ kl 169)))
                                   178) 1))
                     (l (modulo kl 169))
                     (ii 0))
                 (let loop ((ii ii))
                   (if (< ii 97)
                       (let ((s 0.0)
                             (t1 0.5)
                             (jj 0))
                         (let loop ((jj jj))
                           (if (< jj 24)
                               (let ((m (modulo (* (modulo (* i j)  
179) k) 179)))
                                 (set! i j)
                                 (set! j k)
                                 (set! k m)
                                 (set! l (modulo (+ (* 53 l) 1) 169))
                                 (if (not (< (modulo (* l m) 64) 32))
                                     (set! s (+ s t1)))
                                 (set! t1 (* 0.5 t1))
                                 (loop (+ jj 1)))))
                         (vector-set! u ii s)
                         (loop (+ ii 1)))))))))

     (srandom:seed-state)

     (let ((srandom-one
            (lambda ()
              (let ((uni (- (vector-ref u i97)
                            (vector-ref u j97))))
                (if (< uni 0)
                    (set! uni (+ uni 1)))
                (vector-set! u i97 uni)
                (set! i97 (- i97 1))
                (if (< i97 0)
                    (set! i97 96))
                (set! j97 (- j97 1))
                (if (< j97 0) (set! j97 96))
                (set! c (- c cd))
                (if (< c 0) (set! c (+ c cm)))
                (set! uni (- uni c))
                (if (< uni 0) (set! uni (+ uni 1)))
                uni))))

       (set! srandom
             (lambda (n)
               ;; n > 0
               (let ((r (* n (srandom-one))))
                 (if (and (integer? n) (exact? n))
                     (inexact->exact (floor r))
                     r))))))

   )




Posted on the users mailing list.