[Source] x86 Random Generator: Mersenne Twister

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.

Citation:

---

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:

---

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 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464

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 s:DWORD   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_key:DWORD, key_length:DWORD   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 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 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 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: 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

---

Include:

Code:

---

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:DWORD init_by_array PROTO init_key:DWORD,key_length:DWORD   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:

---

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:PTR 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:DWORD 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:DWORD 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:DWORD 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:DWORD   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:DWORD   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:DWORD   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:DWORD 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() ):

Citation:

---

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)
...

---