(c-declare #<<END /* This program implements the Mersenne twister algorithm for generation of pseudorandom numbers. The program returns random integers in the range 0 to 2^32-1 (this holds even if a long int is larger than 32 bits). Timing with gcc indicates that it is about twice as fast as the built in rand function. The original code was written by Michael Brundage and has been placed in the public domain. There are a three minor changes here: (1) This comment has been added to the program. (2) Type specifiers (ul) have been appended to constants. (3) A commented out block near the end has been removed. */ #include <stdlib.h> #define MT_LEN 624 struct mt_state { int mt_index; ___U32 mt_buffer[MT_LEN]; }; void init_mt(struct mt_state *mt,unsigned int seed) { int i; srand(seed); for (i = 0; i < MT_LEN; i++) mt->mt_buffer[i] = rand(); mt->mt_index = MT_LEN; // DJB: set this to MT_LEN rather than zero so the first call to mt_random runs the update pass } struct mt_state *create_mt(unsigned int seed) { struct mt_state *mt = malloc(sizeof(struct mt_state)); init_mt(mt,seed); return mt; } ___SCMOBJ free_mt(void *p) { free(p); return ___FIX(___NO_ERR); } #define MT_IA 397 #define MT_IB (MT_LEN - MT_IA) #define UPPER_MASK 0x80000000 #define LOWER_MASK 0x7FFFFFFF #define MATRIX_A 0x9908B0DF #define TWIST(b,i,j) ((b)[i] & UPPER_MASK) | ((b)[j] & LOWER_MASK) #define MAGIC(s) (((s)&1)*MATRIX_A) ___U32 mt_random(struct mt_state *mt) { ___U32 * b = mt->mt_buffer; int idx = mt->mt_index; ___U32 s; int i; if (idx == MT_LEN) { idx = 0; i = 0; for (; i < MT_IB; i++) { s = TWIST(b, i, i+1); b[i] = b[i + MT_IA] ^ (s >> 1) ^ MAGIC(s); } for (; i < MT_LEN-1; i++) { s = TWIST(b, i, i+1); b[i] = b[i - MT_IB] ^ (s >> 1) ^ MAGIC(s); } s = TWIST(b, MT_LEN-1, 0); b[MT_LEN-1] = b[MT_IA-1] ^ (s >> 1) ^ MAGIC(s); } mt->mt_index = idx + 1; return b[idx]; } // Random from 0 to range-1 inclusive ___U32 rnd_int(struct mt_state *mt, ___U32 range) { // NOTE high CANNOT EQUAL 2^32-1 (see below) //ASSERTP(high != 0xFFFFFFFFu); // divided by 2^32 ___U32 r = (___U32)( ((double)mt_random(mt)) / 4294967296.0 * ((double)range) ); // Thanks to floating point rounding at numbers very near to 4294967296.0 (2^32), it is possible to very occasionally // round up to <range> // So, we check for this... if (r==range) { return range-1; } return r; } // Random u8 ___U8 rnd_u8(struct mt_state *mt) { return mt_random(mt) & 0xFF; } /** generates a random number on [0,1]-real-interval */ double rnd_01(struct mt_state *mt) { // divided by 2^32-1 return ((double)mt_random(mt))/4294967295.0; } /** generates a random number on [0,1)-real-interval */ double rnd_01x(struct mt_state *mt) { // divided by 2^32 return ((double)mt_random(mt))/4294967296.0; } /** generates a random number on (0,1)-real-interval */ double rnd_0x1x(struct mt_state *mt) { // divided by 2^32 return (((double)mt_random(mt))+0.5)/4294967296.0; } /** generates a random number on [0,1) with 53-bit resolution*/ double rnd_01x_53(struct mt_state *mt) { ___U32 a=mt_random(mt)>>5; ___U32 b=mt_random(mt)>>6; return ((double)a*67108864.0+(double)b)/9007199254740992.0; } void rnd_u8buf(struct mt_state *mt,___U8* buf,unsigned int count) { unsigned int i; for (i=0; i!=count; ++i) { buf[i] = mt_random(mt) & 0xFF; } } END ) ;; Make a new random seed value based on the current time. (define (new-random-seed) (inexact->exact (modulo (truncate (time->seconds (current-time))) #xFFFFFFF)) ) ;; Construct a new random number generator with an optional seed value. (define (new-rnd #!optional (seed (new-random-seed))) ((c-lambda (unsigned-int) (pointer (struct "mt_state") #f "free_mt") "create_mt") seed) ) ;; Randomize (reseed) an existing rng. (define (rnd-randomize! rng #!optional (seed (new-random-seed))) ((c-lambda ((pointer (struct "mt_state")) unsigned-int32) void "init_mt") rng seed) ) ;; Generate a random integer in [0,range), ie 0..range-1 inclusive. (define (rnd-int range #!optional (rng default-rnd)) ((c-lambda ((pointer (struct "mt_state")) unsigned-int32) unsigned-int32 "rnd_int") rng range) ) ;; Generate a random integer in [0,255] (define (rnd-u8 #!optional (rng default-rnd)) ((c-lambda ((pointer (struct "mt_state"))) unsigned-int8 "rnd_u8") rng) ) ;; Generate a random flonum in [0,1] (define (rnd-01 #!optional (rng default-rnd)) ((c-lambda ((pointer (struct "mt_state"))) double "rnd_01") rng) ) ;; Generate a random flonum in [0,1) (define (rnd-01x #!optional (rng default-rnd)) ((c-lambda ((pointer (struct "mt_state"))) double "rnd_01x") rng) ) ;; Generate a random flonum in (0,1) (define (rnd-0x1x #!optional (rng default-rnd)) ((c-lambda ((pointer (struct "mt_state"))) double "rnd_01") rng) ) ;; Generate a random flonum in [0,1) with 53-bit precision (ie IEEE double mantissa) ;; The other flonum generators only have 32 bit precision. (define (rnd-01x-53 #!optional (rng default-rnd)) ((c-lambda ((pointer (struct "mt_state"))) double "rnd_01x_53") rng) ) ;; Create and fill a u8vector with random numbers. (define (rnd-u8vector count #!optional (rng default-rnd)) (let ((v (make-u8vector count))) ((c-lambda ((pointer (struct "mt_state")) scheme-object unsigned-int) void #<<END ___U8 *buf = ___CAST(___U8*,___BODY(___arg2)); rnd_u8buf(___arg1,buf,___arg3); END ) rng v count) v ) ) ;; The default rng, randomly initialised from the system time. (define default-rnd (new-rnd))