I'm still using Mersenne Twister B.
Anyone willing to provide a randomgen with better distribution?

(When I coded some 3D particle demos a while ago it became apparent that Mersenne B distributes values badly - when I generated random vectors for newly-emitted particles, they formed a "Maltese cross" pattern)...
Posted on 2005-02-17 03:55:44 by Homer
Have you tried the random generator that comes with the MASM32 lib? It works quite well from my tests (compared to the other choices at the time). I used a similar test with direct X to paint pixel patterns (or the lack there of). This one was best... but perhaps your particle streams will yeild something else?

; #########################################################################

;
; Park Miller random number algorithm.
;
; Written by Jaymeson Trudgen (NaN)
; Optimized by Rickey Bowers Jr. (bitRAKE)
;
; #########################################################################

.486 ; create 32 bit code
.model flat, stdcall ; 32 bit memory model
option casemap :none ; case sensitive

.code

; #########################################################################

nrandom PROC base:DWORD

mov eax, nrandom_seed

xor edx, edx
mov ecx, 127773
div ecx
mov ecx, eax
mov eax, 16807
mul edx
mov edx, ecx
mov ecx, eax
mov eax, 2836
mul edx
sub ecx, eax
xor edx, edx
mov eax, ecx
mov nrandom_seed, ecx
div base

mov eax, edx
ret

nrandom ENDP

; #########################################################################

nseed proc TheSeed:DWORD

.data
nrandom_seed dd 12345678
.code

mov eax, TheSeed
mov nrandom_seed, eax

ret

nseed endp

; #########################################################################

end
Posted on 2005-02-17 19:09:11 by NaN
Interesting - I just attempted merging them in the following way:

The "guts" function of the MersenneB procedures is called TBRandom, it generates binary randoms.
I modified your nrandom code to check the Base parameter is NULL - if so, forego the divide-by-base, call TBRandom, and xor the two randoms (ie, temper YOUR random with MINE). Then I replaced (in other procs) any call to TBRandom with a call to nrandom with NULL as Base param.

The nrandom_seed is now being set up from Mersenne's init proc, thus I ditched your initializer proc.

Just to let you know, your code alone has craptastic distribution, with obvious "spikes" that are far away from the "core" values.
I have some 3D code that is generating a bunch of randomly located and colored spheres, with nrandom alone it creates a "cluster" of these objects near the center, with outer distribution being 1) less dense and 2) marked by "spikes" ie clusters of objects that are remote from the "core" (beyond the edge of the distribution falloff).

I did not notice any artifacts produced by pure Mersenne B in this situation, but since I know for a fact there's a distribution problem I went ahead and merged the random generators.

I'm now getting wonderful distribution and I'd be willing to bet that I've increased the "period" of the random stream by a factor of several units.

Here is my source for anyone who wants/needs a better and "more random" RandomGenerator:



;MERSENNE TWISTER B
;Originally implemented in asm by Agner Fog
;Modified for ParkMiller tempering by Homer

;=======================================================
TRandomInit PROTO :DWORD
TBRandom PROTO ;Binary Random
TIRandom PROTO :DWORD,:DWORD ;Integer Random (Min > X Max)
TFRandom PROTO ;Float Random (0.0 > X <1.0)
nrandom PROTO :DWORD ;Integer Random (0 > X < Base)
;=======================================================
TEMPERING EQU 1 ; comment out this line if no tempering

IF 1
; define constants for MT11213A:
N = 351
M = 175
R = 19
MATRIX_A = 0E4BD75F5H
TEMU = 11
TEMS = 7
TEMT = 15
TEML = 17
TEMB = 655E5280H
TEMC = 0FFD58000H

ELSE
; or constants for MT19937:
N = 624
M = 397
R = 31
MATRIX_A = 09908B0DFH
TEMU = 11
TEMS = 7
TEMT = 15
TEML = 18
TEMB = 9D2C5680H
TEMC = 0EFC60000H

ENDIF

LOWER_MASK = (1 SHL R) - 1 ; lower R bits
UPPER_MASK = -1 SHL R ; upper 32-R bits

;=======================================================

