Random number generators ( MWC256 & CMWC4096 )

[wilbert]

mwc256.pbi :

Code: Select all

; *********************************************
; * A PureBasic version of the                *
; *                                           *
; * MWC256 PRNG algorithm by George Marsaglia *
; * (period = 2^8222)                         *
; *                                           *
; *********************************************

Macro m_mwc256
  !mov eax, 809430660
  !movzx ecx, byte [mwc256index]
  !inc cl
  !mov dword [mwc256index], ecx
  CompilerIf #PB_Compiler_Processor = #PB_Processor_x86
    !mul dword [mwc256state + ecx * 4]
    !add eax, dword [mwc256carry]
    !adc edx, 0
    !mov dword [mwc256state + ecx * 4], eax
  CompilerElse
    !lea r8, [mwc256state]
    !mul dword [r8 + rcx * 4]
    !add eax, dword [mwc256carry]
    !adc edx, 0
    !mov dword [r8 + rcx * 4], eax
  CompilerEndIf
  !mov dword [mwc256carry], edx
EndMacro

DataSection
  !mwc256state: times 256 dd 0
  !mwc256carry: dd 362436
  !mwc256index: dd 255
EndDataSection

; *** mwc256seed(seed.l) ***

Procedure mwc256seed(seed.l)
  EnableASM
  MOV eax, seed
  DisableASM
  CompilerIf #PB_Compiler_Processor = #PB_Processor_x64
    !lea r8, [mwc256state]
  CompilerEndIf
  !mov dword [mwc256state], eax
  !mov ecx, 1
  !mwc256seed_loop:
  !mov edx, eax
  !shr edx, 30
  !xor eax, edx
  !mov edx, 1812433253
  !mul edx
  !add eax, ecx
  CompilerIf #PB_Compiler_Processor = #PB_Processor_x86
    !mov dword [mwc256state + ecx * 4], eax
  CompilerElse
    !mov dword [r8 + rcx * 4], eax
  CompilerEndIf
  !inc cl
  !jnz mwc256seed_loop
  !xor edx, edx
  !mov ecx, 61137367 
  !div ecx
  !mov dword [mwc256carry], edx
  !mov dword [mwc256index], 255
EndProcedure

mwc256seed(Date())

; *** mwc256() - get random long ***

Procedure mwc256()
  m_mwc256
  ProcedureReturn
EndProcedure

; *** mwc256range(min.l, max.l) - random [min - max] ***

; requirements : min >= 0, max <= 2147483647, min < max  
Procedure mwc256range(min.l, max.l)
  m_mwc256
  EnableASM
  MOV ecx, min
  MOV edx, max
  DisableASM
  !sub edx, ecx
  !inc edx
  !mul edx
  !mov eax, edx
  !add eax, ecx
  ProcedureReturn
EndProcedure

; *** mwc256f() - get random float [0 - 1) ***

Procedure.f mwc256f()
  Static result.f
  m_mwc256
  !movd xmm0, eax
  !psrlq xmm0, 9
  !mov eax, 1
  !cvtsi2ss xmm1, eax
  !por xmm0, xmm1
  !subss xmm0, xmm1
  !movd [s_mwc256f.v_result], xmm0
  ProcedureReturn result
EndProcedure

; *** mwc256d() - get random double [0 - 1) ***

Procedure.d mwc256d()
  Static result.d
  m_mwc256
  !movd xmm0, eax
  m_mwc256
  !movd xmm1, eax
  !punpckldq xmm0, xmm1
  !psrlq xmm0, 12
  !mov eax, 1
  !cvtsi2sd xmm1, eax
  !por xmm0, xmm1
  !subsd xmm0, xmm1
  !movq [s_mwc256d.v_result], xmm0
  ProcedureReturn result
EndProcedure

[wilbert]

cmwc4096.pbi :

Code: Select all

; ***********************************************
; * A PureBasic version of the                  *
; *                                             *
; * CMWC4096 PRNG algorithm by George Marsaglia *
; * (period = 2^131104)                         *
; *                                             *
; ***********************************************

Macro m_cmwc4096
  !mov eax, 18782
  !mov ecx, dword [cmwc4096index]
  !inc ecx
  !and ecx, 4095
  !mov dword [cmwc4096index], ecx
  CompilerIf #PB_Compiler_Processor = #PB_Processor_x86
    !mul dword [cmwc4096state + ecx * 4]
  CompilerElse
    !lea r8, [cmwc4096state]
    !mul dword [r8 + rcx * 4]
  CompilerEndIf
  !add eax, dword [cmwc4096carry]
  !adc edx, 0
  !xor ecx, ecx
  !add eax, edx
  !setc cl
  !add eax, ecx
  !add edx, ecx
  !mov dword [cmwc4096carry], edx
  !neg eax
  !add eax, 0xfffffffe
  !mov ecx, dword [cmwc4096index]  
  CompilerIf #PB_Compiler_Processor = #PB_Processor_x86
    !mov dword [cmwc4096state + ecx * 4], eax
  CompilerElse
    !mov dword [r8 + rcx * 4], eax
  CompilerEndIf
EndMacro

DataSection
  !cmwc4096state: times 4096 dd 0
  !cmwc4096carry: dd 362436
  !cmwc4096index: dd 4095
EndDataSection

; *** cmwc4096seed(seed.l) ***

