Bonjour,

Ci-joint une implémentation en assembleur x86 compatible MASM de l'algorithme MT19937 (version 2002) aussi connu sous le nom de Mersenne Twister.

L'implémentation a été testée sur environ 100 000 itérations et donne les résultats attendus. Par contre elle ne fonctionne probablement pas sur x64 (pas testé). Compilation possible avec MASM, ML (Visual Studio) ou JWasm.

La rapidité de cette implémentation est sensiblement égale à celle de référence (en C). Pour des raisons pratique et de lecture du source je n'ai pas "inliné" les appels à la fonction gen_randint32, même si cela aurait pu éviter quelques CALLs.

A noter: il existe une implémentation en C utilisant les jeux SIMD (dite SFMT) bien plus rapide, mais pas 100% assembleur.

Un programme de test pour Windows est fourni. Le code est distribué sous licence BSD.


assemble: ml /c /coff /Cp /nologo "MersenneTwisterTest.asm" "MersenneTwister.asm"
link: link /subsystem:console MersenneTwisterTest.obj MersenneTwister.obj kernel32.lib msvcrt.lib
Mersenne Twister:

[code=asm]
TITLE MersenneTwister.asm

comment @
A x86 MASM compatible program for MT19937, with initialization improved 2002/2/10.
Coded by Neitsa, based on implementation by Takuji Nishimura and Makoto Matsumoto.
This is a faster version by taking Shawn Cokus's optimization,
Matthe Bellew's simplification, Isaku Wada's real version.

Before using, initialize the state by using init_genrand(seed)
or init_by_array(init_key, key_length).

Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
are met:

1. Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.

3. The names of its contributors may not be used to endorse or promote
products derived from this software without specific prior written
permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
@

.686P
.XMM
.model flat, stdcall ;32 bit memory model
option casemap :none ;case sensitive


.const
N equ 624 ; degree of recurrence: number of 32-bit integers in the internal state array.
M equ 397 ; middle word, or the number of parallel sequences, 1 <= m <= n.
MATRIX_A equ 09908b0dfH ; constant vector a: coefficients of the rational normal form twist matrix.
UMASK equ 080000000H ; most significant w-r bits
LMASK equ 07fffffffH ; least significant r bits

.data
TWOPOW32 REAL8 4294967296.0 ; double(2^32)
ONEDIV2POW32 REAL8 03df0000000000000r ; 1.0 / double(2^32)
ONEDIV2POW32M1 REAL8 03df0000000100000r ; 1.0 / (double(2^32) - 1.0)
ONEDIV2 REAL8 0.5 ; 1.0 / 2.0
D1 REAL8 67108864.0 ; used by gen_randres53
D2 REAL8 03ca0000000000000r ; used by gen_randres53

.data?
_state DD N DUP (?); internal: random generator state.
_initf DD ? ; set if the internal state has been initalized.
_left DD ? ; number of generation left before a new internal state is required.
_next DD ? ; holds pointer to the next internal state.


.code
;============================================================
; Function: init_genrand
; Initializes state[N] with a seed.
; ARGS:
; s: The seed. Must be a 32-bit integer.
; Output: None
; Precondition: None
;============================================================
init_genrand PROC uses EDI sWORD

mov eax, s
mov _state, eax

mov edx, 1
lea eax, _state

ALIGN 4
@loop_init_genrand:
mov ecx, DWORD PTR [eax]
mov edi, ecx
shr edi, 30
xor edi, ecx
imul edi, 1812433253
add edi, edx
; opt: not needed if compiling on a 32-bit machine because we use machine word size
;and edi, 0FFFFFFFFH
mov DWORD PTR [eax+4], edi
add eax, 4
inc edx
cmp eax, OFFSET _state + ( (N * sizeof DWORD) - SIZEOF DWORD )
jl @loop_init_genrand

mov _left, 1
mov _initf, 1

ret

init_genrand endp

;============================================================
; Function: init_by_array
; Initializes the generator with an array.
; ARGS:
; init_key: An array of 32-bit integers. [no size limit]
; key_length: Number of 32-bit integers in init_key.
; Output: None
; Precondition: None
;============================================================
init_by_array PROC uses EBX ESI EDI, init_keyWORD, key_lengthWORD

push 19650218
call init_genrand

