;;; EXPERIMENTS WITH LCM DERIVATIONS, BASED ON KNUTH

(declaim (optimize (safety 3) (debug 3) (speed 0)))

(defmacro make-rng (func &key modulus maximum inherit)
  (let ((cmd (gensym)))
	`(lambda (,cmd)
	   (ecase ,cmd
		 (next ,func)
		 ,@(if inherit
			   `(((mod max) (funcall ,inherit ,cmd)))
			   `((mod (eval ,modulus))
				 (max ,(if maximum
						   `(eval ,maximum)
						   `(1- (eval ,modulus))))))))))

(defun rng-next (rng &optional (n 1))
  (loop repeat (1- n) do (funcall rng 'next))
  (funcall rng 'next))

(defun rng-maximum (rng)
  (funcall rng 'max))

(defun rng-modulus (rng)
  (funcall rng 'mod))

 ;;; GENERATORS BASED ON make-general-history-generator

(defun make-general-history-generator (last op factor m)
  (declare (optimize speed (safety 0))
		   (function op) (function factor)
		   (list last) (fixnum m))
  (let* ((vec (coerce last 'vector))
		 (pos 0)
		 (len (length last)))
	(make-rng (let* ((indeces (loop for i from 0 to (length vec)
									collect (mod (+ i pos) len)))
					 (factors (map 'list factor indeces vec))
					 (old (svref vec pos))
					 (next (mod (eval (apply op factors)) m)))
				(setf (svref vec pos) next
					  pos (mod (1+ pos) len))
				old)
			  :modulus m)))

(defun make-linear-comb-generator (last factors m)
  (make-general-history-generator
   last #'+
   (lambda (pos val) (* (elt factors pos) val))
   m))

(defmacro make-term-generator (x0 m (val) &body term)
  `(make-general-history-generator
	(list ,x0) #'identity (lambda (pos ,val)
							(declare (ignore pos))
						   ,@term)
	,m))

(defun make-qcm-generator (x0 d a c m)
  (make-term-generator x0 m (val)
	(+ (* d val val) (* a val) c)))

(defun make-lcm-generator (x0 a c m)
  (make-term-generator x0 m (val)
	(+ (* a val) c)))

(defparameter *default-rng*
  (make-lcm-generator 0 1103515245 12345 (ash 1 31)))

(defun make-add-previous-two-generator (last j k m)
  (assert (< -1 j k (length last)))
  (make-general-history-generator
   last #'+
   (lambda (pos val)
	 (if (or (= pos j) (= pos k)) val 0))
   m))

