note
	description: "A Linear congruential random number generator"
	author: "Colin LeMahieu"
	date: "$Date$"
	revision: "$Revision$"
	quote: "To limit the press is to insult a nation; to prohibit reading of certain books is to declare the inhabitants to be either fools or knaves. -  Claude-Adrien Helvetius"

class
	LINEAR_CONGRUENTIAL_RNG

inherit
	RANDOM_NUMBER_GENERATOR
	LIMB_MANIPULATION
	SPECIAL_UTILITY
	SPECIAL_ARITHMETIC
		rename
			add as add_special,
			sub as sub_special,
			mul as mul_special
		export
			{NONE}
				all
		end
	INTEGER_X_DIVISION

create
	make

feature

	make (size: INTEGER)
		require
			size >= 0
			size <= 128
		local
			a: INTEGER_X
			m2exp: INTEGER
			number: STRING_8
		do
			inspect size
			when 0..16 then
				m2exp := 32
				number := "29cf535"
			when 17 then
				m2exp := 34
				number := "A3D73AD"
			when 18 then
				m2exp := 36
				number := "28f825c5"
			when 19 then
				m2exp := 38
				number := "a3dd4cdd"
			when 20 then
				m2exp := 40
				number := "28f5da175"
			when 21..28 then
				m2exp := 56
				number := "AA7D735234C0DD"
			when 29..32 then
				m2exp := 64
				number := "BAECD515DAF0B49D"
			when 33..50 then
				m2exp := 100
				number := "292787EBD3329AD7E7575E2FD"
			when 51..64 then
				m2exp := 128
				number := "48A74F367FA7B5C8ACBB36901308FA85"
			when 65..78 then
				m2exp := 156
				number := "78A7FDDDC43611B527C3F1D760F36E5D7FC7C45"
			when 79..98 then
				m2exp := 196
				number := "41BA2E104EE34C66B3520CE706A56498DE6D44721E5E24F5"
			when 99..100 then
				m2exp := 200
				number := "4E5A24C38B981EAFE84CD9D0BEC48E83911362C114F30072C5"
			when 101..238 then
				m2exp := 256
				number := "AF66BA932AAF58A071FD8F0742A99A0C76982D648509973DB802303128A14CB5"
			end
			create a.make_from_string_base (number, 16)
			randinit_lc_2exp (a, 1, m2exp)
		end

