;;;; well512.scm
;;;; Kon Lovett, Nov '17

(module well512

(;export
  make-random-source-well512)

(import (except scheme <= inexact->exact exact->inexact number?))

(import chicken foreign)

(import
  srfi-4
  (only numbers <= inexact->exact exact->inexact number?)
  random-source
  entropy-source
  (only srfi-27-numbers
    check-positive-integer
    random-large-integer random-large-real
    native-real-precision?))
(require-library
  srfi-4
  numbers
  random-source entropy-source
  srfi-27-numbers)

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

#>
#include "WELL512a.c"

// uint init[ 16 ]
static void
init_state() {
  InitWELLRNG512a( unsigned int *init );
}

static double
uniformf64( uint32_t *state )
{
  return WELLRNG512a();
}

static uint32_t
randomu32( uint32_t *state, uint32_t m )
{
  return WELLRNG512a (void) ;
}

static void
uniformu64_ith_state( uint32_t *state, uint32_t i )
{
  for( ;i > 0; --i )
    uniformu32( state );
}

static void
uniformu64_jth_offset_state( uint32_t *state, uint32_t j )
{
  uint64_t x = (uint64_t) pow( A, j ) * state[ W ] + state[ C ];
  state[ W ] = low32( x );
  state[ C ] = high32( x );
}

<#

(define init_state (foreign-lambda void "init_state" nonnull-u32vector unsigned-integer64))

(define well512-random-integer (foreign-lambda unsigned-integer32 "randomu32" nonnull-u32vector unsigned-integer32))
(define well512-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 eMAX (inexact->exact fpMAX)) ;Create a "bignum" if necessary

(define-constant INITIAL-SEED maximum-unsigned-integer32-flonum)

(define-constant STATE-LENGTH 10)

(define INTERNAL-ID 'well512)
(define EXTERNAL-ID 'well512)

;;

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

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

(define (well512-unpack-state state)
  (cons EXTERNAL-ID (u32vector->list state)) )

(define (well512-pack-state external-state)
  (unless (well512-external-state? external-state)
    (error 'well512-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) (<= 0 x 4294967295))
          (u32vector-set! state i x)
          (error 'well512-pack-state "illegal value" x) ) ) ) ) )

(define (well512-external-state? obj)
  (and
    (pair? obj)
    (eq? EXTERNAL-ID (car obj))
    (fx= (fx+ STATE-LENGTH 1) (length obj)) ) )

(define (well512-randomize-state state entropy-source)
  (init_state
    state
    (exact->inexact
      (modulo
        (fpabs (entropy-source-f64-integer entropy-source))
        fpMAX)))
  state )

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

(define (well512-random-large state n)
  (random-large-integer well512-random-integer state fpMAX eMAX n) )

(define (well512-random-real-mp state prec)
  (random-large-real well512-random-integer state fpMAX eMAX prec) )

;;;

(define (make-random-source-well512)
  (let ((state (well512-initial-state)))
    (*make-random-source
      ;
      make-random-source-well512
      ;
      EXTERNAL-ID
      ;
      "Well's 512-bit Generator"
      ;
      LOG2-PERIOD
      ;
      fpMAX
      ;
      #f
      ;
      (lambda ()
        (well512-unpack-state state) )
      ;
      (lambda (new-state)
        (set! state (well512-pack-state new-state)) )
      ;
      (lambda (entropy-source)
        (set! state (well512-randomize-state state entropy-source)) )
      ;
      (lambda (i j)
        (set! state (well512-pseudo-randomize-state i j)) )
      ;
      (lambda ()
        (lambda (n)
          (check-positive-integer INTERNAL-ID n 'range)
          (cond-expand
            (64bit
              (cond
                ((and (fixnum? n) (<= n maximum-unsigned-integer32))
                  (well512-random-integer state n))
                (else
                  (well512-random-large state n) ) ) )
            (else ;32bit
              (cond
                ((and (fixnum? n) (<= n maximum-unsigned-integer32))
                  (well512-random-integer state n))
                ;'n' maybe bignum - must be convertable to "unsigned-integer32"
                ((<= n eMAX)
                  (well512-random-integer state (exact->inexact n)))
                (else
                  (well512-random-large state n) ) ) ) ) ) )
      ;
      (lambda (prec)
        (cond
          ((native-real-precision? prec eMAX)
            (lambda ()
              (well512-random-real state) ) )
          (else
            (lambda ()
              (well512-random-real-mp state prec) ) ) ) ) ) ) )

;;;
;;; Module Init
;;;

(register-random-source! INTERNAL-ID make-random-source-well512)

) ;module well512