PRNG conversion from masm32 SDK

[zedd151]

I have converted the pseudo random number generator from the Masm32SDK to 64 bits.
Not for beginners. 

Original algorithm from the masm32 SDK:

Code Select Expand

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

nrandom PROC base:DWORD

    mov eax, nrandom_seed

  ; ****************************************
    test eax, 80000000h
    jz  @F
    add eax, 7fffffffh
  @@: 
  ; ****************************************

    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
Coverted version:

Code Select Expand

; «««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««
;
;                    Park Miller random number algorithm.
;
;                      Written by Jaymeson Trudgen (NaN)
;                  Optimized by Rickey Bowers Jr. (bitRAKE)
;                      converted to 64 bit - zedd151 :D
; «««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««

rand64 proc base:qword
.data
randseed dq 0
.code
  .if randseed == 0
    invoke GetTickCount
    mov randseed, rax
  .endif
    mov rax, randseed
    xor rdx, rdx
    mov rcx, 127773
    div rcx
    mov rcx, rax
    mov rax, 16807
    mul rdx
    mov rdx, rcx
    mov rcx, rax
    mov rax, 2836
    mul rdx
    sub rcx, rax
    xor rdx, rdx
    mov rax, rcx
    mov randseed, rcx
    div base
    mov rax, rdx
    ret
rand64 endp

Testing the randomness using the "ent.exe" program used by hutch-- in the attachment in This Thread.
The testbed for this algorithm:

Code Select Expand

include \masm64\include64\masm64rt.inc

rand64      proto :qword

reps equ 40000000h  ;; 1 GB

.code

start proc
local houtput:qword, bw:qword, memory:qword

    fn CreateFile, "random.bin", GENERIC_WRITE, FILE_SHARE_WRITE, 0, CREATE_ALWAYS, 0, 0
    mov houtput, rax
    invoke GlobalAlloc, GPTR, reps
    mov memory, rax
    mov r10, 0
    mov rsi, memory
  top:
    invoke rand64, 256
    mov byte ptr [rsi+r10], al
    inc r10
    cmp r10, reps
    jnz top
    invoke WriteFile, houtput, memory, reps, addr bw, 0
    invoke CloseHandle, houtput
    invoke GlobalFree, memory



   
  quit:
    invoke ExitProcess, 0
start endp

; «««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««
;
;                    Park Miller random number algorithm.
;
;                      Written by Jaymeson Trudgen (NaN)
;                  Optimized by Rickey Bowers Jr. (bitRAKE)
;                      converted to 64 bit - zedd151 :D
; «««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««

rand64 proc base:qword
.data
randseed dq 0
.code
  .if randseed == 0
    invoke GetTickCount
    mov randseed, rax
  .endif
    mov rax, randseed
    xor rdx, rdx
    mov rcx, 127773
    div rcx
    mov rcx, rax
    mov rax, 16807
    mul rdx
    mov rdx, rcx
    mov rcx, rax
    mov rax, 2836
    mul rdx
    sub rcx, rax
    xor rdx, rdx
    mov rax, rcx
    mov randseed, rcx
    div base
    mov rax, rdx
    ret
rand64 endp

end
Entire project is attached below. Using the batch file "test_the _prng.bat" it will first run the testbed to produce a 1 GB file filled with random byte data.
Then, the program "ent.exe" is executed to evaluate the file "random.bin" for randomness.
Typical result:

Entropy = 8.000000 bits per byte.

Optimum compression would reduce the size
of this 1073741824 byte file by 0 percent.

Chi square distribution for 1073741824 samples is 144.83, and randomly
would exceed this value more than than 99.99 percent of the times.

Arithmetic mean value of data bytes is 127.5015 (127.5 = random).
Monte Carlo value for Pi is 3.141702612 (error 0.00 percent).
Serial correlation coefficient is -0.000031 (totally uncorrelated = 0.0).