.data
MT DD N dup (0) ; history buffer
TEMP DQ 0 ; dual purpose: conversion to float and MT buffer overrun
USE_MMX DB 0 ; 1 if use of MMX registers
PENDXOK DB 0 ; nonzero if number in PENDINGX is OK
MTI DD 0 ; index into MT buffer
LMASK DD LOWER_MASK ; constants
UMASK DD UPPER_MASK
MATA DD MATRIX_A
PENDINGX DQ 0 ; random number ready for submit in upper part
IFDEF TEMPERING
TMB LABEL QWORD ; constants
DD TEMB, TEMB
TMC LABEL QWORD
DD TEMC, TEMC
ENDIF

;===================================================

.code
TRandomInit PROC Seed

MOV EAX, Seed
mov nrandom_seed,eax
XOR ECX, ECX

R10: IMUL EAX, 29943829 ; make random numbers and put them into buffer
DEC EAX
MOV MT[ECX], EAX
ADD ECX, 4
CMP ECX, N*4
JB R10
MOV [MTI], ECX

; check microprocessor support for MMX instructions:
PUSHAD
; detect if CPUID instruction supported by microprocessor:
PUSHFD
POP EAX
MOV EBX, EAX
XOR EAX, 1 SHL 21
PUSH EAX
POPFD
PUSHFD
POP EAX
XOR EAX, EBX
XOR EDX, EDX
AND EAX, 1 SHL 21
JZ SHORT NOCPUID ; CPUID instruction not supported
XOR EAX, EAX
DB 0FH, 0A2H ; CPUID: get number of CPUID functions
TEST EAX, EAX
JZ SHORT NOCPUID ; CPUID function 1 not supported
MOV EAX, 1
DB 0FH, 0A2H ; CPUID
SHR EDX, 23 ; bit 32 = MMX support
AND EDX, 1
NOCPUID:MOV [USE_MMX], DL
POPAD
MOV [PENDXOK], 0
RET 0 ; RET 4 if not _cdecl calling
TRandomInit ENDP

;=======================================================
;Get random integer within given range
TIRandom PROC uses ebx ecx Min:DWORD,Max:DWORD
; CALL TBRandom ; make random number
invoke nrandom, NULL
MOV EDX, Max ; max
MOV ECX, Min ; min
SUB EDX, ECX
JS RERROR ; max < min
INC EDX ; max - min + 1
MUL EDX ; multiply random number by interval and truncate
LEA EAX, [EDX+ECX] ; add min
RET
RERROR: MOV EAX, 80000000H ; error exit
RET
TIRandom ENDP

;=======================================================
TFRandom PROC NEAR ; generate random float
; CALL TBRandom ; random bits
invoke nrandom, NULL
MOV EDX, EAX ; fast conversion to float
SHR EAX, 12
OR EAX, 3FF00000H
SHL EDX, 20
MOV DWORD PTR [TEMP+4], EAX
MOV DWORD PTR [TEMP], EDX
FLD1
FLD QWORD PTR [TEMP] ; partial memory stall here
FSUBR
RET
; This conversion has a partial memory stall for combining two values in TEMP.
; A solution to this problem is to use pre-calculated values, as shown in
; MERSENNEp.ASM, but at the cost of slowing down TIRandom because it must
; prepare a float in case the next call is to TRandom.

TFRandom ENDP

;===============================================================================

TBRandom PROC ; generate random bits
CMP [USE_MMX], 0 ; can we use MMX registers?
MOV ECX, [MTI]
JE R20

; this version uses MMX registers
CMP ECX, N*4
JNB M70 ; buffer is empty, fill it
M40:

IFDEF TEMPERING ; optional tempering
CMP [PENDXOK], 0
JE SHORT M50
MOV EAX, DWORD PTR [PENDINGX+4] ; a tempered random number is allready pending
MOV [PENDXOK], 0
RET

