;;; "prng.scm" schlep compiles to prng.c for comparison with Scheme.

#.
#+schlep
(declare-suffixes!
 '((eed (array "unsigned char"))
   (ate (array "unsigned char"))
   (sta (array "unsigned char"))
   (sra (array unsigned))
   ("rgv" (array (array char)))
   ("ean" double)
   ("iff" double)))

(require 'byte)
(require 'printf)
(require 'logical)
#+schlep(require "math")
#+schlep(require "stdio")
#+schlep(require "stdlib")

#+schlep
(define (integer-length n)
  (case n
    ((0) 0)
    ((1) 1)
    ((2 3) 2)
    ((4 5 6 7) 3)
    (else (+ 4 (integer-length (ash n -4))))))

(define (random:chunk sta)
  (cond ((positive? (array-ref sta 258))
	 (array-set! sta 0 258)
	 (printf "random state called reentrantly\n")
	 (exit 1)))
  (array-set! sta 1 258)
  (let* ((idx (logand #xff (+ 1 (array-ref sta 256))))
	 (xtm (array-ref sta idx))
	 (idy (logand #xff (+ (array-ref sta 257) xtm))))
    (array-set! sta idx 256)
    (array-set! sta idy 257)
    (let ((ytm (array-ref sta idy)))
      (array-set! sta xtm idy)
      (array-set! sta ytm idx)
      (let ((ans (array-ref sta (logand #xff (+ ytm xtm)))))
	(array-set! sta 0 258)
	ans))))

(define (random:random modu state)
  (define bitlen (integer-length (+ -1 modu)))
  (define (rnd)
    (do ((bln bitlen (+ -8 bln))
	 (rbs 0 (+ (ash rbs 8) (random:chunk state))))
	((<= bln 7)
	 (if (positive? bln)
	     (set! rbs (logxor (ash rbs bln)
			       (random:chunk state))))
	 (if (< rbs modu) rbs (rnd)))))
  (rnd))

(define (seed->random-state! sta seed)
					; initialize state
  (do ((idx #xff (+ -1 idx)))
      ((negative? idx))
    (array-set! sta idx idx))
					; merge seed into state
  (do ((i 0 (+ 1 i))
       (j 0 (modulo (+ 1 j) seed-len))
       (seed-len (bytes-length seed))
       (k 0))
      ((>= i 256))
    (let ((swp (array-ref sta i)))
      (set! k (logand #xff (+ k (byte-ref seed j) swp)))
      (array-set! sta (array-ref sta k) i)
      (array-set! sta swp k))))

#-schlep
(define (seed->random-state seed)
  (define sta (create-array (Au8 0) (+ 3 256)))
  (seed->random-state! sta seed)
  sta)

(define (square-diff x y)
  (define z (- x y))
  (* z z))

(define (prng! samples modu sta)
  (define sra (create-array (Au32) samples))
  (do ((cnt (+ -1 samples) (+ -1  cnt))
       (num (random:random modu sta) (random:random modu sta))
       (sum 0 (+ sum num)))
      ((negative? cnt)
       (set! sum (+ sum num))
       (let ((mean (/ sum samples)))
	 (do ((cnt (+ -1 samples) (+ -1 cnt))
	      (var2 0 (+ (square-diff mean (array-ref sra cnt)) var2)))
	     ((negative? cnt)
	      (printf "%d / %d = %g +/- %g\n"
		      sum samples mean (sqrt (/ var2 samples)))))))
    (array-set! sra num cnt)))

(define (main argc argv)
  (define n (if (> argc 1)
		(string->number (array-ref argv 1))
		1000))
  (define sta (create-array (Au8 0) (+ 3 256)))
  (seed->random-state! sta "http://swissnet.ai.mit.edu/~jaffer/SLIB.html")
  (prng! n 999 sta)
  0)