* Updated version rand64a.zip attached
Has issues with anything larger than dword sized base and output. see this post.
Disclaimer: Use at your own risk!!! Not intended for serous work, this is just experimental.

Originally I had thought that hutch himself wrote 'ent.exe', but later found out that it was downloaded from Github. Correction made above.

[sinsi]

Just curious about a few things
 Why change from using ECX to using RSI?
 Why the push/pop RDX?
 Why omit the test eax, 80000000h part? Admittedly, I don't understand why it's in the original.  Overflow?
 Why convert this to 64-bit? This is one algorithm that doesn't seem to benefit from 64-bitness

(this sounds way more confrontational than I mean it to be  )

[zedd151]

Hi sinsi, I had been converting some of my sudoku code to 64 bit, and needed a PRNG. The 32 bit code used the algorithm from masm32 SDK, so I decided to do this conversion.
The push/pop is because some of the code that uses the PRNG uses rdx and depends on it being the same after the call. For everything else, just because I can.    Not a great answer but it's all I got. Much of this was experimental. And I was thrilled that it all worked out for me, so I am sharing it here. Any helpful suggestions for changes are welcome. 

[sinsi]

Fair enough

[zedd151]

To be fair, I had only ran tests on byte values. The algorithm may not be suitable for larger than dword sized values. Further testing will be needed. I am thinking that the 'magic numbers' would have to be changed to work in 64 bit land. But that is beyond my skill set and pay grade.   

[zedd151]

To be fair, I had only ran tests on byte values. The algorithm may not be suitable for larger than dword sized values. Further testing will be needed. I am thinking that the 'magic numbers' would have to be changed to work in 64 bit land. But that is beyond my skill set and pay grade. 

Definitely not suitable for values larger than dword size.

results from ent.exe -- random qword sized data using the converted algo.

Entropy = 5.925595 bits per byte.

Optimum compression would reduce the size
of this 1073741824 byte file by 25 percent.

Chi square distribution for 1073741824 samples is 16860734851.80, and randomly
would exceed this value less than 0.01 percent of the times.

Arithmetic mean value of data bytes is 147.1846 (127.5 = random).
Monte Carlo value for Pi is 2.931844566 (error 6.68 percent).
Serial correlation coefficient is -0.358022 (totally uncorrelated = 0.0).

Generating dword size values does not look all that great either, but it is better nonetheless than generating qwords.
results from ent.exe -- random dword sized data using the converted algo.

Entropy = 7.955143 bits per byte.

Optimum compression would reduce the size
of this 1073741824 byte file by 0 percent.

Chi square distribution for 1073741824 samples is 66066933.33, and randomly
would exceed this value less than 0.01 percent of the times.

Arithmetic mean value of data bytes is 129.9778 (127.5 = random).
Monte Carlo value for Pi is 3.218921442 (error 2.46 percent).
Serial correlation coefficient is -0.002509 (totally uncorrelated = 0.0).

results from ent.exe -- random word sized data using the converted algo. for completeness.

Entropy = 8.000000 bits per byte.

Optimum compression would reduce the size
of this 1073741824 byte file by 0 percent.

Chi square distribution for 1073741824 samples is 185.61, and randomly
would exceed this value 99.96 percent of the times.

Arithmetic mean value of data bytes is 127.5008 (127.5 = random).
Monte Carlo value for Pi is 3.141651314 (error 0.00 percent).
Serial correlation coefficient is -0.000024 (totally uncorrelated = 0.0).

For byte sized data and word sized data,the new algo is totally acceptable.
dword sized data is not real bad, but could use improvement.
qword sized data, could be greatly improved upon. (Need new 'magic' numbers)

Given the discrepancies in the results with larger data sizes, I have asked that this thread be moved to The Workshop.    Not quite 'ready for primetime'. 

[mineiro]

Hi sir zedd151;

Code Select Expand

;ml64 zed151rand.asm
;link /subsystem:console zed151rand.obj /entry:start
;zed151rand

include \masm64\include64\masm64rt.inc
include \masm64\include64\msvcrt.inc