M50: ; tempering of two random numbers in parallel
MTQ = QWORD PTR MT
MOVQ MM0, MTQ[ECX]
MOVQ MM1, MM0
PSRLD MM0, TEMU
PXOR MM0, MM1
MOVQ MM1, MM0
PSLLD MM0, TEMS
PAND MM0, [TMB]
PXOR MM0, MM1
MOVQ MM1, MM0
PSLLD MM0, TEMT
PAND MM0, [TMC]
PXOR MM0, MM1
MOVQ MM1, MM0
PSRLD MM0, TEML
PXOR MM0, MM1
MOVD EAX, MM0 ; return low number
MOVQ [PENDINGX], MM0 ; save high number for next time
EMMS
IF N AND 1 ; check for the last odd one if N is odd
CMP ECX, (N-1)*4
SETB [PENDXOK]
ELSE
MOV [PENDXOK], 1
ENDIF
ADD ECX, 8
MOV [MTI], ECX
RET

ELSE ; no tempering
MOV EAX, MT[ECX]
ENDIF
ADD ECX, 4
MOV [MTI], ECX
RET

M70: ; buffer is empty. Fill it up
MOVD MM3, [UMASK] ; load constants
MOVD MM4, [LMASK]
MOVD MM5, [MATA]
MOV ECX, OFFSET MT ; kk
MOV EDX, OFFSET MT+M*4 ; km
PUNPCKLDQ MM3, MM3 ; duplicate constants
PUNPCKLDQ MM4, MM4
PUNPCKLDQ MM5, MM5

M80: ; kk loop, first part
MOVQ MM1, [ECX] ; mt[kk+1] : mt[kk]
; get misaligned pair (mt[kk+2] : mt[kk+1])
MOVD MM2, [ECX+4] ; 0 : mt[kk+1]
PUNPCKLDQ MM2, [ECX+8] ; mt[kk+2] : mt[kk+1]
IF M AND 1 ; km is misaligned if M is odd. Get misaligned pair
MOVD MM0, [EDX] ; 0 : mt[km]
PUNPCKLDQ MM0, [EDX+4] ; mt[km+1] : mt[km]
ELSE
MOVQ MM0, [EDX] ; mt[km+1] : mt[km]
ENDIF
PAND MM1, MM3 ; & UPPER_MASK
PAND MM2, MM4 ; & LOWER_MASK
POR MM1, MM2 ; y1 : y0
MOVQ MM2, MM1 ; y1 : y0
PSLLD MM1, 31 ; copy bit 0 into all bits
PSRAD MM1, 31 ; -(y & 1)
PAND MM1, MM5 ; & MATRIX_A
PSRLD MM2, 1 ; y >> 1
PXOR MM2, MM1
PXOR MM0, MM2
MOVQ [ECX], MM0 ; result into mt[kk+1] : mt[kk]
ADD EDX, 8
ADD ECX, 8
CMP EDX, OFFSET MT + (N-1)*4
JB M80 ; loop until km wraparound

MOV EAX, [MT] ; copy beginning to end for kk wraparound
IF (M-N) AND 1 ; do the split km pair if M-N is odd
MOV [EDX+4], EAX
MOVD MM0, [EDX]
PUNPCKLDQ MM0, [MT] ; mt[km+1] : mt[km]
SUB EDX, N*4 ; km wraparound
JMP SHORT M95
ELSE
MOV [EDX], EAX
SUB EDX, N*4 ; km wraparound
ENDIF

M90: ; second part of loop, where km has wrapped around
IF (M-N) AND 1 ; km is misaligned if M-N is odd. Get misaligned pair
MOVD MM0, [EDX] ; 0 : mt[km]
PUNPCKLDQ MM0, [EDX+4] ; mt[km+1] : mt[km]
ELSE
MOVQ MM0, [EDX] ; mt[km+1] : mt[km]
ENDIF
M95: MOVQ MM1, [ECX] ; mt[kk+1] : mt[kk]
; get misaligned pair mt[kk+2] : mt[kk+1]
MOVD MM2, [ECX+4] ; 0 : mt[kk+1]
PUNPCKLDQ MM2, [ECX+8] ; mt[kk+2] : mt[kk+1]
PAND MM1, MM3 ; & UPPER_MASK
PAND MM2, MM4 ; & LOWER_MASK
POR MM1, MM2 ; y1 : y0
MOVQ MM2, MM1 ; y1 : y0
PSLLD MM1, 31 ; copy bit 0 into all bits
PSRAD MM1, 31 ; -(y & 1)
PAND MM1, MM5 ; & MATRIX_A
PSRLD MM2, 1 ; y >> 1
PXOR MM2, MM1
PXOR MM0, MM2
MOVQ [ECX], MM0 ; result into mt[kk+1] : mt[kk]
; loop epilog
ADD ECX, 8
ADD EDX, 8
CMP ECX, OFFSET MT + N*4
JB M90

