;;;; moa.scm ;;;; Kon Lovett, Mar '21 ;;;; Kon Lovett, Oct '09 #> #include /* bitsizeof is C++, sigh */ #define bitsizeof( t ) (CHAR_BIT * sizeof( t )) /* Mother ************************************************************** | George Marsaglia's The mother of all random number generators | producing uniformly distributed pseudo random 32 bit values with | period about 2^250. | | The arrays mother1 and mother2 store carry values in their | first element, and random 16 bit numbers in elements 1 to 8. | These random numbers are moved to elements 2 to 9 and a new | carry and number are generated and placed in elements 0 and 1. | The arrays mother1 and mother2 are filled with random 16 bit values. | | Returns: | A 32 bit random number is obtained by combining the output of the | two generators | | Bob Wheeler 8/8/94 */ static uint8_t masktab[ 256 ] = { 0, 1, 3, 3, 7, 7, 7, 7, 15, 15, 15, 15, 15, 15, 15, 15, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255 }; #define M15MASK 0x7FFF /* mask for lower 15 bits */ #define M31MASK 0x7FFFFFFF /* mask for 31 bits */ #define M32DOUBLE 4294967295.0 /* 2^32-1 */ #define STATE_LENGTH 10 #define low16( x ) ((x) & 0xffff) #define high16( x ) (((x) >> 16) & 0xffff) #define make32( x, y ) (((x) << 16) | (y)) #define low32( x ) ((x) & 0xffffffff) #define high32( x ) (((x) >> 32) & 0xffffffff) static void init_state( uint32_t *state, uint64_t seed ) { int i; /* Initialize mother with 9 random values */ uint32_t nseed = low32( seed ); uint64_t nxtval = high32( seed ) & M31MASK; for( i = 0; i < STATE_LENGTH; ++i ) { /* One line multiply-with-carry */ nxtval = 30903 * nseed + high32( nxtval ); state[ i ] = nseed = low32( nxtval ); } /* make carry 31 bits */ state[ 0 ] &= M31MASK; } static uint32_t uniformu32( uint32_t *state ) { # define mom_mov( m, i ) ((m)[ (i) + 1 ] = (m)[ (i) ]) uint64_t nxtval; /* Move elements 1 to 8 to 2 to 9 */ mom_mov( state, 8 ); mom_mov( state, 7 ); mom_mov( state, 6 ); mom_mov( state, 5 ); mom_mov( state, 4 ); mom_mov( state, 3 ); mom_mov( state, 2 ); mom_mov( state, 1 ); /* Form the linear combinations */ nxtval = state[ 0 ] + make32( 1111, 1941) * state[ 2 ] + make32( 2222, 1860 ) * state[ 3 ] + make32( 3333, 1812 ) * state[ 4 ] + make32( 4444, 1776 ) * state[ 5 ] + make32( 5555, 1492 ) * state[ 6 ] + make32( 6666, 1215 ) * state[ 7 ] + make32( 7777, 1066 ) * state[ 8 ] + make32( 9272, 12013 ) * state[ 9 ]; /* Save the high bits of nxtval as the new carry */ state[ 0 ] = high32( nxtval ); /* The low bits of nxtval are the new "seed" */ state[ 1 ] = low32( nxtval ); return state[ 1 ]; # undef mom_mov } static double uniformf64( uint32_t *state ) { return ((((double) uniformu32( state )) / M32DOUBLE) + ((double) uniformu32( state ))) / M32DOUBLE; } static uint32_t randomu32( uint32_t *state, uint32_t m ) { uint32_t r, mask; mask = (m < 256 ? masktab[ m ] : (m < 65536 ? masktab[ m >> 8 ] << 8 | 0xff : (m < 0x1000000 ? masktab[ m >> 16 ] << 16 | 0xffff : masktab[ m >> 24 ] << 24 | 0xffffff))); while( (r = uniformu32( state ) & mask) >= m ); return r; } #undef high32 #undef low32 #undef make32 #undef high16 #undef low16 #undef STATE_LENGTH #undef M32DOUBLE #undef M31MASK #undef M15MASK #undef bitsizeof <# (module moa (;export make-random-source-moa) (import scheme (chicken base) (chicken foreign) (chicken fixnum) (chicken flonum) (srfi 4) (type-checks-numbers integer) random-source entropy-source (only srfi-27-numbers random-large-integer random-large-real native-real-precision?)) ;;; (define init_state (foreign-lambda void "init_state" nonnull-u32vector unsigned-integer64)) (define moa_random_integer (foreign-lambda unsigned-integer32 "randomu32" nonnull-u32vector unsigned-integer32)) (define moa_random_real (foreign-lambda double "uniformf64" nonnull-u32vector)) ;;; ;; (define-constant maximum-unsigned-integer32-flonum 4294967295.0) (cond-expand (64bit (define-constant maximum-unsigned-integer32 4294967295) ) ;MAX (else ;32bit (define-constant maximum-unsigned-integer32 1073741823) ) ) ;32bit most-positive-fixnum (define-constant fpMAX maximum-unsigned-integer32-flonum) ;2^32 - 1 (define-constant LOG2-PERIOD 250) (define-constant eMAX (inexact->exact fpMAX)) ;Create a "bignum" if necessary (define-constant INITIAL-SEED (inexact->exact maximum-unsigned-integer32-flonum)) (define-constant STATE-LENGTH 10) (define-constant INTERNAL-ID 'moa) (define-constant EXTERNAL-ID 'moa) ;; (define (entropy-source-big-integer entropy-source) ;64 bit entropy used as a bit source, not a number source! (modulo (inexact->exact (entropy-source-f64-integer entropy-source)) eMAX) ) (define (make-state) (make-u32vector STATE-LENGTH)) (define (moa-initial-state) (let ((state (make-state))) (init_state state INITIAL-SEED) state ) ) (define (moa-unpack-state state) (cons EXTERNAL-ID (u32vector->list state)) ) (define (moa-pack-state external-state) (unless (moa-external-state? external-state) (error 'moa-pack-state "malformed state" external-state) ) (let ((state (make-state))) (do ((i 0 (fx+ i 1)) (ss (cdr external-state) (cdr ss)) ) ((null? ss) state) (let ((x (car ss))) ;(if {integer?(x) && 0 <= x <= 4294967295} ;(if {integer?(x) && 0 <= x && x <= 4294967295} (if (and (integer? x) (<= 0 x 4294967295)) (u32vector-set! state i x) (error 'moa-pack-state "illegal value" x) ) ) ) ) ) (define (moa-external-state? obj) (and (pair? obj) (eq? EXTERNAL-ID (car obj)) (fx= (fx+ STATE-LENGTH 1) (length obj)) ) ) (define (moa-randomize-state state entropy-source) (init_state state (entropy-source-big-integer entropy-source)) state ) (define (moa-pseudo-randomize-state i j) (let ((state (make-state))) (init_state state (+ i j)) state ) ) (define (moa-random-large state n) (random-large-integer moa_random_integer state fpMAX eMAX n) ) (define (moa-random-real-mp state prec) (random-large-real moa_random_integer state fpMAX eMAX prec) ) ;;; (define (make-random-source-moa) (let ((state (moa-initial-state))) (*make-random-source ; make-random-source-moa ; EXTERNAL-ID ; "George Marsaglia's Mother-Of-All Generator" ; LOG2-PERIOD ; fpMAX ; #f ; (lambda () (moa-unpack-state state) ) ; (lambda (new-state) (set! state (moa-pack-state new-state)) ) ; (lambda (entropy-source) (set! state (moa-randomize-state state entropy-source)) ) ; (lambda (i j) (set! state (moa-pseudo-randomize-state i j)) ) ; (lambda () (lambda (n) (check-positive-integer INTERNAL-ID n 'range) (cond ((and (fixnum? n) (<= n maximum-unsigned-integer32)) (moa_random_integer state n)) (else (moa-random-large state n) ) ) ) ) ; (lambda (prec) (cond ((native-real-precision? prec eMAX) (lambda () (moa_random_real state) ) ) (else (lambda () (moa-random-real-mp state prec) ) ) ) ) ) ) ) ;;; ;;; Module Init ;;; (register-random-source! INTERNAL-ID make-random-source-moa) ) ;module moa