Random Number Generators

My postdetection filter for HD Radio signals uses a random number generator to initially populate its Differential Evolution optimizer and then to update the population by random combination. The program uses sfc32 by Chris Doty-Humphrey. It is fast and has a good reputation for randomness.

The coordinates of each dot are random numbers from sfc32.

This ten-second stereo recording uses a sequence of random numbers. I passed uniformly distributed samples from sfc32 through a function that made them normally distributed like natural random noise.

Spectrum of nearly four minutes of sound.

Code

;	Assemble with NASM and call from FreeBASIC

	global _UN@0, _INI@0
	extern _fb_Timer@0

	section .data align=16

a	dd	0
b	dd	0		; b & c must follow a
c	dd	0
d	dd	1

	section	.text align=16

;	DECLARE FUNCTION un AS LONG
;	sfc32 by Chris Doty-Humphrey.
;	Return 32-bit uniform deviate.

_UN@0:	mov	edx,[b]
	mov	eax,edx
	shr	edx,9
	xor	edx,eax
	add	eax,[a]
	mov	[a],edx
	mov	ecx,[d]
	inc	ecx
	mov	[d],ecx
	add	eax,ecx
	mov	ecx,[c]
	mov	edx,ecx
	lea	ecx,[8*ecx+ecx]
	mov	[b],ecx
	rol	edx,21
	add	edx,eax
	mov	[c],edx
	ret

;	DECLARE SUB ini
;	Call to initialize

_INI@0:	call	_fb_Timer@0	; FreeBASIC TIMER
	fstp	tword [a]	; tbyte for MASM
warmup:	call	_UN@0
	cmp	byte [d],13
	jne	warmup
	ret

More Output Types

These routines provide uniformly distributed floating point numbers and bounded integers.

;	Assemble with NASM and call from FreeBASIC

	global _UNS@0, _UND@0, _UNI@4
	extern _UN@0

	section .data align=16

temp	dw	0,0,0,0,3fffh	; 80-bit floating point

	section	.text align=16

;	DECLARE FUNCTION uns AS SINGLE
;	Return single precision uniform deviate > 0 and < 1.

_UNS@0:	call	_UN@0
	shr	eax,9
	or	eax,3f800000h	; 1->2
	mov	[temp],eax
	fld1
	fld	dword [temp]
	fsubrp			; 0->1
	ret

;	DECLARE FUNCTION und AS DOUBLE
;	Return extended precision uniform deviate > 0 and < 1.

	align 16

_UND@0:	call	_UN@0
	mov	[temp],eax
	call	_UN@0
	or	eax,80000000h	; mantissa MSB (explicit 1)
	mov	[temp+4],eax
	fld1
	fld	tword [temp]	; tbyte for MASM
	fsubrp
	ret

;	DECLARE FUNCTION uni(BYVAL n AS LONG) AS LONG
;	Return integer uniformly distributed over 0 -> n-1.

	align 16

_UNI@4:	call	_UN@0
	mul	dword [esp+4]
	mov	eax,edx
	ret	4

Normal Distribution

This code based on the Leva algorithm described in Numerical Recipes generates a normally distributed random deviate using two or more uniform deviates (2.74 on average).

;	Assemble with NASM and call from FreeBASIC

	global _GAU@0
	extern _UN@0

	section .data align=16

lnc	dq	-2.772588722239781	; -4*ln(2)
a1	dd	1.7156
a2	dd	-0.449871
a3	dd	0.386595
a4	dd	-0.25472
a5	dd	0.196
a6	dd	0.27597
a7	dd	0.27846
b1	dd	1.5
temp	dd	0

	section	.text align=16

;	DECLARE FUNCTION gau AS SINGLE
;	Return Gaussian deviate, mean = 0, std dev = 1