XOR ECX, ECX

IFDEF TEMPERING ; optional tempering
MOV [MTI], ECX
EMMS
JMP M40
ELSE
EMMS
MOV EAX, MT[ECX]
ADD ECX, 4
MOV [MTI], ECX
RET
ENDIF


; this version for old processors without MMX support:
R20: CMP ECX, N*4
JNB SHORT R50 ; buffer is empty, fill it
R40: MOV EAX, MT[ECX]
ADD ECX, 4
MOV [MTI], ECX

IFDEF TEMPERING ; optional tempering
MOV EDX, EAX
SHR EAX, TEMU
XOR EAX, EDX
MOV EDX, EAX
SHL EAX, TEMS
AND EAX, TEMB
XOR EAX, EDX
MOV EDX, EAX
SHL EAX, TEMT
AND EAX, TEMC
XOR EAX, EDX
MOV EDX, EAX
SHR EAX, TEML
XOR EAX, EDX
ENDIF
RET

; fill buffer with random numbers
R50: PUSH EBX
MOV ECX, OFFSET MT
MOV EDX, OFFSET MT + M*4
; kk loop
R60: MOV EAX, [ECX]
MOV EBX, [ECX+4]
AND EAX, UPPER_MASK
AND EBX, LOWER_MASK
OR EAX, EBX
SHR EAX, 1
SBB EBX, EBX
AND EBX, MATRIX_A
XOR EAX, EBX
XOR EAX, [EDX]
MOV [ECX], EAX
ADD EDX, 4
CMP EDX, OFFSET MT + N*4
JB SHORT R70
MOV EAX, [MT]
MOV [EDX], EAX ; copy begin of table to after end to simplify kk+1 wraparound
MOV EDX, OFFSET MT
R70: ADD ECX, 4
CMP ECX, OFFSET MT + N*4
JB R60 ; loop end
XOR ECX, ECX
MOV [MTI], ECX
POP EBX
JMP R40

TBRandom ENDP




; #########################################################################
;
; Park Miller random number algorithm.
;
; Written by Jaymeson Trudgen (NaN)
; Optimized by Rickey Bowers Jr. (bitRAKE)
;
; #########################################################################
.data
nrandom_seed dd 12345678
; #########################################################################
.code

nrandom PROC base:DWORD
mov eax, nrandom_seed
xor edx, edx
mov ecx, 127773
div ecx
mov ecx, eax
mov eax, 16807
mul edx
mov edx, ecx
mov ecx, eax
mov eax, 2836
mul edx
sub ecx, eax
xor edx, edx
mov eax, ecx
mov nrandom_seed, ecx
.if base!=0
div base
mov eax, edx
.else
push eax
CALL TBRandom ; random bits
pop ebx
xor ebx,eax
.endif
ret

nrandom ENDP

; #########################################################################
Posted on 2005-02-17 22:58:47 by Homer
Interesting... I should dig up my old DirectX compare program to see how your new creation weighs...

To be honest, i have no idea what your talking about (spikes etc). But a picture is worth a thousand words (hint hint) ;)

Regards,
:NaN:
Posted on 2005-02-17 23:23:17 by NaN
Yes visually it's easy to see what I mean.
When I refer to "core" values, I mean that MOST of the values generated congregate around zero, which is very obvious when you begin applying min/max constraints to the generated series... the distribution is tighter around zero, with a lower probability of values approaching the threshholds.
When I refer to "spikes", I mean that there are singularities - for certain values near the threshholds, the probability rises sharply.. certain values are generated more often than others near the threshholds.

What we see under your randomgen alone with my 3D visualisation is a high density around zero, with density decreasing towards the threshholds, but with "Artifacts" in the mostly-empy outer region which have consistant values, ie are far from random.