rand64      proto :qword
entropy    proto :qword,:ptr

reps equ 40000  ;; 1 GB

SYMBOLS EQU 256
.data
pfsymbols dq SYMBOLS dup (0)

.code

start proc
local houtput:qword, bw:qword, memory:qword, pfrequency:qword

    fn CreateFile, "random.bin", GENERIC_WRITE, FILE_SHARE_WRITE, 0, CREATE_ALWAYS, 0, 0
    mov houtput, rax
    invoke GlobalAlloc, GPTR, reps
    mov memory, rax
    mov r10, 0
    mov rsi, memory
    lea r11,pfsymbols
  top:
    invoke rand64, 256
    mov byte ptr [rsi+r10], al

;---------------------------------------------------------
    and rax,0ffh        ;get only a byte
    inc qword ptr [r11+rax*sizeof qword]    ;increment symbol frequency
;---------------------------------------------------------

    inc r10
    cmp r10, reps
    jnz top
    invoke WriteFile, houtput, memory, reps, addr bw, 0
    invoke CloseHandle, houtput
    invoke GlobalFree, memory

;---------------------------------------------------------
    invoke entropy,SYMBOLS,addr pfsymbols
    movsd finalent,xmm0

.data

    fmt_msg db "entropy = %.3f bits",0dh,0ah,0
    fmt_msg1 db "entropy = %.3f bytes",0dh,0ah,0
    align 16
    finalent real8 0.0
   
.code
    invoke vc__cprintf,addr fmt_msg,finalent

    movsd xmm1,finalent
    mov rax,8               ;per byte (symbol)
    cvtsi2sd xmm0,rax
    divsd xmm1,xmm0
    movsd finalent,xmm1
    invoke vc__cprintf,addr fmt_msg1,finalent
;---------------------------------------------------------

  quit:
    invoke ExitProcess, 0
start endp

; «««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««
;
;                    Park Miller random number algorithm.
;
;                      Written by Jaymeson Trudgen (NaN)
;                  Optimized by Rickey Bowers Jr. (bitRAKE)
;                      converted to 64 bit - zedd151 :D
; «««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««

rand64 proc base:qword
.data
randseed dq 0
.code
  .if randseed == 0
    invoke GetTickCount
    mov randseed, rax
  .endif
    push rsi
    push rdx
    mov rax, randseed
    xor rdx, rdx
    mov rsi, 127773
    div rsi
    mov rsi, rax
    mov rax, 16807
    mul rdx
    mov rdx, rsi
    mov rsi, rax
    mov rax, 2836
    mul rdx
    sub rsi, rax
    xor rdx, rdx
    mov rax, rsi
    mov randseed, rsi
    div base
    mov rax, rdx
    pop rdx
    pop rsi
    ret
rand64 endp






align 16
;return in xmm0 entropy of symbols frequency table
;entropy measure in bits, so if you want result in bytes should divide by 8 return value (xmm0)
;input:
;_symbolscounter == 256 for bytes, ...
;_psymbolsfrequency == pointer to symbols frequency table of size qword each

;THIS FUNCTION OVERWRITE VALUES (FREQUENCY) IN ORIGINAL FREQUENCY TABLE
entropy proc uses r12 r14 r15 _symbolscounter:qword,_psymbolsfrequency:ptr
    local symbol_counter:qword
    local symbol_frequency_table:ptr
    local sum:real8
   
    mov symbol_counter,rcx
    mov symbol_frequency_table,rdx
   
    mov sum,0
   
    xor r12,r12
    mov r14,symbol_frequency_table
    xor r15,r15

@@:
    cmp r15,symbol_counter            ;.while r15 != symbol_counter
    je @F
        mov rax,qword ptr [r14]       ;symbol frequency
        add r12,rax                   ;accumulator

        add r14,sizeof qword          ;next symbol
        inc r15
        jmp @B
