;;;; moa.scm
;;;; Kon Lovett, Oct '09

(module moa

  (;export
    make-random-source-moa)

  (import
    (except scheme <= - / inexact->exact exact->inexact number? floor)
    chicken
    foreign
    srfi-4
    (only numbers <= - / inexact->exact exact->inexact number? floor)
    random-source
    entropy-source
    (only srfi-27-numbers
      check-positive-integer
      random-large-integer random-large-real))

  (require-library
    srfi-4
    numbers
    random-source entropy-source
    srfi-27-numbers)

  (declare (not usual-integrations <= - / exact->inexact inexact->exact))

#>
#include <math.h>

/* 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 MOTHER_SIZE 10

#define low16( x ) ((x) & 0xffff)
#define high16( x ) (((x) >> 16) & 0xffff)

static uint32_t
uniformu32( uint16_t *state )
{
# define mom_mov( m, i ) ((m)[ (i) + 1 ] = (m)[ (i) ])

  uint32_t number1, number2;
  uint16_t *mother1 = state;
  uint16_t *mother2 = &state[ MOTHER_SIZE ];

  /* Move elements 1 to 8 to 2 to 9 */
  mom_mov( mother1, 8 );
  mom_mov( mother1, 7 );
  mom_mov( mother1, 6 );
  mom_mov( mother1, 5 );
  mom_mov( mother1, 4 );
  mom_mov( mother1, 3 );
  mom_mov( mother1, 2 );
  mom_mov( mother1, 1 );
  mom_mov( mother2, 8 );
  mom_mov( mother2, 7 );
  mom_mov( mother2, 6 );
  mom_mov( mother2, 5 );
  mom_mov( mother2, 4 );
  mom_mov( mother2, 3 );
  mom_mov( mother2, 2 );
  mom_mov( mother2, 1 );

  /* Put the carry values in numberi */
  number1 = mother1[ 0 ];
  number2 = mother2[ 0 ];

  /* Form the linear combinations */
  number1 +=    1941 * mother1[ 2 ] + 1860 * mother1[ 3 ]
              + 1812 * mother1[ 4 ] + 1776 * mother1[ 5 ]
              + 1492 * mother1[ 6 ] + 1215 * mother1[ 7 ]
              + 1066 * mother1[ 8 ] + 12013 * mother1[ 9 ];

  number2 +=    1111 * mother2[ 2 ] + 2222 * mother2[ 3 ]
              + 3333 * mother2[ 4 ] + 4444 * mother2[ 5 ]
              + 5555 * mother2[ 6 ] + 6666 * mother2[ 7 ]
              + 7777 * mother2[ 8 ] + 9272 * mother2[ 9 ];

  /* Save the high bits of numberi as the new carry */
  mother1[ 0 ] = high16( number1 );
  mother2[ 0 ] = high16( number2 );

  /* Put the low bits of numberi into motheri[ 1 ] */
  mother1[ 1 ] = low16( number1 );
  mother2[ 1 ] = low16( number2 );

  /* Combine the two 16 bit random numbers into one 32 bit */
  return (((uint32_t) mother1[ 1 ]) << 16) | ((uint32_t) mother2[ 1 ]);

# undef mom_mov
}

static double
uniformf64( uint16_t *state )
{
  return ((((double) uniformu32( state )) / M32DOUBLE) + ((double) uniformu32( state ))) / M32DOUBLE;
}

static uint32_t
randomu32( uint16_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;
}

static void
init_state( uint16_t *state, double seed )
{
  int n;
  uint16_t *p;
  uint16_t *mother1 = state;
  uint16_t *mother2 = &state[ MOTHER_SIZE ];

  /* Initialize mother 1 & 2 with 9 random values */

  uint16_t sNumber = low16( (uint32_t) seed );
  uint32_t number = (((uint32_t) seed) >> 16) & M31MASK;

  p = mother1;
  for( n = MOTHER_SIZE - 1; n--; ) {
    /* One line multiply-with-carry */
    number = 30903 * sNumber + (number >> 16);
    *p++ = sNumber = low16( number );
  }
  p = mother2;
  for( n = MOTHER_SIZE - 1; n--; ) {
    /* One line multiply-with-carry */
    number = 30903 * sNumber + (number >> 16);
    *p++ = sNumber = low16( number );
  }

  /* make carry 15 bits */
  mother1[ 0 ] &= M15MASK;
  mother2[ 0 ] &= M15MASK;
}

#undef MOTHER_SIZE
#undef M32DOUBLE
#undef M31MASK
#undef M15MASK

#undef bitsizeof
<#

(define init_state (foreign-lambda void "init_state" nonnull-u16vector double))

(define moa-random-integer (foreign-lambda unsigned-integer32 "randomu32" nonnull-u16vector unsigned-integer32))
(define moa-random-real (foreign-lambda double "uniformf64" nonnull-u16vector))

;;; fp stuff

(include "fp-extn")

;;;

;;

(define-constant fpMAX 4294967295.0)  ; 2^32 - 1 
(define-constant LOG2-PERIOD 250)

(define MAX (inexact->exact fpMAX)) ; Create a "bignum" if necessary

(define-constant INITIAL-SEED 4294967295.0)

(define-constant STATE-LENGTH 20)

(define-constant EXTERNAL-ID ''marsaglia-moa)

;;

(define (make-state) (make-u16vector STATE-LENGTH))

(define (moa-initial-state)
  (let ((state (make-state)))
    (init_state state INITIAL-SEED)
    state ) )

(define (moa-unpack-state state)
  (cons EXTERNAL-ID (u16vector->list state)) )

(define (moa-pack-state external-state)
  (if (not (and (pair? external-state)
                (eq? EXTERNAL-ID (car external-state))
                (fx= (fx+ STATE-LENGTH 1) (length 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 (and (integer? x) (exact? x) (<= 0 x 65535)) (u16vector-set! state i x)
                (error 'moa-pack-state "illegal value" x) ) ) ) ) ) )

(define (moa-randomize-state state entropy-source)
  (init_state state ((@entropy-source-f64 entropy-source)))
  state )

(define (moa-pseudo-randomize-state i j)
  (let ((state (make-state)))
    (init_state state (exact->inexact (+ i j)))
    state ) )

(define (moa-random-large state n)
  (random-large-integer moa-random-integer state fpMAX MAX n) )

(define (moa-random-real-mp state prec)
  (random-large-real moa-random-integer state fpMAX MAX prec) )

;;;

(define (make-random-source-moa)
  (let ((state (moa-initial-state)))
    (*make-random-source
      ;
      'MOA
      ;
      "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 'moa n 'range)
          (cond ((fixnum? n)  (moa-random-integer state n))
                ; 'n' maybe bignum - must be convertable to "unsigned-integer32"
                ((<= n MAX)   (moa-random-integer state (exact->inexact n)))
                (else         (moa-random-large state n)))))
      ;
      (lambda (prec)
        (cond ((or (not prec)
                   (<= (- (floor (/ 1 prec)) 1) MAX))
                (lambda () (moa-random-real state)))
              (else
                (lambda () (moa-random-real-mp state prec)))))) ) )

;;;
;;; Module Init
;;;

(register-random-source! 'MOA make-random-source-moa)

) ;module moa