issue532.reds

Red/System []
	   
int-to-float: func [
	n [integer!]
	return: [float!]
	/local sign shifts less? z
][
	sign: 0
	shifts: 0
	less?: true

;    1. If N is negative, negate it in two's complement. Set the high bit (2^31) of the result.
	if n < 0 [
		; note: FIXME: once bug #224 is fixed
		n: 0 or not n - 1
		sign: 1 << 31
	]

;    2. If N < 2^23, left shift it (multiply by 2) until it is greater or equal to.
	while [n < 00800000h] [
		less?: true
		shifts: shifts + 1
		n: n << 1
	]

;    3. If N >= 2^24, right shift it (unsigned divide by 2) until it is less.
	while [n >= 01000000h] [
		less?: false
		shifts: shifts + 1
		n: n >> 1
	]

;    4. Bitwise AND with ~2^23 (one's complement).
	n: n and not 00800000h

;    5. If it was less, subtract the number of left shifts from 150 (127+23).
	if all [less? 0 < shifts][
		shifts: 150 - shifts
	]

;    6. If it was more, add the number of right shifts to 150.
	if all [not less? 0 < shifts][
		shifts: 150 + shifts
	]
;    7. This new number is the exponent. Left shift it by 23 and add it to the number from step 3.
	shifts: shifts << 23
	;n + shifts

;	hack to convert float32! to float64!
	0.0 + as float32! sign or n + shifts
]

random-mt19937: context [

	;; Period parameters  
	#define STATE-SIZE      624       ;; length of the state-array 
	#define STATE-HALF-SIZE 397       ;; half state-array length
	#define MATRIX_A        9908B0DFh ;; constant vector a 
	#define UPPER_MASK      80000000h ;; most significant w-r bits 
	#define LOWER_MASK      7FFFFFFFh ;; least significant r bits 

	;; the array for the state vector  
	state-array: as int-ptr! allocate STATE-SIZE * size? integer!
	idx: STATE-SIZE + 2  ;; index for state-array results
	                     ;; idx = STATE-SIZE + 2 means array mt is not initialized

	;; initialize array state-array with a seed
	init-genrand: func [ 
		  seed [integer!] 
		  /local i j
	][
		state-array/1: seed and FFFFFFFFh 
		i: 1
		j: 1
		until [
			i: i + 1
			state-array/i: state-array/j >>> 30 xor state-array/j * 1812433253 + j
			;; See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier.
			;; In the previous versions, MSBs of the seed affected  
			;; only MSBs of the array.                        
			;; 2002/01/09 modified by Makoto Matsumoto   
			state-array/i: state-array/i and FFFFFFFFh ;; for >32 bit machines 
			j: j + 1 
			i = STATE-SIZE
		]
		idx: STATE-SIZE + 1
	]

	magic-array: as int-ptr! allocate 2 * size? integer! 
	magic-idx:   1 

	;; generates a random number on [0,FFFFFFFF]-interval
	genrand-int32: func [
		return: [integer!] 
		/local p r s y 
	][
		magic-array/1: 00000000h 
		magic-array/2: MATRIX_A

		if  idx > STATE-SIZE [            ;; generate STATE-SIZE words at one time 
			if  idx = (STATE-SIZE + 2) [  ;; if init-genrand has not been called,
				init-genrand 5489         ;; then a default initial seed is used
			]
			p: 1 
			r: 1 
			while [p < (STATE-SIZE - STATE-HALF-SIZE + 1)][ 
				p: p + 1
				y: (state-array/r and UPPER_MASK) or (state-array/p and LOWER_MASK) 
				s: r + STATE-HALF-SIZE
				magic-idx: (y and 00000001h) + 1
				state-array/r: (state-array/s xor (y >>> 1)) xor magic-array/magic-idx
				r: r + 1 ;; r follows p 
			]
		
			while [p < STATE-SIZE][            
				p: p + 1
				y: (state-array/r and UPPER_MASK) or (state-array/p and LOWER_MASK)
				s: r + STATE-HALF-SIZE - STATE-SIZE
				magic-idx: (y and 00000001h) + 1
				state-array/r: state-array/s xor (y >>> 1) xor magic-array/magic-idx
				r: r + 1
			]
		
			;y: (state-array/STATE-SIZE and UPPER_MASK) or (state-array/1 and LOWER_MASK) 
			y: (state-array/624 and UPPER_MASK) or (state-array/1 and LOWER_MASK) 
			magic-idx: (y and 00000001h) + 1
			;state-array/STATE-SIZE: state-array/STATE-HALF-SIZE xor (y >>> 1) xor magic-array/magic-idx
			state-array/624: state-array/397 xor (y >>> 1) xor magic-array/magic-idx

			idx: 1
		]
  
		y: state-array/idx
		idx: idx + 1

		;; Tempering 
		y: y xor (y >>> 11)
		y: y xor ((y <<  7) and 9D2C5680h)
		y: y xor ((y << 15) and EFC60000h)
		y: y xor (y >>> 18)

		return y
	]

	;; generates a random number on [0,1)-real-interval
	genrand-real2: func [
		return: [float!] 
		/local result [float!] betw [float!]
	][  
	    betw: int-to-float genrand-int32
		((1.0 / 4294967296.0) * betw) ;; divided by 2^32 
	]
]

random-mt19937/init-genrand 1234567890

i: 0
while [i < 10][ 
	print-line random-mt19937/genrand-real2		;; crashes on second call
	i: i + 1
]