mov esi, init_key

mov eax, 1
xor ecx, ecx

mov edi, N
cmp edi, key_length
jb @set_key_length
jmp @loop1_init_by_array
@set_key_length:
mov edi, key_length

ALIGN 4
@loop1_init_by_array:
mov edx, DWORD PTR _state[eax*4 - SIZEOF DWORD]
mov ebx, edx
shr ebx, 30
xor ebx, edx
imul ebx, 1664525
xor ebx, DWORD PTR _state[eax*4]
add ebx, DWORD PTR [esi+ecx*4]
add ebx, ecx
; opt: not needed if compiling on a 32-bit machine because we already use machine word size
;and ebx, 0FFFFFFFFH
mov DWORD PTR _state[eax*4], ebx
inc eax
inc ecx

cmp eax, N
jl @i_below_n1
mov eax, DWORD PTR _state + ( (N * sizeof DWORD) - SIZEOF DWORD )
mov DWORD PTR _state, eax
mov eax, 1

@i_below_n1:
mov edx, key_length
cmp ecx, edx
jl @j_below_key_length
xor ecx, ecx

@j_below_key_length:
dec edi
jne @loop1_init_by_array
;----------


mov edx, N-1

ALIGN 4
@loop2_init_by_array:
mov ecx, DWORD PTR _state[eax*4-SIZEOF DWORD]
mov esi, ecx
shr esi, 30
xor esi, ecx
imul esi, 1566083941
xor esi, DWORD PTR _state[eax*4]
sub esi, eax
mov DWORD PTR _state[eax*4], esi
inc eax

cmp eax, N
jl @i_below_n2
mov ecx, DWORD PTR _state + ( (N * sizeof DWORD) - SIZEOF DWORD )
mov DWORD PTR _state, ecx
mov eax, 1

@i_below_n2:
dec edx
jne @loop2_init_by_array

mov DWORD PTR _state, 080000000h
mov DWORD PTR _left, 1
mov DWORD PTR _initf, 1

ret

init_by_array endp

;============================================================
; Function: next_state
; Generates the next internal state for the random generator.
; ARGS: None
; Output: None
; Precondition: None
;============================================================
next_state PROC uses ESI

lea esi, _state

cmp DWORD PTR _initf, 0
jne @dont_call_initgenrand
push 5489
call init_genrand

@dont_call_initgenrand:
mov _left, N
mov _next, esi

mov ecx, N-M+1-1

ALIGN 4
@next_state_loop1:

mov edx, DWORD PTR [esi]
mov eax, DWORD PTR [esi+4]

and edx, UMASK
and eax, LMASK
or edx, eax

shr edx, 1

and eax, 1
neg eax
sbb eax, eax
and eax, MATRIX_A
xor edx, eax
xor edx, DWORD PTR [esi + (M * SIZEOF DWORD)]
mov DWORD PTR [esi], edx
add esi, 4
dec ecx
jne @next_state_loop1

mov ecx, M-1

ALIGN 4
@next_state_loop2:
mov edx, DWORD PTR [esi]
mov eax, DWORD PTR [esi+4]
and edx, UMASK
and eax, LMASK
or edx, eax
shr edx, 1
and eax, 1
neg eax
sbb eax, eax
and eax, MATRIX_A
xor edx, eax
xor edx, DWORD PTR [esi + ( (M-N)*SIZEOF DWORD )] ; note: M-N will be negative, so that gives [ESI+(-908)] -> [ESI-908] with standard M and N
mov DWORD PTR [esi], edx
add esi, 4
dec ecx
jne @next_state_loop2


mov edx, DWORD PTR [esi]
mov eax, DWORD PTR _state
and edx, UMASK
and eax, LMASK
or edx, eax
shr edx, 1
and eax, 1
neg eax
sbb eax, eax
and eax, MATRIX_A
xor edx, eax
xor edx, DWORD PTR [esi+( (M-N)*SIZEOF DWORD )] ; see note above
mov DWORD PTR [esi], edx

ret

next_state endp

;============================================================
; Function: gen_randint32
; generates a 32-bit random number on [0,0xffffffff] interval.
; ARGS: None
; Output: an integer, the result is placed in the EAX register.
; Precondition: One of the initialization functions must have been called.
;============================================================
gen_randint32 PROC