feature

	randinit_lc_2exp (a: READABLE_INTEGER_X; c: NATURAL_32; m2exp: INTEGER)
		local
			seedn: INTEGER
		do
			seedn := bits_to_limbs (m2exp)
			create seed.make (m2exp, a, seedn, c)
		end

	seed: RAND_LC_STRUCT

	randseed (seed_a: READABLE_INTEGER_X)
		local
			seedz: READABLE_INTEGER_X
			seedn: INTEGER
		do
			seedz := seed.seed
			seedn := bits_to_limbs (seed.m2exp)
			fdiv_r_2exp (seedz, seed_a, seed.m2exp)
			seedz.item.fill_with (0, seedz.count, seedz.count - (seedn - seedz.count))
			seedz.count := seedn
		end

	randget (target: SPECIAL [NATURAL_32]; target_offset: INTEGER; nbits: INTEGER)
		local
			rbitpos: INTEGER
			chunk_nbits: INTEGER
			tp: SPECIAL [NATURAL_32]
			tp_offset: INTEGER
			tn: INTEGER
			r2p: SPECIAL [NATURAL_32]
			r2p_offset: INTEGER
			rcy: NATURAL_32
			junk: INTEGER
			last_nbits: INTEGER
			savelimb: NATURAL_32
			carry: CELL [NATURAL_32]
		do
			create carry.put (0)
			chunk_nbits := seed.m2exp // 2
			tn := bits_to_limbs (chunk_nbits)
			create tp.make_filled (0, tn)
			from
				rbitpos := 0
			until
				rbitpos + chunk_nbits > nbits
			loop
				r2p := target
				r2p_offset := target_offset + rbitpos // limb_bits
				if rbitpos \\ limb_bits /= 0 then
					junk := lc (tp, tp_offset)
					savelimb := r2p [r2p_offset]
					lshift (r2p, r2p_offset, tp, tp_offset, tn, rbitpos \\ limb_bits, carry)
					rcy := carry.item
					r2p [r2p_offset] := r2p [r2p_offset].bit_or (savelimb)
					if (chunk_nbits \\ limb_bits + rbitpos \\ limb_bits) > limb_bits then
						r2p [r2p_offset + tn] := rcy
					end
				else
					junk := lc (r2p, r2p_offset)
				end
				rbitpos := rbitpos + chunk_nbits
			end
			if rbitpos /= nbits then
				r2p := target
				r2p_offset := target_offset + rbitpos // limb_bits
				last_nbits := nbits - rbitpos
				tn := bits_to_limbs (last_nbits)
				junk := lc (tp, tp_offset)
				if rbitpos \\ limb_bits /= 0 then
					savelimb := r2p [r2p_offset]
					lshift (r2p, r2p_offset, tp, tp_offset, tn, rbitpos \\ limb_bits, carry)
					rcy := carry.item
					r2p [r2p_offset] := r2p [r2p_offset].bit_or (savelimb)
					if rbitpos + tn * limb_bits - rbitpos \\ limb_bits < nbits then
						r2p [r2p_offset + tn] := rcy
					end
				else
					r2p.copy_data (tp, tp_offset, r2p_offset, tn)
				end
				if nbits \\ limb_bits /= 0 then
					target [target_offset + nbits // limb_bits] := target [target_offset + nbits // limb_bits].bit_and (((0).to_natural_32.bit_not |<< (nbits \\ limb_bits)).bit_not)
				end
			end
		end

	lc (target: SPECIAL [NATURAL_32]; target_offset: INTEGER): INTEGER
		local
			tp: SPECIAL [NATURAL_32]
			tp_offset: INTEGER
			seedp: SPECIAL [NATURAL_32]
			seedp_offset: INTEGER
			ap: SPECIAL [NATURAL_32]
			ap_offset: INTEGER
			ta: INTEGER
			tn: INTEGER
			seedn: INTEGER
			an: INTEGER
			m2exp: INTEGER
			bits: INTEGER
			cy: INTEGER
			xn: INTEGER
			tmp: INTEGER
			i: INTEGER
			x: NATURAL_32
			limb: NATURAL_32
			count: INTEGER
			junk: NATURAL_32
			carry: CELL [NATURAL_32]
		do
			create carry.put (0)
			m2exp := seed.m2exp
			seedp := seed.seed.item
			seedn := seed.seed.count
			ap := seed.a.item
			an := seed.a.count
			ta := an + seedn + 1
			tn := bits_to_limbs (m2exp)
			if ta <= tn then
				tmp := an + seedn
				ta := tn + 1
			end
			create tp.make_filled (0, ta)
			mul_special (tp, tp_offset, seedp, seedp_offset, seedn, ap, ap_offset, an, carry)
			junk := carry.item
			i := seed.cn
			if i /= 0 then
				add_n (tp, tp_offset, tp, tp_offset, seed.cp, 0, i, carry)
				if carry.item /= 0 then
					from
						cy := 0
						limb := 1
					until
						cy /= 0 or limb = 0
					loop
						if i >= tn then
							cy := 1
						else
							x := tp [tp_offset + i]
							limb := x + 1
							tp [tp_offset + i] := limb
							i := i + 1
						end
					end
				end
			end
			tp [tp_offset + m2exp // limb_bits] := tp [tp_offset + m2exp // limb_bits].bit_and ((integer_to_limb (1) |<< m2exp \\ integer_to_limb (limb_bits)) - 1)
			seed.seed.item.copy_data (tp, tp_offset, 0, tn)
			bits := m2exp // 2
			xn := bits // limb_bits
			tn := tn - xn
			if tn > 0 then
				count := bits \\ limb_bits
				if count /= 0 then
					rshift (tp, tp_offset, tp, tp_offset + xn, tn, count, carry)
					junk := carry.item
					target.copy_data (tp, tp_offset, target_offset, xn + 1)
				else
					target.copy_data (tp, tp_offset + xn, target_offset, tn)
				end
			end
			Result := (m2exp + 1) // 2
		end
end