I'll take some screenies on Monday and post them (I'm away for the weekend), but yes by all means dig your 2D visualisation tool out and compare the results yourself.
Posted on 2005-02-18 00:33:02 by Homer
Regarding the Park-Miller and m32lib nrandom functions, this thread might be of interest.
Posted on 2005-02-18 02:33:51 by Jibz
Hmm I didnt realize it was this "crappy"... I cant find my old DirectX test examples. But i did find my origional source (before it was cleaned up).

I know at the time i was using this in my Direct X visual test for developing paterns. Again, I got a very even distribution of pixels and colors, resulting in an overall grey image with no noticable patterns.

Here is the very first source (badly written, but it was 4 years ago ;) )
RAND32 MACRO base:REQ

; Random number generator based on the Real time clock
; and the Park, Miller random number algorithm
;
; Coded by NaN for WIN32ASM
; May 5, 2001
; rev 2.


push ecx
push edx

ifndef __RAND_BY_NAN__
__RAND_BY_NAN__ equ 1

.data?
NaNRand dd ?
.code

db 0fh,31h
shr eax, 2
add eax, 1
mov NaNRand, eax
endif

mov eax, NaNRand
mov edx,0
mov ecx, 127773 ;q
div ecx ; eax == floor( seed / q)
; edx == remainder
SWAP eax, edx
push edx
mov ecx, 16807
mul ecx ; eax = mul of remainder * a
pop edx ; edx == floor of seed/q

SWAP eax, edx
push edx
mov ecx, 2836
mul ecx
pop edx ; edx == mull of rem * a
; eax == mull of seed/q * r

sub edx, eax
mov eax, edx
mov NaNRand, eax ; save next seed
mov ecx, base
mov edx, 0
div ecx
mov eax, edx
pop edx
pop ecx
ENDM


Give it a test if you dont mind and tell me if its better or worse than the MASM32 package. I've been taking it for granted the 'optomizations' had not operational effect on its quality. Perhaps I was wrong to assume this?

Regards,
:NaN:
Posted on 2005-02-18 23:49:47 by NaN
I've uploaded some screen shots to http://www.homer.ultrano.com/Upload/CombinedVersusNRand.rtf

The upper shot shows the combined generators, the lower shows nrandom alone .. camera positions are not the same, and static images don't do justice but all the same..

What I didn't show is a shot of mersenne by itself. I didn't bother because this visualisation tool doesn't show up its poor distribution.
Posted on 2005-02-20 01:41:28 by Homer
I've re-uploaded the RTF :(
Posted on 2005-02-20 07:38:31 by Homer
Man, Thanks for the pic, but i have no idea which is which.. nor how it prooves anything.

Big dots dont tell me much.. but perhaps im missing somthing.

My test generator had hundreds single pixels, randomly gererated x, then y, then a color (out of 16). The overall screen appeared to be static like with no observable patters. It seems to shift around, vibrantly. If you let your eyes go, it would give you a sorta over all greyish look. Which was enough to satisfy me, as grey is the expected out come of a random distribution of all the colors in the spectrum (over a large sampled area such as the screen space).

Thanks anyways...
Regards,
:NaN:
Posted on 2005-02-20 22:03:32 by NaN
evilhomer,

i don't know how you are using them, but if you are happy to gather up a large amount, and read them from a file, for whatever purposes, have a look at:

http://michaelsilk.blogspot.com/2004/10/article-really-random-numbers.html

it's simply really, using "real" random sources from sites around the net to create a nice stream for you.

the code is written in java, but is trivial to adapt to what you like.

hope it helps.