dec DWORD PTR _left
jne @L1_gen_randint32
call next_state

@L1_gen_randint32:


mov eax, DWORD PTR _next
mov ecx, DWORD PTR [eax]
add eax, 4
mov DWORD PTR _next, eax

mov eax, ecx
shr eax, 11
xor ecx, eax

mov edx, ecx
shl edx, 7
and edx, 09d2c5680H
xor ecx, edx

mov eax, ecx
shl eax, 15
and eax, 0efc60000H
xor ecx, eax

mov eax, ecx
shr eax, 18
xor eax, ecx

ret

gen_randint32 endp

;============================================================
; Function: gen_randint31
; generates a 32-bit random number on [0,0x7fffffff] interval.
; ARGS: None
; Output: an integer, the result is placed in the EAX register.
; Precondition: One of the initialization functions must have been called.
;============================================================
gen_randint31 PROC

call gen_randint32
shr eax, 1

ret

gen_randint31 endp

;============================================================
; Function: gen_real1
; Generates a random number on [0,1) real interval ; ARGS: None ; Output: a double precision (REAL8 / QWORD) number. ; On function exit, the output is placed on top ; of the FPU stack, st(0). ; Precondition: One of the initialization functions must have been called. ;============================================================ gen_real1 PROC LOCAL Temp:DWORD call gen_randint32 mov Temp, eax fild Temp test eax, eax jns @integer_not_signed fadd TWOPOW32 ; make sure it's >= 0 @integer_not_signed: fmul ONEDIV2POW32M1 ret gen_real1 endp ;============================================================ ; Function: gen_real2 ; Generates a random number on [0,1] real interval.
; ARGS: None
; Output: a double precision (REAL8 / QWORD) number.
; On function exit, the output is placed on top
; of the FPU stack, st(0).
; Precondition: One of the initialization functions must have been called.
;============================================================
gen_real2 PROC
LOCAL TempWORD

call gen_randint32
mov Temp, eax
fild Temp
test eax, eax
jns @integer_not_signed
fadd TWOPOW32 ; make sure it's >= 0
@integer_not_signed:
fmul ONEDIV2POW32

ret

gen_real2 endp

;============================================================
; Function: gen_real3
; Generates a random number on (0,1) real interval.
; ARGS: None
; Output: a double precision (REAL8 / QWORD) number.
; On function exit, the output is placed on top
; of the FPU stack, st(0).
; Precondition: One of the initialization functions must have been called.
;============================================================
gen_real3 PROC
LOCAL TempWORD

call gen_randint32
mov Temp, eax
fild Temp
test eax, eax
jns @integer_not_signed
fadd TWOPOW32 ; make sure it's >= 0
@integer_not_signed:
fadd ONEDIV2
fmul ONEDIV2POW32

ret

gen_real3 endp

;============================================================
; Function: gen_randres53
; Generates a random number on [0,1) with 53-bit resolution. ; ARGS: None ; Output: a double precision (REAL8 / QWORD) number. ; On function exit, the output is placed on top ; of the FPU stack, st(0). ; Precondition: One of the initialization functions must have been called. ; PostCondition: st(0) and st(1) are destroyed. ;============================================================ gen_randres53 PROC uses ESI LOCAL Temp1:DWORD LOCAL Temp2:DWORD call gen_randint32 mov esi, eax shr esi, 5 mov Temp2, esi call gen_randint32 shr eax, 6 mov Temp1, eax fild Temp1 test eax, eax jns @L1 fadd TWOPOW32 ; make sure it's >= 0 @L1: fild Temp2 test esi, esi jns @L2 fadd TWOPOW32 ; make sure it's >= 0 @L2: fmul D1 faddp st(1), st(0) fmul D2 ret gen_randres53 endp end [/code]

Include:

Code asm :Sélectionner tout -Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
 
;==================================================
;
; PROTOTYPES
;
;==================================================	
 
; initialization functions
init_genrand PROTO s<img src="images/smilies/icon_biggrin.gif" border="0" alt="" title=":D" class="inlineimg" />WORD
init_by_array PROTO init_key<img src="images/smilies/icon_biggrin.gif" border="0" alt="" title=":D" class="inlineimg" />WORD,key_length<img src="images/smilies/icon_biggrin.gif" border="0" alt="" title=":D" class="inlineimg" />WORD
 