@@:
    mov sum,r12                        ;symbols frequency sum


    mov r12,0
    mov r14,symbol_frequency_table
    mov rax,sum
    cvtsi2sd xmm1,rax                  ;sum.0
    finit
    pxor xmm2,xmm2                     ;0.0

@@:
    cmp r12,symbol_counter                 ;.while r12 != symbol_counter
    je @F
   
    cmp qword ptr [r14+r12*sizeof qword],0     ;don't count symbols with zero frequency
    je end_if

;(symbol frequency / sum frequency) * log2(symbol frequency / sum frequency)
;        .if qword ptr [r14+r12*8] != 0
            mov rax,qword ptr [r14+r12*sizeof qword]     ;get unsigned integer symbol frequency
            cvtsi2sd xmm0,rax                            ;convert to real8
            divsd xmm0,xmm1                              ;symbol frequency / frequency sum
            movsd real8 ptr [r14+r12*sizeof qword],xmm0  ;store in frequency table
           
            fld qword ptr [r14+r12*sizeof qword]         ;calculate x*log2(x)
            fld qword ptr [r14+r12*sizeof qword]
            fyl2x                        ;x * log2 (x)
            fstp qword ptr [r14+r12*sizeof qword]        ;store in frequency table
           
            addsd xmm2,qword ptr [r14+r12*sizeof qword]  ;partial entropy sum
;        .endif
end_if:
        inc r12           ;next symbol frequency
jmp @B
@@:                                 ;.endw

;final entropy sum * frequency sum
    mulsd xmm1,xmm2
    movsd sum,xmm1

    movsd xmm0,sum

ret
entropy endp

end


reps equ 40000 bytes

Code Select Expand

C:\masm64>zed151rand
entropy = -319808.597 bits
entropy = -39976.075 bytes


[zedd151]

Hola mineiro, how many bits per byte entropy? Judging by your output, seems very close to 8 bits per byte, which is ideal.  I don't know where you are getting the '-' sign from, but that looks wrong.

~320000 bits (rounded) divided by
~40000 bytes  (rounded)
= ~ 8 bits per byte. Where '~' means "around". I didn't make proper calculation, as it is almost impossible on my iPad. But that is close enough, I guess.

Why are you only testing with 40000 bytes? A larger data set should be used, imo for better accuracy.

[zedd151]

Why omit the test eax, 80000000h part? Admittedly, I don't understand why it's in the original.

I have been doing a little research on the park-miller implementation. I believe that the original was only used for the range 0-7FFFFFFFh. Hence why using dword values (full range 0-FFFFFFFFh) is less than optimal. Forget about using it for qwords
For my purposes, I only needed it for 0-46656D range. So it is very adequate.

I also looked up where the 'magic' numbers came from. Will be looking for something (magic numbers) that would make the algo suitable for qword values. At some time in the distant future. 

[mineiro]

Hi sir zedd151;

reps equ 40000 bytes
Code Select
C:\masm64>zed151rand
entropy = -319808.597 bits
entropy = -39976.075 bytes

319808.597 / 40000 == 7,995214925 bits per symbol.

It's not necessary to write that file in disk if we have a frequency table. Ent.exe give us more information, only in this situation. Now you can try that algo with 1 terabyte random digits and check entropy results.
It's also possible to try that with 2 symbols (0 and 1), with nibbles (4 bits), 65536 symbols ... . Random digits generally create a frequency table where each symbol frequency is very close to others symbols frequency, so we can't predict what's next symbol, generating poor compression (result in bytes).

The result in bytes it's just for convenience, so you can check what is supposed to be the "compressed" file size.
An use example is: open a txt file (like an .inc file), create a frequency table of all symbols (chars(bytes)) inside that file, apply entropy in that frequency table and the result in bytes tell's you the compressed file size. The result is "order 0" model (we are just counting symbols frequency without any other correlation (like digraphs, ...).

The minus sign is right, the entropy result is negative, but you can convert that to positive, the result will be the same.

