Any good random number algorithm?

[xanatose]

I basically need a pseudo random number generator that is fair.  By fair I mean, if the random number is between 0 and 1. Values below 0.5 will have the same chance of appearing than values above 0.5.

A way I found to test this (in pseudo C code) is:

Code Select Expand


double x = 0;
for(;;) {
x += (2.0 * RandomDoubleBetween0and1()) - 1.0;
printf("%f\n", x);
}


Using this is easy to see if a value is growing forever or going down forever.

If x seem to be increasing forever, then the random generator favors values higher than 0.5. While if x seem to be decreasing forever, the algorithm favors values lower than 0.5. If half of the time x seem to be increasing and half of the time seem to be decreasing (of course randomly), then the algorithm is a fair one, if it does not stays in the same area forever.

A perfect one would look like a stock chart. Sometimes going up a long time, sometimes going down a long time, sometimes being undecided (going up and down). but never staying on one of this states forever.

What pseudo random algorithm would you recommend with this property?

[oex]

Antariy has an excellent one here:
http://www.masm32.com/board/index.php?topic=11679.msg122749#msg122749
(There might be a more up to date version in that thread)

You could weight numbers to get more of a 'stock chart feel'.... ie use the random number to signal the change from 0.5 and if number is bit 0-1 (bit 1) then have this have an equal effect to 2147483648 -> 4294967296 (bit 32)....

ie there is a 50% chance of bit 32 being set so if this is set have this 0.5 +/- 1/32 (ie 0.484375 to 0.515625)
However if bit 1 is set have this 0.5 +/- 32/32 (ie 0 or 1)

(This centers on 0.5 however you could just as easily reset center after each change)

EDIT: Note I did not adjust for sign....

Second EDIT: You could use more 'natural' numbers ie Fibonachi sequence to enhance the change range....

[jj2007]

Antariy has an excellent one here:
http://www.masm32.com/board/index.php?topic=11679.msg122749#msg122749

You can find Alex' algo also as Rand() in MasmBasic.

[FORTRANS]

A perfect one would look like a stock chart. Sometimes going up a long time, sometimes going down a long time, sometimes being undecided (going up and down). but never staying on one of this states forever.

What pseudo random algorithm would you recommend with this property?

Hi,

   An algorithm to do that could be:  You generate a random number,
test it to see if you want to continue in the same direction or maybe
change.  If continure use the last direction variable.  If change, get
another random number and use it to set the direction variable.  So,
then adjust the continue/change test to what looks good to you.
90% continue will look quite a bit different than a 10% continue.

   A refinement might be an outer loop that uses a random number
to see if you want to change the continue/change test point.

HTH,

Steve N.

Code Select Expand


Label1
    Get random number
    IF # > Change1 THEN
        Get random number
        Change2 = #
    END IF
    Get random number
    IF # > Change2 THEN
        Get random number
        Direction = # converted to +1 or -1
    END IF
    Current Point = Current Point + Direction
LOOP Label1


[FORTRANS]

Hi,

   Okay, I coded up what I suggested.  Putting a random number
into Change2 tends to mask other effects.  One level of testing
(no "outer loop") with no change of the initial test value, thus shows
more change.

   Oh well.  It needs a bit more control.

Regards,

Steve N.

[MichaelW]

I basically need a pseudo random number generator that is fair.  By fair I mean, if the random number is between 0 and 1. Values below 0.5 will have the same chance of appearing than values above 0.5.

Any "good" pseudo random number generator should be able to do this, at least reasonably. I don't know if a generator that produces numbers between 0 and 1 was part of the goal here, but that's what I used for my test. And as my primary test for "fairness" I used the arithmetic mean of a large number of values

Park-miller-carta.asm:

Code Select Expand


;==============================================================================
    include \masm32\include\advapi32.inc
    includelib \masm32\lib\advapi32.lib
;==============================================================================
    .data
        rand32_seed dd 1
        randfp_div  dq 2147483648
    .code
;==============================================================================

;------------------------------------------------------------------------
; This is Abel's version of a Park-Miller-Carta generator, details here:
; http://www.masm32.com/board/index.php?topic=6558.0
;------------------------------------------------------------------------

OPTION PROLOGUE:NONE
OPTION EPILOGUE:NONE

align 4
rand32 proc base:DWORD

    mov eax, rand32_seed
    mov ecx, 16807              ; a = 7^5
    mul ecx                     ; edx:eax == a*seed == D:A
    mov ecx, 7fffffffh          ; ecx = m