-- Michael
Posted on 2005-02-21 05:49:26 by abc123
In my applications, I fetch "real randoms" from a web source , I hash the values together, then I ror it with the tick count, and use the result as my Seed for the generator code posted earlier.
I'm trying hard to make sure the Seed is truly random, so then we don't have to worry about the properties of random generators quite as much (we can't predict the series anymore).
Posted on 2005-06-01 08:51:48 by Homer
Homer,

Its been a while but I did some work on random sequence generation some time ago using a 256 pixel square pad and used a mouse manually and with a little care it produces reliable random number sequences that appear to pass the ENT random tests written years ago by John Walker. It takes a while to create a sequence of any size but it was no big deal to use a shorter sequence as a seed to create far larger sequences that also test up very well. Back then I used two seperate random generator, NaN's version and a floating point version reseeding each other in a loop.

The testing showed the method to be as good as atomic decay recordings and similar real world sources.
Posted on 2005-06-01 09:09:08 by hutch--
Considering the "Park Miller" algo, isn't it possible that it outputs zero ? In that case it would be very problematic because it would output zero in every next calls.

EDIT : To be more explicit,  if the result of "sub ecx, eax" is equal to zero, the seed "nrandom_seed" and thus, the result, would never vary and stay equal to zero...

; #########################################################################
;
;                    Park Miller random number algorithm.
;
;                      Written by Jaymeson Trudgen (NaN)
;                  Optimized by Rickey Bowers Jr. (bitRAKE)
;
; #########################################################################

      .486                      ; create 32 bit code
      .model flat, stdcall      ; 32 bit memory model
      option casemap :none      ; case sensitive

    .code

; #########################################################################

nrandom PROC base:DWORD

    mov eax, nrandom_seed

    xor edx, edx
    mov ecx, 127773
    div ecx
    mov ecx, eax
    mov eax, 16807
    mul edx
    mov edx, ecx
    mov ecx, eax
    mov eax, 2836
    mul edx
    sub ecx, eax      <- could this be equal to ZERO ??????
    xor edx, edx
    mov eax, ecx
    mov nrandom_seed, ecx
    div base

    mov eax, edx
    ret

nrandom ENDP

; #########################################################################

nseed proc TheSeed:DWORD

    .data
      nrandom_seed dd 12345678
    .code

    mov eax, TheSeed
    mov nrandom_seed, eax

    ret

nseed endp

; #########################################################################

    end


Posted on 2006-04-20 08:28:37 by Axial
Not at all.
I'm only using Park-Miller to "mutate the seed".
Even if the Park-Miller outputs NULL in a given iteration, its merely used as an input permutation for the mersenne matrix.
Posted on 2006-04-20 09:01:39 by Homer
Hi Homer,
I agree but do you think that this equation is mathematically possible if initial SEED != 0 ?

((SEED / 127773) * 2836) - ((SEED mod 127773) * 16807) = 0


PS : I did ask because my maths starts to get *really* crappy :lol:
Posted on 2006-04-20 09:27:59 by Axial
I suppose its mathematically possible, yes, although I didn't check your example.
The selection of the "initial random seed" is something I mention in the comments which accompany that sourcecode. I mention it as a potentially useful feature, not a potentially desastrous condition..
There's roughly a 1 in 4.3 billion chance that the same initial seed will be used, which will yield the same output series, regardless of the fact that I have a feedback mechanism in place.

This belies the fact that even if the initial seed handed to the park-miller stage produces zero as the output seed which it feeds to the mersenne stage, zero is a perfectly legal seed for the mersenne algorithm, it's not a simple linear function like parkmiller is.
Posted on 2006-04-20 10:48:45 by Homer
If seed = 2^31-1

you will get an output of 0. I wonder if that answer the qn.
Posted on 2006-04-22 12:02:04 by roticv
You answered perfectly ! Thank you Rotciv ! ;)
But may I ask you how you did compute that and if your sure about the fact that 7fffffffh is the only seed that returns 0 ?
Posted on 2006-04-25 11:07:32 by Axial
Sorry for the late reply...

Based on your equation of
((SEED / 127773) * 2836) - ((SEED mod 127773) * 16807) = 0


It looks like there is infinite number of solutions for SEED, but of course we have one constrict which is SEED must be 32bit. However, the answer to your question whether there are other solutions, I would say that I am not too sure. I need to do further testing to solve this question.

Furthermore your equation is not exactly correct, as you fail to take note that in the algorithm we only use only the 32bit part of the product.

So I think your equation should have been

((SEED / 127773) * 2836) mod 2^32- ((SEED mod 127773) * 16807) mod 2^32 = 0


As to how I solve it, I started off with


(SEED / 127773) * 2836 = (SEED mod 127773) * 16807

Looking at it carefully, the easiest solution is of course when

it's 16807 * 2836 and hence you have simultaneous equations where

SEED / 127773 = 16807

and

SEED mod 127773 = 2836

Working it out SEED has to be 16807 * 127773 + 2836 which gives you 7fffffffh.
Posted on 2006-04-28 07:53:44 by roticv