An example done by hands and calculator:
message = BANANA
Symbol frequency:
A = 3
B = 1
N = 2
SUM = 3+1+2 = 6

= (3/6)*log2(3/6) + (1/6)*log2(1/6) + (2/6)*log2(2/6)
= 0,5 *log2(0,5)  + (0,166666667)*log2(0,166666667) + (0,333333333)*log2(0,333333333)
= 0,5 * −1        + 0,166666667*−2,584962498        +  0,333333333*−1,584962502
= −0,5            + −0,430827084                    +  −0,528320834
= −1,459147918 bits
So, this is the medium result of bits per symbol, but we have 3 symbols that appear 6 times(SUM), we need multiply:
= −1,459147918 bits * 6 symbols == −8,754887502 bits
The result −8,754887502 bits means that we can send message BANANA using ~9 bits.

Let's check: If I "encode" symbols like A=0, symbol B=10 and symbol N=11 (without ambiguity) I can send the info in bits:
B    A   N   A  N   A
10   0   11  0  11  0 
So, the message BANANA was sent as 100110110 bits, counting bits we have exactly 9 bits.

[mineiro]

If you want result in bits per symbols, I do this simple modification:

Code Select Expand

;ml64 zed151rand.asm
;link /subsystem:console zed151rand.obj /entry:start
;zed151rand

include \masm64\include64\masm64rt.inc
include \masm64\include64\msvcrt.inc

rand64      proto :qword
entropy    proto :qword,:ptr

reps equ 40000  ;; 1 GB

SYMBOLS EQU 256
.data
pfsymbols dq SYMBOLS dup (0)



.code

start proc
local houtput:qword, bw:qword, memory:qword, pfrequency:qword


    fn CreateFile, "random.bin", GENERIC_WRITE, FILE_SHARE_WRITE, 0, CREATE_ALWAYS, 0, 0
    mov houtput, rax
    invoke GlobalAlloc, GPTR, reps
    mov memory, rax
    mov r10, 0
    mov rsi, memory
    lea r11,pfsymbols
  top:
    invoke rand64, 256
    mov byte ptr [rsi+r10], al

;---------------------------------------------------------
    and rax,0ffh        ;get only a byte
    inc qword ptr [r11+rax*sizeof qword]    ;increment symbol frequency
;---------------------------------------------------------


    inc r10
    cmp r10, reps
    jnz top
    invoke WriteFile, houtput, memory, reps, addr bw, 0
    invoke CloseHandle, houtput
    invoke GlobalFree, memory


;---------------------------------------------------------
    invoke entropy,SYMBOLS,addr pfsymbols
    movsd finalent,xmm0

.data

    fmt_msg db "entropy = %.3f bits per symbols",0dh,0ah,0
    align 16
    finalent real8 0.0
   
.code
    invoke vc__cprintf,addr fmt_msg,finalent

  quit:
    invoke ExitProcess, 0
start endp

; «««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««
;
;                    Park Miller random number algorithm.
;
;                      Written by Jaymeson Trudgen (NaN)
;                  Optimized by Rickey Bowers Jr. (bitRAKE)
;                      converted to 64 bit - zedd151 :D
; «««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««

rand64 proc base:qword
.data
randseed dq 0
.code
  .if randseed == 0
    invoke GetTickCount
    mov randseed, rax
  .endif
    push rsi
    push rdx
    mov rax, randseed
    xor rdx, rdx
    mov rsi, 127773
    div rsi
    mov rsi, rax
    mov rax, 16807
    mul rdx
    mov rdx, rsi
    mov rsi, rax
    mov rax, 2836
    mul rdx
    sub rsi, rax
    xor rdx, rdx
    mov rax, rsi
    mov randseed, rsi
    div base
    mov rax, rdx
    pop rdx
    pop rsi
    ret
rand64 endp






align 16
;return in xmm0 entropy of symbols frequency table
;entropy measure in bits, so if you want result in bytes should divide by 8 return value (xmm0)
;input:
;_symbolscounter == 256 for bytes, ...
;_psymbolsfrequency == pointer to symbols frequency table of size qword each