pp:	fcompp
_GAU@0:	call	_UN@0
	shr	eax,9
	or	eax,3f800000h
	mov	[temp],eax
	call	_UN@0
	fld1			; 1
	fld	dword [temp]	; t 1		1->2
	fsubrp			; u		0->1
	shr	eax,9
	or	eax,3f800000h
	mov	[temp],eax
	fld	dword [temp]	; t u		1->2
	fsub	dword [b1]	; t-1.5 u	0.5
	fmul	dword [a1]	; v u
	fld	st1		; u v u
	fadd	dword [a2]	; x v u
	fld	st1		; v x v u
	fabs			; |v| x v u
	fadd	dword [a3]	; y x v u
	fld	st1		; x y x v u
	fmul	st0,st0		; x y x v u
	fxch	st2		; x y x v u
	fmul	dword [a4]	; b y x v u
	fld	st1		; y b y x v u
	fmul	dword [a5]	; a b y x v u
	faddp			; a+b y x v u
	fmulp			; y*(a+b) x v u
	faddp			; q v u
	fcom	dword [a6]
	fnstsw	ax
	fcomp	dword [a7]	; v u
	sahf
	jbe	ok
	fnstsw	ax
	sahf
	ja	pp
	fld	qword [lnc]	; -4*ln(2) v u		this path taken
	fld	st2		; u -4*ln(2) v u	1.2% of the time
	fyl2x			; -4*ln(u) v u		so the slow fyl2x
	fld	st2		; u -4*ln(u) v u	is infrequently
	fmul	st0,st0		; u -4*ln(u) v u	executed
	fld	st2		; v u -4*ln(u) v u
	fmul	st0,st0		; v u -4*ln(u) v u
	fxch			; u v -4*ln(u) v u
	fmulp	st2,st0		; v -4*ln(u)*u v u
	fcompp			; v u
	fnstsw	ax
	sahf
	ja	pp
ok:	fdivrp			; v/u
        ret

Marsaglia MWC

The multiply-with-carry random number generator by George Marsaglia that Numerical Recipes calls B1 is simpler, smaller, and somewhat slower than sfc32. Although I can't distinguish the dot patterns by eye or the white noise by ear, the B1 numbers are likely less random than those from sfc32 when subjected to sophisticated statistical tests. Among the 64-bit Numerical Recipes generators, only MWC has a fast 32-bit implementation.

;	Assemble with NASM and call from FreeBASIC

	global _UN@0, _INI@0
	extern _fb_Timer@0

	section .data align=16

x	dq	0

	section	.text align=16

;	DECLARE FUNCTION un AS LONG
;	64-bit B1 from Numerical Recipes, 3rd Ed, p 347.
;	Lower 32 bits passes Diehard. Period ~ 10^19.
;	Return 32-bit uniform deviate.

_UN@0:	mov	eax,4294957665
	mul	dword [x]
	add	eax,[x+4]
	adc	edx,0
	mov	[x],eax
	mov	[x+4],edx
	ret

;	DECLARE SUB ini
;	Call to initialize

_INI@0:	call	_fb_Timer@0	; FreeBASIC TIMER
	fstp	dword [x]
	ret

I use B1 for short, in-line code when speed doesn't matter and highly random output isn't needed. For example, my harmonic distortion generator uses it to select distorted or undistorted sound while inhibiting long runs with no change that always occur for random binary sequences.

;	Assemble with NASM and call from FreeBASIC

	global _BIT@0, _INI@0
	extern _fb_Timer@0

	section .data align=16

x	dq	0
lo	dd	20000000h
mid	dd	80000000h
thr	dd	80000000h	; first two outputs
last	dd	1		; are unbiased

	section	.text align=16

;	DECLARE FUNCTION bit AS LONG
;	Return semirandom 0 or -1 (no long runs).

_BIT@0:	mov	eax,4294957665
	mul	dword [x]
	add	eax,[x+4]
	adc	edx,0
	mov	[x],eax
	mov	[x+4],edx
	mov	edx,0e0000000h	; high thr
	cmp	eax,[thr]
	cmovb	edx,[lo]	; low thr
	sbb	eax,eax
	cmp	[last],eax
	mov	[last],eax
	cmovne	edx,[mid]	; unbiased thr if
	mov	[thr],edx	; output changes
	ret

;	DECLARE SUB ini
;	Call to initialize

_INI@0:	call	_fb_Timer@0	; FreeBASIC TIMER
	fstp	dword [x]
	ret

June 14, 202188108 MHz