Page 1 of 1

Random number generators ( MWC256 & CMWC4096 )

Posted: Wed Feb 22, 2012 4:12 pm
by 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

Re: Random number generators ( MWC256 & CMWC4096 )

Posted: Sat Feb 25, 2012 6:39 pm
by 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

Re: Random number generators ( MWC256 & CMWC4096 )

Posted: Sat Feb 25, 2012 6:41 pm
by 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

Re: Random number generators ( MWC256 & CMWC4096 )

Posted: Sat Feb 25, 2012 9:57 pm
by IdeasVacuum
Nice. :)

Re: Random number generators ( MWC256 & CMWC4096 )

Posted: Sun Feb 26, 2012 2:16 pm
by wilbert
IdeasVacuum wrote:Nice. :)
Thanks :D

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.

Re: Random number generators ( MWC256 & CMWC4096 )

Posted: Sun Feb 26, 2012 2:49 pm
by 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 ?

Re: Random number generators ( MWC256 & CMWC4096 )

Posted: Sun Feb 26, 2012 3:05 pm
by 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.

Re: Random number generators ( MWC256 & CMWC4096 )

Posted: Sun Feb 26, 2012 4:58 pm
by Thorium
wilbert wrote: 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.

Re: Random number generators ( MWC256 & CMWC4096 )

Posted: Sun Feb 26, 2012 5:12 pm
by wilbert
Thorium wrote: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);
}

Re: Random number generators ( MWC256 & CMWC4096 )

Posted: Sun Feb 26, 2012 6:36 pm
by Thorium
Ah, your right, sorry.