;**** Park-Miller mod
;   div     ecx
;   mov     eax, edx            ; (a*seed) % 7fffffff
;**** Carta alternative mod

    add edx, edx                ; edx = 2*D
    cmp eax, ecx                ; eax = A
    jna @F
    sub eax, ecx                ; if A>m, A = A - m
  @@:
    add eax, edx                ; eax = A + 2*D
    jns @F
    sub eax, ecx                ; If (A + 2*D)>m
;****
  @@:
    mov rand32_seed, eax        ; save new seed

    xor edx, edx
    div DWORD PTR [esp+4]       ; base
    mov eax, edx

    ret 4

rand32  endp

OPTION PROLOGUE:PrologueDef
OPTION EPILOGUE:EpilogueDef

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

;------------------------------------------------------------------
; This is Abel's code modified to return a floating-point value in
; the interval [0,1), at the top of the FPU stack in ST(0), as per
; the normal convention.
;------------------------------------------------------------------

OPTION PROLOGUE:NONE
OPTION EPILOGUE:NONE

align 4
randfp proc

    mov eax, rand32_seed
    mov ecx, 16807              ; a = 7^5
    mul ecx                     ; edx:eax == a*seed == D:A
    mov ecx, 7fffffffh          ; ecx = m

;**** Park-Miller mod
;   div     ecx
;   mov     eax, edx            ; (a*seed) % 7fffffff
;**** Carta alternative mod

    add edx, edx                ; edx = 2*D
    cmp eax, ecx                ; eax = A
    jna @F
    sub eax, ecx                ; if A>m, A = A - m
  @@:
    add eax, edx                ; eax = A + 2*D
    jns @F
    sub eax, ecx                ; If (A + 2*D)>m
;****
  @@:
    mov rand32_seed, eax        ; save new seed

    fild rand32_seed
    fild randfp_div
    fdiv

    ret

randfp  endp

OPTION PROLOGUE:PrologueDef
OPTION EPILOGUE:EpilogueDef


My test app:

Code Select Expand


;==============================================================================
    include \masm32\include\masm32rt.inc
    .686
    include park-miller-carta.asm
;==============================================================================

printf MACRO format:REQ, args:VARARG
    IFNB <args>
        invoke crt_printf, cfm$(format), args
    ELSE
        invoke crt_printf, cfm$(format)
    ENDIF
    EXITM <>
ENDM

;==============================================================================
    .data
        r8 REAL8 ?
        lo REAL8 1.0
        hi REAL8 0.0
    .code
;==============================================================================

;--------------------------------------------------------
; This procedure calculates the arithmetic mean for the
; specified REAL8 array and leaves the result in ST(0).
;--------------------------------------------------------

amean proc pArray:DWORD, count:DWORD

    mov edx, pArray
    mov ecx, count
    dec count
    fldz
  @@:
    fadd REAL8 ptr [edx+ecx*8]
    dec ecx
    jns @B
    fild count
    fdiv

    ret

amean endp

;==============================================================================
start:
;==============================================================================

    N = 100000000

    REPEAT 200
      invoke randfp
      fstp r8
      printf( "%.15f\n", r8 )
    ENDM

    mov esi, alloc( N*8 )

    ;------------------------------------------------------------
    ; Fill a large array with the value 0.5 and use it to verify
    ; that the amean procedure is working correctly.
    ;------------------------------------------------------------

    xor ebx, ebx
    .WHILE ebx < N
        fld10 0.5
        fstp REAL8 ptr [esi+ebx*8]
        inc ebx
    .ENDW
    invoke amean, esi, N
    fstp r8
    printf( "\ntest mean   : %.15f\n", r8 )

    ;---------------------------------------------
    ; Use the randfp procedure to fill the array,
    ; and display the low, high, and mean values.
    ;---------------------------------------------

    xor ebx, ebx
    .WHILE ebx < N
        invoke randfp
        fcom lo
        fstsw ax
        sahf
        jnb @F
        fst lo
      @@:
        fcom hi
        fstsw ax
        sahf
        jna @F
        fst hi
      @@:
        fstp REAL8 ptr [esi+ebx*8]
        inc ebx
    .ENDW
    printf( "\nrandfp low  : %.15f", lo )
    printf( "\nrandfp high : %.15f", hi )
    invoke amean, esi, N
    fstp r8
    printf( "\nrandfp mean : %.15f\n\n", r8 )

    free esi

    inkey "Press any key to exit..."
    exit