PUBLIC init_genrand
PUBLIC init_by_array
 
; generating functions
gen_randint32 PROTO
gen_randint31 PROTO
gen_real1 PROTO
gen_real2 PROTO
gen_real3 PROTO
gen_randres53 PROTO
 
PUBLIC gen_randint32
PUBLIC gen_randint31
PUBLIC gen_real1
PUBLIC gen_real2
PUBLIC gen_real3
PUBLIC gen_randres53


Programme de test:

Code :Sélectionner tout -Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
 
    TITLE    MersenneTwisterTest.asm
 
    comment $ 
    assemble: ml /c /coff /Cp /nologo "MersenneTwisterTest.asm" "MersenneTwister.asm"
    link: link /subsystem:console MersenneTwisterTest.obj MersenneTwister.obj kernel32.lib msvcrt.lib
    $
 
    .686P
    .XMM
    .model flat, stdcall  ;32 bit memory model
    option casemap :none  ;case sensitive
    include MersenneTwister.inc
 
 
    ; printf    
    c_msvcrt typedef PROTO C :VARARG
    externdef _imp__printf<img src="images/smilies/icon_razz.gif" border="0" alt="" title=":P" class="inlineimg" />TR c_msvcrt
    crt_printf equ <_imp__printf>
    ; exitprocess
    ExitProcess   proto :dword
 
;==================================================    
.const
    ; number of test loops
    MAX_LOOP_NUMBER equ 1000
 
;==================================================    
.data
    ; output format for integers
    FMT_INT        BYTE "Generated Number [%d]: %d (%#08lx)", 0ah, 0
 
    ; output format for double (REAL8) numbers
    FMT_DOUBLE    BYTE "Generated Number [%d]: %10.8f", 0ah, 0
 
    ;initializer array used by init_by_array()
    InitKey     DWORD 123h, 234h, 345h, 456h
    InitKeyLength DWORD SIZEOF InitKey SHR 2 
 
;==================================================
;==================================================    
.code
start:
    call Main
    ret
Main PROC
 
    finit   ; ensure FPU is initialized correctly.
 
    ; note: use only ONE test at a time (each test contain the state initialization)
 
    call Test_initgenrand_randint32 ; -> TEST: ok
 
    ;call Test_initbyarray_randint32 ; -> TEST: ok
    ;call Test_initgenrand_randint31 ; -> TEST: ok
    ;call Test_initgenrand_real1 ;-> TEST: ok
    ;call Test_initgenrand_real2 ;-> TEST: ok
    ;call Test_initgenrand_real3 ;-> TEST: ok
    ;call Test_initgenrand_genrandres53 ; -> TEST: ok    
 
    push 0
    call ExitProcess
 
    ret    
Main endp
;==================================================
;
; Test with init_genrand() + gen_real1
;
;==================================================    
Test_initgenrand_real1 PROC
    LOCAL counter<img src="images/smilies/icon_biggrin.gif" border="0" alt="" title=":D" class="inlineimg" />WORD
    LOCAL double:REAL8
 
    push 987654
    call init_genrand
 
    mov counter, 0
 
@loop:    
    call gen_real1
    ;on exit from gen_real1() function the result is placed on top of the FPU stack: st(0).
    fstp [double]    
 
    invoke crt_printf, OFFSET FMT_DOUBLE, counter, double
 
    inc counter
    cmp counter, MAX_LOOP_NUMBER
    jne @loop    
 
    ret    
 
Test_initgenrand_real1 endp
 
;==================================================
;
; Test with init_genrand() + gen_real2
;
;==================================================    
Test_initgenrand_real2 PROC
    LOCAL counter<img src="images/smilies/icon_biggrin.gif" border="0" alt="" title=":D" class="inlineimg" />WORD
    LOCAL double:REAL8
 
    push 424242
    call init_genrand
 
    mov counter, 0
 
@loop:    
    call gen_real2    
    ;on exit from gen_real2() function the result is placed on top of the FPU stack: st(0).
    fstp [double]    
 
    invoke crt_printf, OFFSET FMT_DOUBLE, counter, double
 
    inc counter
    cmp counter, MAX_LOOP_NUMBER
    jne @loop    
 
    ret    
 
Test_initgenrand_real2 endp
 