;THIS FUNCTION OVERWRITE VALUES (FREQUENCY) IN ORIGINAL FREQUENCY TABLE
entropy proc uses r12 r14 r15 _symbolscounter:qword,_psymbolsfrequency:ptr
local symbol_counter:qword
local symbol_frequency_table:ptr
local sum:real8

mov symbol_counter,rcx
mov symbol_frequency_table,rdx

mov sum,0

xor r12,r12
mov r14,symbol_frequency_table
xor r15,r15

@@:
    cmp r15,symbol_counter            ;.while r15 != symbol_counter
    je @F
mov rax,qword ptr [r14]       ;symbol frequency
add r12,rax                   ;accumulator

add r14,sizeof qword          ;next symbol
inc r15
jmp @B
@@:
mov sum,r12                        ;symbols frequency sum


mov r12,0
mov r14,symbol_frequency_table
mov rax,sum
cvtsi2sd xmm1,rax                  ;sum.0
finit
pxor xmm2,xmm2                     ;0.0

@@:
cmp r12,symbol_counter                 ;.while r12 != symbol_counter
je @F

cmp qword ptr [r14+r12*sizeof qword],0     ;don't count symbols with zero frequency
je end_if

;(symbol frequency / sum frequency) * log2(symbol frequency / sum frequency)
; .if qword ptr [r14+r12*8] != 0
mov rax,qword ptr [r14+r12*sizeof qword]     ;get unsigned integer symbol frequency
cvtsi2sd xmm0,rax                            ;convert to real8
divsd xmm0,xmm1                              ;symbol frequency / frequency sum
movsd real8 ptr [r14+r12*sizeof qword],xmm0  ;store in frequency table

fld qword ptr [r14+r12*sizeof qword]         ;calculate x*log2(x)
fld qword ptr [r14+r12*sizeof qword]
fyl2x                        ;x * log2 (x)
fstp qword ptr [r14+r12*sizeof qword]        ;store in frequency table

addsd xmm2,qword ptr [r14+r12*sizeof qword]  ;partial entropy sum
; .endif
end_if:
inc r12           ;next symbol frequency
jmp @B
@@:                                 ;.endw

;------------------------------------------
;final entropy sum * frequency sum
; mulsd xmm1,xmm2
; movsd sum,xmm1
;------------------------------------------
    movsd sum,xmm2
;------------------------------------------

movsd xmm0,sum

ret
entropy endp

end

medium of bits per symbols

Code Select Expand

C:\masm64>zed151rand
entropy = -7.995 bits per symbols

It's also possible to modify code to give bits per symbol instead of medium result of bits per symbolS.

---------edited--------
check link with C source code of ent.exe https://github.com/ProhtMeyhet/ent-random-entropy-test

[zedd151]

Okay.   
At any rate, usually I seldom have a need for 'random' numbers exceeding 1000000, which the algo does well for, statically speaking.

Also, I recommended using the 'ent' program so that the results can more easily be compared with what I had posted, or what others may post.

Just as in the Laboratory, when algorithms are tested for speed, it is recommended that all participants use the exact same testing methods. Very similar reasoning. 

[zedd151]

Just curious about a few things
 Why change from using ECX to using RSI?

Registers now restored to (almost) original design. Version rand64a.zip in first post.

[zedd151]

Why omit the test eax, 80000000h part? Admittedly, I don't understand why it's in the original.  Overflow?

I have been scouring through some papers regarding the original Park-Miller algorithm. Apparently the seed value should be an integer less than or equal to 2147483647  (2^31)-1.   
I have been looking at some other prng algos, but specifically for 64 bit use. The Park-Miller algo is useful, but its not for 64 bit. Anyway I will be experimenting with some new 'magic' numbers from some of the 64 bit algos I have been reading about, attempting to adapt the Park-Miller algo further.