I'm new to LISP, and I have a project where it is absolutely necessary
that we have a fast and accurate Random Number Generator. My
supervisor wants me to use the Mersenne Twister, but he wants me to
port it to HIS favorite LISP programming environment (found at
www.lispworks.com).
I've downloaded the CMU MT code, but it doesn't port straight across.
I've having problems with their packages ("EXT" and "Kernel") and with
some other flags?? that I don't really understand (:make-load-form-fun
and :just-dump-it-normally).
I'm really excited to be working in LISP, but I just don't have the
experience to know what changes to make to the code to help it work in
this situation.
Does anyone have any good ideas??
ronny
"ronny bj" <·······@hotmail.com> wrote in message
·································@posting.google.com...
> I'm new to LISP, and I have a project where it is absolutely necessary
> that we have a fast and accurate Random Number Generator. My
> supervisor wants me to use the Mersenne Twister, but he wants me to
> port it to HIS favorite LISP programming environment (found at
> www.lispworks.com).
>
> I've downloaded the CMU MT code, but it doesn't port straight across.
> I've having problems with their packages ("EXT" and "Kernel") and with
> some other flags?? that I don't really understand (:make-load-form-fun
> and :just-dump-it-normally).
>
> I'm really excited to be working in LISP, but I just don't have the
> experience to know what changes to make to the code to help it work in
> this situation.
>
> Does anyone have any good ideas??
>
> ronny
I got an implementation of the Mersenne Twister for Lispworks from Xanalys
(it cost me a support incident), so I know they have it. They may give/sell
it to you, or they may be planning to include it in 4.3, shortly to be
released.
Jim Bushnell
·······@hotmail.com (ronny bj) writes:
> I'm new to LISP, and I have a project where it is absolutely necessary
> that we have a fast and accurate Random Number Generator.
I wonder if "accurate" means "deterministic" in this context.
>>>>> "ronny" == ronny bj <·······@hotmail.com> writes:
ronny> I'm new to LISP, and I have a project where it is absolutely necessary
ronny> that we have a fast and accurate Random Number Generator. My
ronny> supervisor wants me to use the Mersenne Twister, but he wants me to
ronny> port it to HIS favorite LISP programming environment (found at
ronny> www.lispworks.com).
ronny> I've downloaded the CMU MT code, but it doesn't port straight across.
ronny> I've having problems with their packages ("EXT" and "Kernel") and with
ronny> some other flags?? that I don't really understand (:make-load-form-fun
ronny> and :just-dump-it-normally).
[snip]
ronny> Does anyone have any good ideas??
Did you just take rand-mt19937.lisp and try to use that? Well, that
probably won't work. You'll need to extract out the parts of mt19937
yourself.
You'll probably need make-random-state, rand-mt19937-initializer, the
defconstants, random-mt19937-update, and random-chunk. Perhaps a
bunch of other support functions too.
Ray
In article <····························@posting.google.com>,
ronny bj <·······@hotmail.com> wrote:
>I'm new to LISP, and I have a project where it is absolutely necessary
>that we have a fast and accurate Random Number Generator. My
>supervisor wants me to use the Mersenne Twister, but he wants me to
>port it to HIS favorite LISP programming environment (found at
>www.lispworks.com).
I would recommend the multiple recursive random number generator (MRG) of
L'Ecuyer, which is used as the basis for the random number package in
the Scheme SRFI-27 (http://srfi.schemers.org/srfi-27/). It has very good
distribution properties, allows parallel independent streams of random
numbers, and is very fast (up to 3,000,000 random fixnums/second and
4,000,000 random flonums/second on a 1200MHz Athlon, programmed in Scheme).
Brad Lucier
>>>>> "Bradley" == Bradley J Lucier <···@cs.purdue.edu> writes:
Bradley> In article <····························@posting.google.com>,
Bradley> ronny bj <·······@hotmail.com> wrote:
>> I'm new to LISP, and I have a project where it is absolutely necessary
>> that we have a fast and accurate Random Number Generator. My
>> supervisor wants me to use the Mersenne Twister, but he wants me to
>> port it to HIS favorite LISP programming environment (found at
>> www.lispworks.com).
Bradley> I would recommend the multiple recursive random number generator (MRG) of
Bradley> L'Ecuyer, which is used as the basis for the random number package in
Bradley> the Scheme SRFI-27 (http://srfi.schemers.org/srfi-27/). It has very good
Bradley> distribution properties, allows parallel independent streams of random
Bradley> numbers, and is very fast (up to 3,000,000 random fixnums/second and
Bradley> 4,000,000 random flonums/second on a 1200MHz Athlon, programmed in Scheme).
Hmm, CMUCL with mt-19937 will generate about 3.5 million fixnums/second and
5.5 million single-floats/second.
On a fairly slow 300 MHz Ultra 30.
Just a data point. I don't know how mt-19937 compares against MRG.
Ray
In article <··············@edgedsp4.rtp.ericsson.se>,
Raymond Toy <···@rtp.ericsson.se> wrote:
>>>>>> "Bradley" == Bradley J Lucier <···@cs.purdue.edu> writes:
>
> Bradley> In article <····························@posting.google.com>,
> Bradley> ronny bj <·······@hotmail.com> wrote:
> >> I'm new to LISP, and I have a project where it is absolutely necessary
> >> that we have a fast and accurate Random Number Generator. My
> >> supervisor wants me to use the Mersenne Twister, but he wants me to
> >> port it to HIS favorite LISP programming environment (found at
> >> www.lispworks.com).
>
> Bradley> I would recommend the multiple recursive random number
>generator (MRG) of
> Bradley> L'Ecuyer, which is used as the basis for the random number
>package in
> Bradley> the Scheme SRFI-27 (http://srfi.schemers.org/srfi-27/). It
>has very good
> Bradley> distribution properties, allows parallel independent
>streams of random
> Bradley> numbers, and is very fast (up to 3,000,000 random
>fixnums/second and
> Bradley> 4,000,000 random flonums/second on a 1200MHz Athlon,
>programmed in Scheme).
>
>Hmm, CMUCL with mt-19937 will generate about 3.5 million fixnums/second and
>5.5 million single-floats/second.
>
>On a fairly slow 300 MHz Ultra 30.
>
>Just a data point. I don't know how mt-19937 compares against MRG.
Here's the basic code for the uniform-(0,1) RNG I recommended:
; the actual generator
(define (mrg32k3a-random-m1 state)
(declare (flonum))
(let* ((x10 (- (* 1403580.0 (f64vector-ref state 1))
(* 810728.0 (f64vector-ref state 2))))
(y10 (- x10 (* (floor (/ x10 4294967087.0)) 4294967087.0)))
(x20 (- (* 527612.0 (f64vector-ref state 3))
(* 1370589.0 (f64vector-ref state 5))))
(y20 (- x20 (* (floor (/ x20 4294944443.0)) 4294944443.0)))
(dx (- y10 y20))
(dy (- dx (* (floor (/ dx 4294967087.0)) 4294967087.0))))
(f64vector-set! state 5 (f64vector-ref state 4))
(f64vector-set! state 4 (f64vector-ref state 3))
(f64vector-set! state 3 y20)
(f64vector-set! state 2 (f64vector-ref state 1))
(f64vector-set! state 1 (f64vector-ref state 0))
(f64vector-set! state 0 y10)
(f64vector-set! state 9 dy)))
The following paper compares some random number generaters; both MRG32K3A and
MT-19937 pass all tests:
www.iro.umontreal.ca/~lecuyer/myftp/papers/wsc01rng.pdf
Other papers of L'Ecuyer found at
http://www.iro.umontreal.ca/~lecuyer/myftp/papers/
generally state that MT-19937 is a good generator.
I don't know anything about MT-19937; one of the advantages of the generator
above is that it doesn't require 32-bit integers to calculate.
Brad
>>>>> "Bradley" == Bradley J Lucier <···@cs.purdue.edu> writes:
Bradley> Here's the basic code for the uniform-(0,1) RNG I recommended:
Bradley> ; the actual generator
Bradley> (define (mrg32k3a-random-m1 state)
Bradley> (declare (flonum))
Bradley> (let* ((x10 (- (* 1403580.0 (f64vector-ref state 1))
Bradley> (* 810728.0 (f64vector-ref state 2))))
Bradley> (y10 (- x10 (* (floor (/ x10 4294967087.0)) 4294967087.0)))
Bradley> (x20 (- (* 527612.0 (f64vector-ref state 3))
Bradley> (* 1370589.0 (f64vector-ref state 5))))
Bradley> (y20 (- x20 (* (floor (/ x20 4294944443.0)) 4294944443.0)))
Bradley> (dx (- y10 y20))
Bradley> (dy (- dx (* (floor (/ dx 4294967087.0)) 4294967087.0))))
Bradley> (f64vector-set! state 5 (f64vector-ref state 4))
Bradley> (f64vector-set! state 4 (f64vector-ref state 3))
Bradley> (f64vector-set! state 3 y20)
Bradley> (f64vector-set! state 2 (f64vector-ref state 1))
Bradley> (f64vector-set! state 1 (f64vector-ref state 0))
Bradley> (f64vector-set! state 0 y10)
Bradley> (f64vector-set! state 9 dy)))
Here's a Lisp version. To get good speed, I assumed that the arg to
floor is a (double-float 0d0 1d0), but I don't know if that's correct
or not. I also did not look how to initialize the state vector. The
constants above were assumed to be double-floats.
This is hack, so be warned! The inlining is to reduce consing of the
result for the test, since CMUCL has compiler support for mt-19937 and
this makes the comparison a bit more even.
(declaim (optimize (speed 3) (safety 0)))
(declaim (type (simple-array double-float (10)) *s*))
(defvar *s* (make-array 10 :element-type 'double-float :initial-element 0.5d0))
(declaim (inline mrg32k3a-random-m1))
(defun mrg32k3a-random-m1 ()
(declare (optimize (speed 3) (safety 0)))
(let* ((state *s*)
(x10 (- (* 1403580d0 (aref state 1))
(* 810728d0 (aref state 2))))
(y10 (- x10 (* (floor (the (double-float 0d0 1d0) (/ x10 4294967087d0))) 4294967087d0)))
(x20 (- (* 527612d0 (aref state 3))
(* 1370589d0 (aref state 5))))
(y20 (- x20 (* (floor (the (double-float 0d0 1d0) (/ x20 4294944443d0))) 4294944443d0)))
(dx (- y10 y20))
(dy (- dx (* (floor (the (double-float 0d0 1d0) (/ dx 4294967087d0))) 4294967087d0))))
(declare (type (simple-array double-float (10)) state))
(setf (aref state 5) (aref state 4))
(setf (aref state 4) (aref state 3))
(setf (aref state 3) y20)
(setf (aref state 2) (aref state 1))
(setf (aref state 1) (aref state 0))
(setf (aref state 0) y10)
(setf (aref state 9) dy)
(aref state 9)))
(defun gen-ran-mrg (n )
(declare (fixnum n) (optimize (speed 3) (safety 0)))
(let ((y 0d0))
(declare (double-float y))
(dotimes (k n)
(setf y (mrg32k3a-random-m1 )))
(* 3 y)))
On my lowly Ultra-30, this generates about 1.8 million
double-float/sec. mt-19937 generates about 3 million
double-float/sec.
FWIW.
Ray
In article <··············@edgedsp4.rtp.ericsson.se>,
Raymond Toy <···@rtp.ericsson.se> wrote:
>Here's a Lisp version. To get good speed, I assumed that the arg to
>floor is a (double-float 0d0 1d0), but I don't know if that's correct
>or not.
I don't know precisely what this means. If it means that the argument is
between 0d0 and 1d0, then that's not right (and it implies that the result is
always 0 ;-).
> I also did not look how to initialize the state vector. The
>constants above were assumed to be double-floats.
It doesn't matter how the state vector is initialized, as long as they are
not all zero.
>This is hack, so be warned! The inlining is to reduce consing of the
>result for the test, since CMUCL has compiler support for mt-19937 and
>this makes the comparison a bit more even.
I wrote a Gambit-C Scheme version below, but I couldn't get mrg32k3a-random-m1
inlined; of course, all the homogeneous vectors and declarations are not
standard scheme (but the homogeneous vectors are part of a final SRFI).
>On my lowly Ultra-30, this generates about 1.8 million
>double-float/sec. mt-19937 generates about 3 million
>double-float/sec.
Here are my results on a 500 MHz UltraSPARC:
> (time (gen-ran-mrg 1000000))
(time (gen-ran-mrg 1000000))
867 ms real time
870 ms cpu time (840 user, 30 system)
3 collections accounting for 18 ms real time (20 user, 0 system)
32000016 bytes allocated
no minor faults
no major faults
9814170891.
So your compiler is better for this, I couldn't get Gambit-C to eliminate
all the conses and I got only about 1.15 million double-float/sec.
BTW, to get a uniform [0,1] random number requires at least another
multiplication for normalization for mrg32k3a.
Brad
(declare (standard-bindings)(extended-bindings)(block)(not safe)(not interrupts-enabled))
(define *s* (make-f64vector 10 0.5))
(define-macro (FIX . body)
`(let ()
(declare (fixnum))
,@body))
(define-macro (FLOAT . body)
`(let ()
(declare (flonum))
,@body))
(define (mrg32k3a-random-m1 state)
(declare (flonum))
(let* ((x10 (- (* 1403580d0 (f64vector-ref state 1))
(* 810728d0 (f64vector-ref state 2))))
(y10 (- x10 (* (floor (/ x10 4294967087d0)) 4294967087d0)))
(x20 (- (* 527612d0 (f64vector-ref state 3))
(* 1370589d0 (f64vector-ref state 5))))
(y20 (- x20 (* (floor (/ x20 4294944443d0)) 4294944443d0)))
(dx (- y10 y20))
(dy (- dx (* (floor (/ dx 4294967087d0)) 4294967087d0))))
(f64vector-set! state 5 (f64vector-ref state 4))
(f64vector-set! state 4 (f64vector-ref state 3))
(f64vector-set! state 3 y20)
(f64vector-set! state 2 (f64vector-ref state 1))
(f64vector-set! state 1 (f64vector-ref state 0))
(f64vector-set! state 0 y10)
(f64vector-set! state 9 dy)
(f64vector-ref state 9)))
(define (gen-ran-mrg n)
(declare (fixnum))
(do ((i (- n 1) (- i 1))
(y 0.0 (mrg32k3a-random-m1 *s*)))
((< i 0) (FLOAT (* 3.0 y)))))
In article <··········@arthur.cs.purdue.edu>,
Bradley J Lucier <···@cs.purdue.edu> wrote:
>> I also did not look how to initialize the state vector. The
>>constants above were assumed to be double-floats.
>
>It doesn't matter how the state vector is initialized, as long as they are
>not all zero.
Sorry, that's not quite right, they need to be positive integers. All the
requirements are in the papers I cited earlier.
Brad
>>>>> "Bradley" == Bradley J Lucier <···@cs.purdue.edu> writes:
Bradley> In article <··············@edgedsp4.rtp.ericsson.se>,
Bradley> Raymond Toy <···@rtp.ericsson.se> wrote:
>> Here's a Lisp version. To get good speed, I assumed that the arg to
>> floor is a (double-float 0d0 1d0), but I don't know if that's correct
>> or not.
Bradley> I don't know precisely what this means. If it means that the argument is
Bradley> between 0d0 and 1d0, then that's not right (and it implies that the result is
Yes, that's right.
Bradley> always 0 ;-).
Oops. :-)
Ok, I'll go look more closely at the paper and figure this out before
I post any more code.
But it's pretty important to know the range of the floor operation so
that the compiler can generate the best code. And even ffloor might
be more appropriate.
Ray
In article <··············@edgedsp4.rtp.ericsson.se>,
Raymond Toy <···@rtp.ericsson.se> wrote:
>Here's is a Lisp implementation, after reading the paper and looking
>at the reference C code.
>
>On a 400 MHz Ultrasparc, I get 10 million double-float generated in
>3.78 sec for about 2.6 million numbers/sec. On this same machine,
>CMUCL's mt19937 generator gets about 3.2 million.
Very nice. A modified version of the implementation of
http://srfi.schemers.org/srfi-27/
that has all the multiple stream business and a few more abstraction
layers (so it's a bit slower) can be found at
http://www.math.purdue.edu/~lucier/srfi-27
With that implementation I get on my 500 MHz UltraSPARC
banach-129% gsi
loading gambc.scm
Gambit Version 4.0 alpha 8
> (load "test-random.o24")
(time (test-integer 1000000))
1492 ms real time
1460 ms cpu time (1460 user, 0 system)
no collections
no bytes allocated
no minor faults
no major faults
(time (test-real 1000000))
1133 ms real time
1130 ms cpu time (1100 user, 30 system)
192 collections accounting for 125 ms real time (130 user, 0 system)
32000000 bytes allocated
no minor faults
no major faults
Those extra abstraction layers probably aren't worth the slowdown for simpler
applications.
Brad
>>>>> "Bradley" == Bradley J Lucier <···@cs.purdue.edu> writes:
Bradley> (define (mrg32k3a-random-m1 state)
Bradley> (declare (flonum))
Bradley> (let* ((x10 (- (* 1403580d0 (f64vector-ref state 1))
Bradley> (* 810728d0 (f64vector-ref state 2))))
Bradley> (y10 (- x10 (* (floor (/ x10 4294967087d0)) 4294967087d0)))
Bradley> (x20 (- (* 527612d0 (f64vector-ref state 3))
Bradley> (* 1370589d0 (f64vector-ref state 5))))
Bradley> (y20 (- x20 (* (floor (/ x20 4294944443d0)) 4294944443d0)))
Bradley> (dx (- y10 y20))
Bradley> (dy (- dx (* (floor (/ dx 4294967087d0)) 4294967087d0))))
Bradley> (f64vector-set! state 5 (f64vector-ref state 4))
Bradley> (f64vector-set! state 4 (f64vector-ref state 3))
Bradley> (f64vector-set! state 3 y20)
Bradley> (f64vector-set! state 2 (f64vector-ref state 1))
Bradley> (f64vector-set! state 1 (f64vector-ref state 0))
Bradley> (f64vector-set! state 0 y10)
Bradley> (f64vector-set! state 9 dy)
Bradley> (f64vector-ref state 9)))
This differs a bit from the reference code in the paper. Why is y10
this way? In the paper y10 would be test for negative values and made
positive if necessary.
Same for y20.
Ray
Hi ronny bj,
> I'm new to LISP, and I have a project where it is absolutely necessary
> that we have a fast and accurate Random Number Generator. My supervisor
> wants me to use the Mersenne Twister, but he wants me to port it to HIS
> favorite LISP programming environment (found at www.lispworks.com).
Fast and true random number generation is pretty mutually exclusive.
Depending upon your environment (e.g. Linux) you could probably do no
better than to read data in from /dev/random or /dev/urandom. /dev/random
is an extremely high quality but slow source of random data. /dev/urandom
is much faster but still almost certainly better than any pseudo-random
generator you could ever devise.
Regards,
Adam