;==================================================
;
; Test with init_genrand() + gen_real2
;
;==================================================    
Test_initgenrand_real3 PROC
    LOCAL counter<img src="images/smilies/icon_biggrin.gif" border="0" alt="" title=":D" class="inlineimg" />WORD
    LOCAL double:REAL8
 
    push 123456
    call init_genrand
 
    mov counter, 0
 
@loop:    
    call gen_real3    
    ;on exit from gen_real3() function the result is placed on top of the FPU stack: st(0).
    fstp [double]    
 
    invoke crt_printf, OFFSET FMT_DOUBLE, counter, double
 
    inc counter
    cmp counter, MAX_LOOP_NUMBER
    jne @loop    
 
    ret    
 
Test_initgenrand_real3 endp
 
;==================================================
;
; Test with init_by_array() + gen_randint32
;
;==================================================    
Test_initbyarray_randint32 PROC
    LOCAL counter<img src="images/smilies/icon_biggrin.gif" border="0" alt="" title=":D" class="inlineimg" />WORD
 
    push InitKeyLength
    push OFFSET InitKey
    call init_by_array    
 
    mov counter, 0
 
@loop:    
    call gen_randint32
 
    invoke crt_printf, OFFSET FMT_INT, counter, eax, eax    
    inc counter
    cmp counter, MAX_LOOP_NUMBER
    jne @loop    
 
    ret    
Test_initbyarray_randint32 endp
 
;==================================================
;
; Test with init_genrand() + gen_randint32
;
;==================================================
Test_initgenrand_randint32 PROC
    LOCAL counter<img src="images/smilies/icon_biggrin.gif" border="0" alt="" title=":D" class="inlineimg" />WORD
 
    push 58
    call init_genrand
 
    mov counter, 0
 
@loop:    
    call gen_randint32
 
    invoke crt_printf, OFFSET FMT_INT, counter, eax, eax
 
    inc counter
    cmp counter, MAX_LOOP_NUMBER
    jne @loop    
 
    ret
Test_initgenrand_randint32 endp
 
;==================================================
;
; Test with init_genrand() + gen_randint31
;
;==================================================
Test_initgenrand_randint31 PROC
    LOCAL counter<img src="images/smilies/icon_biggrin.gif" border="0" alt="" title=":D" class="inlineimg" />WORD
 
    push 456321
    call init_genrand
 
    mov counter, 0
 
@loop:    
    call gen_randint31
 
    invoke crt_printf, OFFSET FMT_INT, counter, eax, eax
 
    inc counter
    cmp counter, MAX_LOOP_NUMBER
    jne @loop    
 
    ret
Test_initgenrand_randint31 endp
 
;==================================================
;
; Test with init_genrand() + gen_real2
;
;==================================================    
Test_initgenrand_genrandres53 PROC
    LOCAL counter<img src="images/smilies/icon_biggrin.gif" border="0" alt="" title=":D" class="inlineimg" />WORD
    LOCAL double:REAL8
 
    push 357951
    call init_genrand
 
    mov counter, 0
 
@loop:    
    call gen_randres53    
    ;on exit from gen_randres53() function the result is placed on top of the FPU stack: st(0).
    fstp [double]    
 
    invoke crt_printf, OFFSET FMT_DOUBLE, counter, double
 
    inc counter
    cmp counter, MAX_LOOP_NUMBER
    jne @loop    
 
    ret    
 
Test_initgenrand_genrandres53 endp
 
end start


Example d'output du programme de test ( initgenrand(58) + randint32() ):


Generated Number [0]: 1568116515 (0x5d778f23)
Generated Number [1]: -715449824 (0xd55b1a20)
Generated Number [2]: 1937914647 (0x73823b17)
Generated Number [3]: -466592743 (0xe4305c19)
Generated Number [4]: 2130562958 (0x7efdcf8e)
Generated Number [5]: -161059947 (0xf6666b95)
Generated Number [6]: 324794275 (0x135bf7a3)
Generated Number [7]: -150049036 (0xf70e6ef4)
Generated Number [8]: -1839269991 (0x925ef799)
Generated Number [9]: -1985513311 (0x89a778a1)
Generated Number [10]: 1243968602 (0x4a25745a)
Generated Number [11]: 886994981 (0x34de7825)
...