;==============================================================================
end start


And the results:

Code Select Expand


0.000007826369256
0.131537788081914
0.755605321843177
0.458650131709874
0.532767237164080
0.218959186226130
0.047044616192579
0.678864716552198
0.679296405520290
0.934692895505577
0.383502077311277
0.519416371826082
0.830965345725417
0.034572110511363
0.053461635019630
0.529700193088502
0.671149383764714
0.007698186207563
0.383415650576353
0.066842237487435
0.417485974263400
0.686772712040693
0.588976642582566
0.930436494294554
0.846166890114546
0.526928777340800
0.091964890714735
0.653918961994350
0.415999356657267
0.701190594118088
0.910320830065757
0.762198039330542
0.262452993541956
0.047464513685554
0.736081884242594
0.328234225977212
0.632638567592949
0.756410485599190
0.991037385072559
0.365338670555502
0.247038885485381
0.982550285756588
0.722660400439054
0.753355834633112
0.651518574450165
0.072685882914811
0.631634717807174
0.884707128163427
0.272709966637194
0.436411405447870
0.766494777519256
0.477731764782220
0.237774433568120
0.274906840175390
0.359264979138970
0.166507200337946
0.486517382785678
0.897656286135316
0.909208101220429
0.060564327519387
0.904653091914952
0.504522894509137
0.516291963402182
0.319032941013575
0.986642111558467
0.493976684752852
0.266144507098943
0.090732894837856
0.947764249518514
0.073749075178057
0.500707094557583
0.384142147842795
0.277081800159067
0.913817441556603
0.529747393447906
0.464445824734867
0.940979953389615
0.050083983689547
0.761514261830598
0.770204546395689
0.827817299868912
0.125365375541151
0.015867701265961
0.688455300871283
0.868247131351382
0.629543417599052
0.736224513966590
0.725411998108029
0.999457878526300
0.888572213239968
0.233194878324866
0.306321830954403
0.351015247870237
0.513273702003062
0.591113582253456
0.845981559716165
0.412080770358443
0.841510639060289
0.269317272119224
0.415394615381956
0.537303975317627
0.467917368281633
0.287212371360511
0.178327703848481
0.153719977010041
0.571654810570180
0.802405726630241
0.033053754363209
0.534449840895832
0.498480118811131
0.955360759515315
0.748292650561780
0.554583847988397
0.890737480949610
0.624849291052669
0.842039612121880
0.159767522476614
0.212751514744014
0.714709967374802
0.130427261814475
0.090990336611867
0.274588147643954
0.002999600954354
0.414293263107538
0.026876290794462
0.709819592535496
0.937897298950702
0.239910804666579
0.180895908735693
0.317539536394179
0.886990661732852
0.652058685664088
0.150335059501231
0.681346213445067
0.385814703535289
0.387725336942822
0.499741032253951
0.147533003240824
0.587186622899026
0.845575659070164
0.590108609758317
0.955408826004714
0.556146138347685
0.148151562083513
0.983305096626282
0.408766693435609
0.141819771379232
0.564898680429906
0.252126406412572
0.488514549098909
0.464030528441072
0.961095140315592
0.126030805986375
0.199757199268788
0.319249673746526
0.629269156139344
0.126712158787996
0.651253741234541
0.621634025592357
0.803072995506227
0.247841758187860
0.476431802846491
0.389314169529825
0.203250334598124
0.028375181369483
0.901673498563468
0.426497412845492
0.142021032050252
0.947486779652536
0.410313035361469
0.131188531406224
0.885648370720446
0.092173629906029
0.162198551930487
0.071063565090299
0.365339028649032
0.253057363443077
0.135109368246049
0.783153168391436
0.455307283904403
0.349524144548923
0.452300169039518
0.808944586664438
0.931674399878830
0.651646054815501
0.215248384047300
0.679592367261648
0.908921884838492
0.250125593971461
0.860859835520387
0.471262328326702
0.505955874919891
0.600393738131970
0.817561482544988
0.755843531806022
0.462244979105890
0.951367449946702
0.632738699670881
0.439330320339650
0.824697386473417

test mean   : 0.500000005000000

randfp low  : 0.000000020023435
randfp high : 0.999999992549419
randfp mean : 0.500029647771980


I started out trying to use the cryptographic service provider random number generator, because it is supposed to generate very good random numbers, but it turned out to be too slow for me to reasonably test large numbers of values.