Procedure cmwc4096seed(seed.l)
  EnableASM
  MOV eax, seed
  DisableASM
  CompilerIf #PB_Compiler_Processor = #PB_Processor_x64
    !lea r8, [cmwc4096state]
  CompilerEndIf
  !mov dword [cmwc4096state], eax
  !mov ecx, 1
  !cmwc4096seed_loop:
  !mov edx, eax
  !shr edx, 30
  !xor eax, edx
  !mov edx, 1812433253
  !mul edx
  !add eax, ecx
  CompilerIf #PB_Compiler_Processor = #PB_Processor_x86
    !mov dword [cmwc4096state + ecx * 4], eax
  CompilerElse
    !mov dword [r8 + rcx * 4], eax
  CompilerEndIf
  !inc ecx
  !and ecx, 4095
  !jnz cmwc4096seed_loop
  !xor edx, edx
  !mov ecx, 61137367 
  !div ecx
  !mov dword [cmwc4096carry], edx
  !mov dword [cmwc4096index], 4095
EndProcedure

cmwc4096seed(Date())

; *** cmwc4096() - get random long ***

Procedure cmwc4096()
  m_cmwc4096
  ProcedureReturn
EndProcedure

; *** cmwc4096range(min.l, max.l) - random [min - max] ***

; requirements : min >= 0, max <= 2147483647, min < max  
Procedure cmwc4096range(min.l, max.l)
  m_cmwc4096
  EnableASM
  MOV ecx, min
  MOV edx, max
  DisableASM
  !sub edx, ecx
  !inc edx
  !mul edx
  !mov eax, edx
  !add eax, ecx
  ProcedureReturn
EndProcedure

; *** cmwc4096f() - get random float [0 - 1) ***

Procedure.f cmwc4096f()
  Static result.f
  m_cmwc4096
  !movd xmm0, eax
  !psrlq xmm0, 9
  !mov eax, 1
  !cvtsi2ss xmm1, eax
  !por xmm0, xmm1
  !subss xmm0, xmm1
  !movd [s_cmwc4096f.v_result], xmm0
  ProcedureReturn result
EndProcedure

; *** cmwc4096d() - get random double [0 - 1) ***

Procedure.d cmwc4096d()
  Static result.d
  m_cmwc4096
  !movd xmm0, eax
  m_cmwc4096
  !movd xmm1, eax
  !punpckldq xmm0, xmm1
  !psrlq xmm0, 12
  !mov eax, 1
  !cvtsi2sd xmm1, eax
  !por xmm0, xmm1
  !subsd xmm0, xmm1
  !movq [s_cmwc4096d.v_result], xmm0
  ProcedureReturn result
EndProcedure

[wilbert]

demo :

Code: Select all

XIncludeFile "mwc256.pbi"
   
If OpenWindow(0, 0, 0, 512, 512, "Random generator", #PB_Window_SystemMenu | #PB_Window_ScreenCentered)
  If CreateImage(0, 512, 512) And StartDrawing(ImageOutput(0))
    
    For i = 0  To 131072
      R.l = mwc256()
      X.l = R & 511
      Y.l = (R >> 16) & 511
      Plot(X, Y)
    Next
    
    StopDrawing() 
    ImageGadget(0, 0, 0, 200, 200, ImageID(0))
  EndIf
  
  Repeat
    Event = WaitWindowEvent()
  Until Event = #PB_Event_CloseWindow
EndIf

Update 2012/02/25 :
Added MWC256 generator

[IdeasVacuum]

Nice.

[wilbert]

Nice.

Thanks

I like MWC256 myself.
It seems to give good random numbers and is very fast.
When I compare mwc256range(0, 500) with Random(500) , mwc is about 30% faster on my computer.

[luis]

Thanks

I don't know about the speed (I take your word on it!) but I didn't see an important difference in the distribution of the values compared with the PB random().
I tried to generate from 10-30 numbers to some thousands and the distribution seems comparable between PB and your routines.

Is it true and the advantage is only in speed or am I missing something ?

[wilbert]

From what I have read, different algorithms have different advantages.
The PB generator seems to give good random values but I have no idea what algorithm it is using.

I think maybe the biggest advantage is not depending on the PB generator.
A little while ago, I was surprised the PB compression algorithm isn't cross platform compatible.
When you are integrating a random number generator in your code, you can be sure the generated sequence with a specific seed is the same on all platforms.

[Thorium]

When you are integrating a random number generator in your code, you can be sure the generated sequence with a specific seed is the same on all platforms.

Only if your code is cross platform, which it isnt.
Would be nice if you add code for x64.

[wilbert]

Only if your code is cross platform, which it isnt.
Would be nice if you add code for x64.

On my computer results are the same for both x86 and x64.
Try for yourself

Code: Select all

XIncludeFile "mwc256.pbi"

mwc256seed(0)

Debug mwc256range(0, 999)
Debug mwc256range(0, 999)
Debug mwc256range(0, 999)

End

It's true I didn't use a lot of 64 bit registers.
That's because the algorithm works mainly with 32 bit values.
If for example I would use RAX for the multiplication instead of EAX, I would have to shift the result to get the high dword while it is already available directly as EDX when I use 32 bit multiplication.
I tried it but it would make the code longer and on my computer (Core2Duo) it wasn't faster.
In case you want to try for yourself, here's the C code

/* MWC256 from Usenet posting by G. Marsaglia - Period 2^8222 */ static unsigned int Q[256], c=362436;
unsigned int MWC256(void)
{
unsigned long long t;
static unsigned char i=255;
t = 809430660ULL * Q[++i] + c;
c = (t>>32);
return (Q=t);
}

[Thorium]

Ah, your right, sorry.