(defun make-fib-generator (m)
  (make-add-previous-two-generator '(1 1) 0 1 m))

(defun make-smith-klem-generator (last m)
  (make-add-previous-two-generator last 0 (1- (length last)) m))

(defun make-michael-moore-generator (last m)
  (assert (<= 55 (length last)))
  (make-add-previous-two-generator last 24 55 m))

 ;;; HIGHER ORDER GENERATORS

(defun make-bitwise-generator (bit rng)
  (make-rng (loop for i below bit
				  collect (ash (rng-next rng) i)
					into bits
				  finally (return (apply #'logxor bits)))
			:maximum (ash 1 bit)))

(defun make-algo-m-generator (rng-x rng-y k)
  (let ((v (coerce (loop repeat k collect (rng-next rng-x)) 'vector))
		(rng-y* (make-norming-generator (1- k) rng-y)))
	(make-rng (let ((j (rng-next rng-y*)))
				(prog1 (svref v j)
				  (setf (svref v j) (rng-next rng-x))))
			  :inherit rng-x)))

(defun make-algo-m*-generator (rng-x rng-y rng-z k1 k2)
  (let ((mat (make-array (list k1 k2) :element-type 'fixnum)))
	(dotimes (i k1)
	  (dotimes (j k2)
		(setf (aref mat i j) (rng-next rng-x))))
	(make-rng (let* ((i (floor (* k1 (rng-next rng-y))
							   (rng-modulus rng-y)))
					 (j (floor (* k2 (rng-next rng-z))
							   (rng-modulus rng-z)))
					 (old (aref mat i j)))
				(setf (aref mat i j) (rng-next rng-x))
				old)
			  :inherit rng-x)))

(defun make-algo-b-generator (rng k)
  (assert (rng-modulus rng))
  (let* ((v (loop repeat k collect (rng-next rng) into list
				  finally (return (coerce list 'vector))))
		 (y (rng-next rng 2)))
	(make-rng
	 (let* ((j (floor (* k y) (rng-modulus rng)))
			(old (svref v j)))
	   (setf (svref v j) (rng-next rng))
	   old)
	 :inherit rng)))

(defun make-difference-generator (rng-x rng-y m)
  (make-rng (mod (- (rng-next rng-x)
					(rng-next rng-y))
				 m)
			:modulus m))

(defun make-fraction-lcm-generator (u0 a c m)
  (make-rng (setf u0 (mod (+ (* a u0) (/ c m)) 1))
			:modulus 1
			:maximum 1))

(defun make-factored-bit-generator (z rng)
  (make-rng (mod (* z (rng-next rng)) 2)
			:modulus 2))

(defun make-bit-generator (rng)
  (make-factored-bit-generator 1 rng))

 ;;; DISPENSING GENERATOR

(defun make-general-dispensing-generator (rng pred)
  (make-rng (loop for x = (rng-next rng)
				  when (funcall pred) return x)
			:inherit rng))

(defun make-skip-generator (j rng)
  (let ((i 0))
	(make-general-dispensing-generator
	 rng (lambda () (= 0 (setf i (mod (1+ i) j)))))))

(defun make-luescher-generator (rng)
  (let ((j 0))
	(make-general-dispensing-generator
	 rng (lambda () (< (setf j (mod (1+ j) 500)) 55)))))

 ;;; HELPER GENERATORS

(defun next-uniform-val (rng)
  (/ (rng-next rng) (rng-maximum rng)))

(defun make-uniform-generator (rng)
  (make-rng (next-uniform-val rng)))

(defun next-normed-value (max rng)
  (floor (* max (next-uniform-val rng))))

(defun make-norming-generator (max rng)
  (make-rng (next-normed-value max rng)
			:maximum max))

(defun choose-pos (seq rng)
  (next-normed-value (1- (length seq)) rng))

(defun choose-with (seq rng)
  (elt seq (choose-pos seq rng)))

(defun make-bit-chooser-generator (rng-src &optional rng-ch)
  (unless (integerp (rng-maximum rng-src))
	(error "RNG has no known maximum"))
  (let ((rng-ch (or rng-ch rng-src))
		(bits (ceiling (log (rng-maximum rng-src) 2))))
	(make-rng (if (= 0 (logand (ash 1 (next-normed-value bits rng-ch))
							   (rng-next rng-src)))
				  1 0)
			  :maximum 1)))

 ;;; RANDOM RANDOM NUMBER GENERATOR

(defun make-markov-chain-generator (opts state rng)
  (let ((sums (loop for (name . next) in opts
					for sum = (1+ (apply #'+ (mapcar #'cdr next)))
					collect (cons name sum))))
	(lambda ()
	  (when state
		(let* ((next (cdr (assoc state opts :test 'eq)))
			   (chance (next-normed-value (cdr (assoc state sums)) rng)))
		  (loop for (new . weight) in next
				sum weight into prob
				when (>= prob chance)
				  return (prog1 state (setf state new))))))))

(defun run-markov-chain (opts state rng)
  (loop with gen = (make-markov-chain-generator opts state rng)
		for next = (funcall gen)
		while next collect next))

(run-markov-chain '((state-a (state-a . 10) (nil . 1)))
				  'state-a *default-rng*)

(defun generate-random-term (terms cplx var rng)
  (let ((urng (make-uniform-generator rng))
		(terms (cons (list var) terms)))
	(labels ((make-args (args depth)
			   (cond ((null args) nil)
					 ((car args)
					  (cons (loop with arg = (car args)
								  for term = (make-term depth)
								  if (or (eq arg t)
										 (funcall arg term))
									return term)
							(make-args (cdr args) (1- depth))))
					 (t (cons (make-term (1- depth))
							  (make-args (cdr args) (1- depth))))))
			 (make-term (depth)
			   (when (< (rng-next urng) (/ depth cplx))
				 (let* ((term (choose-with terms rng))
						(op (car term))
						(args (cdr term)))
				   (if args (cons op (make-args args depth))
					   op))))
			 (tree-count (item list)
			   (if (atom list) (if (eq item list) 1 0)
				   (+ (tree-count item (car list))
					  (tree-count item (cdr list))))))
	  (loop for term = (make-term (* 2 cplx))
			do (format nil "~A~%" term)
			when (< cplx (tree-count var term)) return term))))

;; (flet ((positive (x) (and (numberp x) (< 0 x))))
;;   (generate-random-term
;;    `((+ t nil) (* nil nil) (mod listp ,#'positive)
;; 	 (-1) (1) (2) (5))
;;    1 'x *default-rng*))

(defun make-random-term-generator (terms cplx x0 m rng)
  (let ((var (gensym)))
	(make-term-generator x0 m (var)
	  (generate-random-term terms cplx var rng))))

 ;;; UTILITY FUNCTIONS

(defmacro with-output-to-file (file &body body)
  `(with-open-file (*standard-output* ,file :direction :output
											:if-exists :supersede
											:if-does-not-exist :create)
	 ,@body))

(defun collect-n-numbers (rng n)
  (loop repeat n collect (rng-next rng)))

(defun count-different-occurances (rng  &optional (n 1))
  (loop with ht = (make-hash-table)
		repeat (calculate-period rng n)
		do (setf (gethash (rng-next rng) ht) t)
		finally (return (hash-table-count ht))))

(defun count-occurances (list)
  (let ((table (make-hash-table)))
	(loop for item in list
		  do (incf (gethash item table 0)))
	(loop for key being the hash-keys of table
			using (hash-value value)
		  collect (cons key value))))

(defun collect-occurances (rng n)
  (sort (count-occurances (collect-n-numbers rng n))
		#'< :key #'car))

(defun calculate-1-period (rng)
  (loop with start = (rng-next rng)
		for i from 0
		when (= start (rng-next rng))
		  return i))

(defun calculate-period (rng &optional (n 1))
  (loop with start = (loop repeat n
						   collect (rng-next rng))
		for i upfrom 0
		for last = (nconc (cdr (copy-list start)) (list (rng-next rng)))
		  then (nconc (cdr last) (list (rng-next rng)))
		when (equal start last)
		  return i))

(defun count-pred-on-k (rng n k pred)
  (loop repeat n
		for last = (loop repeat k collect (rng-next rng))
		  then (nconc (cdr last) (list (rng-next rng)))
		count (apply pred last) into hit
		finally (return (/ hit n))))

(defun write-n-to-stream (n rng)
  (let ((rng (make-norming-generator 255 rng)))
	(loop repeat n do (format t "~d~%" (rng-next rng)))))

(defun write-n-to-file (file n rng)
  (with-output-to-file file
	(write-n-to-stream n rng)))

(defun write-n-bytes-to-file (file n rng)
  (let ((rng* (make-norming-generator (1- (ash 1 8)) rng)))
	(with-open-file (s file :direction :output
							:if-exists :supersede
							:if-does-not-exist :create
							:element-type 'unsigned-byte)
	  (loop repeat n do (write-byte (rng-next rng*) s)))))

(defun write-rgb-image-with-generator (file width height rng)
  (with-output-to-file file
	(format t "P3~%~d ~d~%255~%" width height)
	(write-n-to-stream (* 3 width height) rng)))

(defun write-grayscale-image-with-generator (file width height rng)
  (with-output-to-file file
	(format t "P2~%~d ~d~%255~%" width height)
	(write-n-to-stream (* width height) rng)))

(defun write-bitmap-image-with-generator (file width height rng)
  (with-output-to-file file
	(format t "P1~%~d ~d~%" width height)
	(loop repeat (* width height) do (format t "~d" (rng-next rng)))))

(defparameter *primes* nil)
(defun give-prime-over (n)
  (declare (optimize speed (safety 0))
		   (fixnum n) (list *primes*))
  (cond ((< n 2) 2)
		((and (first *primes*) (< n (first *primes*)))
		 (loop for (p q . rest) on *primes*
			   when (or (not rest) (> p n (1- q))) return p))
		(t (loop for i from 2
				 unless (some (lambda (p) (= (mod i p) 0)) *primes*)
				   do (push i *primes*)
				   and when (< n i) return i))))

(defun spectral-graph (name n dim rng)
  (with-output-to-file (format nil "spectral/~a.dat" name)
	(loop repeat n
		  for last = (loop repeat dim collect (rng-next rng))
			then (nconc (cdr last) (list (rng-next rng)))
		  do (format t "~{~d ~}~%" last))))

(spectral-graph "fib-100-3d" 1000 3 (make-fib-generator 100))

(defun mean (n rng)
  (/ (loop repeat n sum (rng-next rng)) n))

(defun drain-generator (rng n)
  (rng-next rng n) rng)

(defun domain-fill-check (rng &optional (n 1))
  (when (rng-maximum rng)
	(/ (count-different-occurances rng n)
	   (+ 1 (rng-maximum rng)))))

(defun shuffle-with (type seq rng)
  (loop with seq* = (coerce seq 'vector)
		with nrng = (make-norming-generator (1- (length seq)) rng)
		for i downfrom (1- (length seq)) to 1
		for j = (rng-next nrng)
		for x = (svref seq* i)
		for y = (svref seq* j)
		do (setf (svref seq* i) y)
		   (setf (svref seq* j) x)
		finally (return (coerce seq* type))))

;; (shuffle-with 'list '(1 2 3 4 5 6 7 8 9 10) (make-fib-generator 102))

(defun monte-carlo (tries rng)
  (loop with rng* = (make-uniform-generator rng)
		repeat tries
		for x = (rng-next rng*)
		for y = (rng-next rng*)
		count (< (+ (* x x) (* y y)) 1)
		  into hits
		finally (return (* 4 (/ hits tries)))))

(defun monte-carlo-test (tries rng name)
  (loop for n in tries
		for rng* = (make-uniform-generator rng)
		for hits = (with-output-to-file (format nil "monte-carlo/~a-~d.dat" name n)
					 (loop repeat n
						   for x = (rng-next rng*)
						   for y = (rng-next rng*)
						   do (format t "~f ~f~%" x y)
						   count (< (+ (* x x) (* y y)) 1)))
		collect (list n hits (abs (- pi (* 4 (/ hits n)))))))

 ;;; TESTS

;;; LCM
(monte-carlo-test '(10 100 1000 10000 100000) (make-lcm-generator 7 7 7 10) "lcm-eg")
((10 0 3.141592653589793d0)
 (100 0 3.141592653589793d0)
 (1000 0 3.141592653589793d0)
 (10000 0 3.141592653589793d0)
 (100000 0 3.141592653589793d0))

(write-grayscale-image-with-generator "lcm-eg.pgm" 256 256 (make-lcm-generator 7 7 7 10))

(spectral-graph "lcm-eg" 1000 2 (make-lcm-generator 7 7 7 10))

(monte-carlo-test '(10 100 1000 10000 100000) *default-rng* "glibc-rng")
((10 7 0.3415926535897933d0)
 (100 79 0.018407346410207026d0)
 (1000 797 0.04640734641020705d0)
 (10000 7885 0.012407346410206799d0)
 (100000 78627 0.0034873464102069818d0))

(write-grayscale-image-with-generator "lcm-glibc.pgm" 256 256 *default-rng*)

(spectral-graph "lcm-glibc" 1000 2 *default-rng*)

(spectral-graph "fib-100" 1000 3 (make-fib-generator 100))


;;; Fib
(monte-carlo-test '(10 100 1000 10000 100000) (make-fib-generator 100) "fib")
((10 7         0.3415926535897933d0)
 (100 78       0.02159265358979301d0)
 (1000 750     0.14159265358979312d0)
 (10000 7534   0.12799265358979328d0)
 (100000 75336 0.128152653589793d0))

(spectral-graph "fib-10007" 1000 2 (make-fib-generator 10007))
(spectral-graph "fib-100" 1000 2 (make-fib-generator 100))

(write-grayscale-image-with-generator "fib-100.pgm" 256 256 (make-fib-generator 100))
(write-grayscale-image-with-generator "fib-10007.pgm" 256 256 (make-fib-generator 10007))

(calculate-period (make-fib-generator 10007) 2)

;;; MM
(monte-carlo-test '(10 100 1000 10000 100000)
				  (lambda ()
					(make-michael-moore-generator
					 (loop repeat 56 for i from 1 by 2 collect i)
					 1007))
				  "mm-lcm")
((10 10 0.8584073464102069d0)
 (100 87 0.33840734641020687d0)
 (1000 794 0.03440734641020704d0)
 (10000 7966 0.04480734641020678d0)
 (100000 78655 0.00460734641020677d0))

(write-grayscale-image-with-generator "mm-lcm.pgm" 256 256
									  (make-michael-moore-generator
									   (loop repeat 56 for i from 1 by 2 collect i)
									   1007))

(spectral-graph "mm-lcm" 1000 2 (make-michael-moore-generator
								 (loop repeat 56 for i from 1 by 2 collect i)
								 1007))

;;; Algo m
(monte-carlo-test '(10 100 1000 10000 100000)
				  (make-algo-m-generator
				   (make-lcm-generator 5772156649 3141592653 2718281829 (expt 2 35))
				   (make-lcm-generator 1781072418 2718281829 3131592653 (expt 2 35))
				   64)				  
				  "algo-m-lcm")
((10     5     1.1415926535897931d0)
 (100    82    0.1384073464102067d0)
 (1000   788   0.010407346410207019d0)
 (10000  7803  0.02039265358979314d0)
 (100000 78692 0.006087346410206695d0))

(write-grayscale-image-with-generator
 "algo-m-lcm.pgm" 256 256
 (make-algo-m-generator
  (make-lcm-generator 5772156649 3141592653 2718281829 (expt 2 35))
  (make-lcm-generator 1781072418 2718281829 3131592653 (expt 2 35))
  64))

(write-grayscale-image-with-generator "algo-m-glibc.pgm" 256 256
									  (make-algo-m-generator
									   *default-rng*
									   *default-rng*
									   64))

(spectral-graph "algo-m-lcm" 1000 2
				(make-algo-m-generator
				 (make-lcm-generator 5772156649 3141592653 2718281829 (expt 2 35))
				 (make-lcm-generator 1781072418 2718281829 3131592653 (expt 2 35))
				 64))

(monte-carlo-test '(10 100 1000 10000 100000 1000000)
				  (make-algo-m-generator
				   (make-lcm-generator 5772156649 3141592653 2718281829 (expt 2 35))
				   (make-lcm-generator 1781072418 2718281829 3131592653 (expt 2 35))
				   64)
				  "algo-m-lcm")

(write-grayscale-image-with-generator "algo-m-fib.pgm" 60 60
									  (make-algo-m-generator
									   (make-lcm-generator 7 7 7 10)
									   (make-fib-generator 101)
									  50))

(spectral-graph "algo-m-fib" 1000 2
				(make-algo-m-generator
				 (make-fib-generator 100)
				 (make-fib-generator 101)
				 23))

(monte-carlo-test '(10 100 1000 10000 100000)
				  (make-algo-m-generator
				   *default-rng* *default-rng* 25)
				  "algo-m-glibc")
((10 7 0.3415926535897933d0)
 (100 80 0.05840734641020706d0)
 (1000 794 0.03440734641020704d0)
 (10000 7852 7.926535897930798d-4)
 (100000 78479 0.0024326535897931656d0))

(write-grayscale-image-with-generator "algo-m-glibc" 256 256
									  (make-algo-m-generator
									   *default-rng* *default-rng* 25))

(spectral-graph "algo-m-glibc" 1000 2
				(make-algo-m-generator
				 *default-rng* *default-rng* 25))


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;;; Sound
(write-n-bytes-to-file "lcm-simple.wav" 10000 (make-lcm-generator 7 7 7 10))
(write-n-bytes-to-file "lcm-bad.wav" 10000 (make-lcm-generator 0 1 1 1000))
(write-n-bytes-to-file "lcm-glibc.wav" 10000 *default-rng*)
(write-n-bytes-to-file "fib-100.wav" 10000 (make-fib-generator 100))
(write-n-bytes-to-file "fib-101.wav" 10000 (make-fib-generator 101))
(write-n-bytes-to-file "fib-10007.wav" 10000 (make-fib-generator 10007))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(with-output-to-file "fib-period"
  (loop for i from 1 to 2000
		for p = (calculate-period (make-fib-generator i) 2)
		do (format t "~d ~d~%" i p)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(write-to-file "fib-accuracy" s
  (loop for i from 1 to 500
		for diff = (- (/ i 2) (mean 1000 (make-fib-generator i)))
		do (format s "~f~%" diff)))

(let ((rng (make-lcm-generator 0 5 1 8))
	  (test "lcm-m8-2"))
  (spectral-graph (concatenate 'string "spectral-" test ".dat")
				  1000 2 rng))

(spectral-graph "spectral-lcm-2.dat" 1000 2
				(make-norming-generator 100 *default-rng*))
 
(loop for i from 1 to 255
	  do (write-rgb-image-with-generator
		  (format nil "fib-rgb/~d.pbm" i)
		  512 512 i (make-fib-generator i)))

(write-grayscale-image-with-generator
 "glibc-lcm.pgm" 256 256 *default-rng*)

;;; sufficiently random, high period
(write-grayscale-image-with-generator
 "algo-mat-grayscale.pgm" 512 512
 (make-algo-m*-generator
  (make-fib-generator 100)
  (make-smith-klem-generator
   (collect-n-numbers (make-lcm-generator 7 7 7 101) 100)
   10000)
  (calculate-period
   (make-michael-moore-generator
	(collect-n-numbers (make-lcm-generator 7 7 7 101) 100)
	10000)
   100)
  64 64))

;;; very regular pattern
(write-grayscale-image-with-generator
 "algo-mat-fib-grayscale.pgm" 512 512
 (make-algo-m*-generator
  (make-fib-generator 200)
  (make-fib-generator 300)
  (make-fib-generator 100)
  64 64))

;; (monte-carlo-pi-test 1000000 *default-rng*)

;; (monte-carlo-pi-test 1000000
;; 					 (make-algo-m*-generator
;; 					  (make-fib-generator (give-prime-over 300))
;; 					  (make-fib-generator (give-prime-over 200))
;; 					  (make-fib-generator (give-prime-over 100))
;; 					  64 64))

(write-to-file "fib-prime-count" s
  (loop for p in (reverse *primes*)
		do (format s "~d ~d~%" p (count-different-occurances (make-fib-generator p) 2))))

;;; bit-chooser has great domain fill
(domain-fill-check (make-bitwise-generator 8 (make-bit-chooser-generator (make-lcm-generator 7 7 7 10))))
;;; periods >=2 are excelent
(calculate-period (make-bitwise-generator 8 (make-bit-chooser-generator *default-rng*)) 2)

(collect-n-numbers (make-bitwise-generator 5 (make-bit-chooser-generator *default-rng*)) 10)
(collect-n-numbers (make-bitwise-generator 5 (make-factored-bit-generator 3 (make-lcm-generator 7 7 7 13))) 10)
(calculate-period (make-bitwise-generator 5 (make-factored-bit-generator 7 *default-rng*)))

(count-different-occurances (make-fib-generator (give-prime-over 1000)) 2)
(count-different-occurances (make-fib-generator (give-prime-over 2000)) 2)
(count-different-occurances (make-fib-generator (give-prime-over 3000)) 2)

(let ((rng (make-fib-generator (give-prime-over 2018)x)))
  (collect-n-numbers rng (calculate-period rng 2)))

(domain-fill-check (make-fib-generator 100) 2)

(domain-fill-check
 (make-algo-m-generator
  (make-fib-generator (give-prime-over 2000))
  (make-fib-generator (give-prime-over 3000))
  10)
 5)

(write-grayscale-image-with-generator
 "fib-p2000-grayscale.pgm" 54 32
 (make-fib-generator (give-prime-over 50)))

(write-grayscale-image-with-generator
 "algo-m-grayscale.pgm" 512 512
 (make-algo-m-generator
  (make-fib-generator 100)
  (make-michael-moore-generator
   (collect-n-numbers (make-lcm-generator 7 7 7 10) 100)
   10000)
  64))

(write-grayscale-image-with-generator
 "michael-moore-grayscale.pgm" 512 512
 (make-michael-moore-generator
  (collect-n-numbers *default-rng* 100) 27))

(write-bitmap-image-with-generator
 "michael-moore-bitwise.pgm" 512 512
 (make-bit-chooser-generator
  (make-michael-moore-generator
   (shuffle-with 'list (loop for i from 1 to 512 collect i)
				 *default-rng*)
   2)))

(write-grayscale-image-with-generator
 "bit-chooser-grayscale.pgm" 512 512
 (make-bitwise-generator 8 (make-bit-chooser-generator *default-rng*)))

;;; NOTE: "good" generator, still just produces zeros when transformed
;;; into bit generator, and every second value is taken
(write-grayscale-image-with-generator
 "bitwise-grayscale.pgm"
 512 512
 (make-bitwise-generator 8 (make-bit-generator *default-rng*)))

(collect-n-numbers
 (make-skip-generator 1 (make-bit-generator *default-rng*))
 10)

(loop for i from 1 to 100
	  for rng = (make-smith-klem-generator
				 (collect-n-numbers (make-lcm-generator 7 7 7 10) i)
				 i)
	  do (write-rgb-image-with-generator
		  (format nil "sk-rgb/~d.pbm" i)
		  512 512 i rng))

(calculate-period (make-lcm-generator 1 7 7 1000000000000))

(write-rgb-image-with-generator
 "algo-m-rgb.pgm" 512 512 
 (make-algo-m-generator
  (make-smith-klem-generator
   (collect-n-numbers (make-lcm-generator 7 7 7 (give-prime-over 5000)) 50)
   10000)
  (make-smith-klem-generator
   (collect-n-numbers (make-lcm-generator 13 13 13 (give-prime-over 10000)) 50)
   10000)
  64))

(write-rgb-image-with-generator
 "rrng.pgm" 64 64
 (flet ((positive (x) (and (numberp x) (< 0 x))))
   (make-random-term-generator
	`((+ t nil) (* nil nil) (mod listp ,#'positive)
	  (-1) (1) (2) (5))
	4 0 256 *default-rng*)))

(write-rgb-image-with-generator
 "algo-b-rgb.pgm" 512 512 7
 (make-algo-b-generator (make-lcm-generator 1 5 3 8)
						8 4))

(monte-carlo-pi-test 10000 (give-prime-over 200000)
					 (make-michael-moore-generator
					  (loop for i from 1 to 500 by 5 collect i)
					  (give-prime-over 100000)))

(write-to-file "fib-pi-approx.dat"
  (loop for i from 10 to (give-prime-over 512)
		do (format t "~d ~f~%" i
				   (- (monte-carlo-pi-test 1000000 (make-fib-generator i))
					  pi))))

(- (monte-carlo-pi-test 1000000 (make-fib-generator 503)) pi)

;; (calculate-period *default-rng*)

(- (monte-carlo-pi-test 1000000 *default-rng*) pi)


;; (monte-carlo-test '(1000000)
;; 				  (make-michael-moore-generator (loop repeat 56
;; 															for i from 0 by 30
;; 															collect i)
;; 												50300)
;; 				  "mm-hard")

;; (write-to-file "pi-approx" s
;;   (loop for i from 1 to 100
;; 		do (format s "~f~%"
;; 				   (monte-carlo-pi-test
;; 					10000 10000
;; 					(make-algo-m-generator
;; 					 (make-smith-klem-generator
;; 					  (collect-n-numbers (make-lcm-generator 7 7 7 10000) 50)
;; 					  10000)
;; 					 (make-smith-klem-generator
;; 					  (collect-n-numbers (make-lcm-generator 5 5 5 10000) 30)
;; 					  10000)
;; 					 10000 i)))))

(calculate-period
 (make-algo-m-generator
  (make-smith-klem-generator
   (collect-n-numbers (make-lcm-generator 7 7 7 10000) 50)
   10000)
  (make-smith-klem-generator
   (collect-n-numbers (make-lcm-generator 7 7 7 10000) 30)
   10000)
  64) 2) 								; WARNING: takes a while

(mean 10 (make-uniform-generator
		   (make-smith-klem-generator
			(collect-n-numbers (make-lcm-generator 7 7 7 10000) 50)
			100)))

(calculate-period (make-smith-klem-generator
				   (collect-n-numbers (make-lcm-generator 7 7 7 10000) 50)
				   10000))

(write-to-file "less-than-fib" s
  (loop for i from 1 to 100
		do (format s "~f~%" (count-pred-on-k (make-fib-generator 12700) 100000 i #'<))))

(defun iterate-over (list)
  (lambda () (pop list)))
(defvar *prime-iterator* (iterate-over (reverse *primes*)))

(let* ((mod (funcall *prime-iterator*))
	   (rng (make-algo-m-generator
			 (make-lcm-generator 0 5 3 8) ; period 8
			 (make-lcm-generator 0 5 1 17) ; period 8
			 mod)))	
  (cons mod (calculate-period rng)))		; period 15!

(count-different-occurances
 (make-algo-m-generator
  (make-lcm-generator 0 5 3 8)			; period 8
  (make-lcm-generator 0 5 1 8)			; period 8
  4))

(calculate-period
 (make-algo-m-generator
  (make-lcm-generator 0 5 3 8)			   ; period 8
  (make-lcm-generator 0 5 1 8)			   ; period 8
  4)
 4)

(collect-n-numbers (make-lcm-generator 0 5 1 8) 20)

(calculate-period
 (make-algo-b-generator
  (make-lcm-generator 0 5 3 8)			   ; period 8
  4))


(count-different-occurances
 (make-algo-b-generator
  (make-lcm-generator 0 5 3 8)			; period 8
  4))

(calculate-period (make-lcm-generator 0 5 3 8))

;; (let ((vec (vector 3 2 1)))
;;   (eq vec (coerce (coerce vec 'list) 'vector)))


(collect-n-numbers (make-uniform-generator (make-lcm-generator 7 7 7 10)) 10)
(collect-n-numbers (make-fraction-lcm-generator 7 7 7 10) 100)  ; shorter period (7) than lcm
(collect-n-numbers (make-fraction-lcm-generator 7.0 7 7 10) 100)  ; floats change this!

(let ((rng (make-fraction-lcm-generator 7.0 7 7 10)))
  (write-to-file "fraction-test.dat" 
	(loop repeat 1000 do (format t "~a~%" (rng-next rng)))))


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defvar *mitchel-moore* (make-michael-moore-generator (loop repeat 56
															for i from 0 by 3
															collect i)
													  503))


(defvar *algo-m* (make-algo-m-generator (make-fib-generator 100) *mitchel-moore* 60))



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;; MM

(write-grayscale-image-with-generator "mm-1.pgm" 256 256
									  *mitchel-moore*)

;;; MM

(monte-carlo-test '(10 100 1000 10000 100000)
				  *mitchel-moore* "mitchel-moore")

;;; warm
((10 8 0.05840734641020706d0)
 (100 80 0.05840734641020706d0)
 (1000 794 0.03440734641020704d0)
 (10000 7851 0.0011926535897930357d0)
 (100000 78465 0.0029926535897932816d0))

;;; cold 
((10 10 0.8584073464102069d0)
 (100 73 0.2215926535897932d0)
 (1000 773 0.049592653589793034d0)
 (10000 7875 0.008407346410206795d0)
 (100000 78229 0.012432653589792952d0))

;;; ALGO M (warm)

(monte-carlo-test '(10 100 1000 10000 100000)
				  *algo-m* "algo-m")
((10)
 (100)
 (1000)
 (10000)
 (100000))

;;; FIB

(monte-carlo-test '(10 100 1000 10000 100000)
				  (make-fib-generator 100) "fib-100")
((10)
 (100)
 (1000)
 (10000)
 (100000))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(spectral-graph "mm-cold" 150 2 *mitchel-moore*)
(spectral-graph "mm-warm" 150 2 *mitchel-moore*)

(defun best-monte-carlo (rng &optional n)
  (loop for i from 1 to (calculate-period rng n)
		minimizing (abs (- pi (monte-carlo i rng)))))

(loop for a from 1 to 100
	  minimizing
	  (loop for c from 1 to 100
			minimizing
			(loop for i from 1 to 100
				  minimizing (abs (- pi (monte-carlo i (make-lcm-generator 0 a c 100)))))))

(best-monte-carlo (make-lcm-generator 0 23 52 100))