15. 乱数生成 -- Mersenne Twister --
前回のサンプルでは，/dev/urandom から取得した乱数を 使ってフレームバッファにランダムに点を描画しましたが，非常に遅いことに気がつきます． 今回は疑似乱数生成アルゴリズムとして Mersenne Twister の C による32ビット整数版 mt19937int.c をアセンブリに移植します．Mersenne Twisterは 周期が非常に長く (2^19937-1)，高次元(623次元)で均等に分布し，生成速度が速いという特徴をもつ 非常に高性能な疑似乱数生成アルゴリズムです．
mt19937int.c をコンパイラに翻訳させてみます．
$ gcc -O2 -o mt19937int mt19937int.c $ ls -l mt19937int -rwxr-xr-x 1 jun users 13939 Mar 11 16:08 mt19937int $ ndisasm -b 32 mt19937int >mt19937int.lst
逆アセンブルしたリスト (mt19937int.lst) を読むと(void lsgenrand(seed_array) を除いて) 400バイト強のコードが生成されています．このままでは面白くないので 小さいコードに挑戦してみました． 速度が若干を犠牲になりましたが，256 バイトまで小さくなりました．コードは読みにくく なっています．興味のある方は C のオリジナル mt19937int.c と読み比べてください．
;---------------------------------------------------------------------
; Mersenne Twister
; file : mt19937.inc
; Rewritten in Assembly by Jun Mizutani 2001/03/11.
; From original code in C by Takuji Nishimura(mt19937int.c).
;---------------------------------------------------------------------
; This library is free software; you can redistribute it and/or
; modify it under the terms of the GNU Library General Public
; License as published by the Free Software Foundation; either
; version 2 of the License, or (at your option) any later
; version.
; This library is distributed in the hope that it will be useful,
; but WITHOUT ANY WARRANTY; without even the implied warranty of
; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
; See the GNU Library General Public License for more details.
; You should have received a copy of the GNU Library General
; Public License along with this library; if not, write to the
; Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
; 02111-1307 USA
;-- The original C code contained the following notice:
; Copyright (C) 1997, 1999 Makoto Matsumoto and Takuji Nishimura.
; Any feedback is very welcome. For any question, comments,
; see https://www.math.keio.ac.jp/matumoto/emt.html or email
; matumoto@math.keio.ac.jp
%assign N 624
%assign M 397
section .text
;---------------------------------------------------------------------
; Initialize Mersenne Twister
; enter eax : seed
;---------------------------------------------------------------------
sgenrand:
pusha
mov esi, N
mov ebx, 69069
mov ebp, 0xffff0000
xor edx, edx
.loop:
mov edi, eax ; eax : seed
and edi, ebp ; 0xffff0000
imul eax, ebx ; seed * 69069
inc eax
mov ecx, eax
and ecx, ebp ; 0xffff0000
shr ecx, 16
or ecx, edi
mov [mt+edx*4], ecx ; mt[ebx]
imul eax, ebx ; 69069
inc eax ; seed
inc edx
cmp edx, esi
jl .loop
mov dword [mti], esi ; N:624
popa
ret
;---------------------------------------------------------------------
; Generate Random Number
; return eax : random number
;---------------------------------------------------------------------
genrand:
pusha
mov esi, mt
mov edi, mti
mov eax, [edi]
mov ecx, N-1
cmp eax, ecx ; 623
jle .genrand2
mov ebp, mag01
shl ecx, 2 ; (N-1)*4
xor ebx, ebx
.loop1: call .common
xor eax, [esi+ebx+M*4] ; mt[kk+M]
call .common2
cmp ebx, (N-M-1)*4 ; N-M-1
jle .loop1
cmp ebx, ecx ; (N-1)*4
jge .next
.loop2: call .common
xor eax, [esi+ebx+(M-N)*4] ; (M-N=-227)
call .common2
cmp ebx, ecx ; (N-1)*4
jl .loop2
.next:
mov edx, [esi+ecx] ; (N-1)*4
mov eax, [esi]
call .common1
xor eax, [esi+(M-1)*4]
call .common2 ; ebx = ecx
xor eax, eax
mov [edi], eax ; mti=0
.genrand2: mov eax, [edi]
mov edx, [esi+eax*4] ; mt[mti]
inc dword [edi] ; mti++
mov eax, edx
shr eax, 11
xor edx, eax
mov eax, edx
shl eax, 7
and eax, 0x9d2c5680 ; TEMPERING_MASK_B
xor edx, eax
mov eax, edx
shl eax, 15
and eax, 0xefc60000 ; TEMPERING_MASK_C
xor edx, eax
mov eax, edx
shr eax, 18
xor edx, eax
mov [esp+28], edx ; return eax
popa
ret
.common: mov edx, [esi+ebx]
mov eax, [esi+ebx+4]
.common1: and edx, 0x80000000
and eax, 0x7fffffff
or edx, eax
mov eax, edx
shr eax, 1
ret
.common2: and edx, byte 1
xor eax, [ebp+edx*4]
mov [esi+ebx], eax
add ebx, byte 4
ret
;==============================================================
section .data
mag01 dd 0x00000000
dd 0x9908b0df
mti dd N+1
;==============================================================
section .bss
mt resd N
実際に乱数を生成して表示してみます．コードが小さくなっても間違った結果 ではシャレになりませんからオリジナルの出力結果 mt19937int.out と同じフォーマットで出力しています．
;---------------------------------------------------------------------
; Mersenne Twister
; 2001/03/11 Jun Mizutani
; testmt.asm
;---------------------------------------------------------------------
%include "stdio.inc"
%include "mt19937.inc"
;==============================================================
section .text
global _start
_start:
mov eax, 0x1105 ; 4357(seed)
call sgenrand
mov ecx, 1000
xor ebx, ebx
.loop
call genrand
push ecx
mov ecx, 10
call PrintRightU
pop ecx
mov eax, ' '
call OutChar
inc ebx
cmp ebx, 5
jne .skip
xor ebx, ebx
call NewLine
.skip:
loop .loop
call Exit
初めに eax にシード(ここでは 4357，0は不可)を設定して call sgenrand として準備します． 後は call genrand で eax に 32ビットの乱数が得られます．eax 以外のレジスタは変化しません． 実行結果の一部を示します．
$ ./testmt 2867219139 1585203162 3113124129 2953900839 2463794868 3482265796 1164297043 3598195569 589972756 4112233867 767115311 4093075447 1322433849 3357085324 3300048468 3649464345 3676604632 1475054104 2601934239 3420804864 2492391180 28597038 1901037238 1209433535 3580317774 2488297452 79873538 3308484072 2913896343 4166196021 1930853421 3313543893 2603730014 2827553081 1952080899 1405101208 1959413290 2221997165 4110132150 1025637693
1000個の乱数についてオリジナルの出力結果 mt19937int.out と比較してみます．
$ ./testmt >testmt.out $ diff mt19937int.out testmt.out
何も表示されなければ一致しています． で，一致しました(^^;
前回のサンプルを，移植した mt19937.inc を使って 書き直します．1,000 個の点の表示では速すぎるので 1万倍して 10,000,000 個を 表示しています．また点の色も乱数で決めているため，前回のサンプルの 2万倍の 数の乱数を求めています．
;---------------------------------------------------------------------
; Frame buffer 2
; 2001/03/11 Jun Mizutani
; fbtest2.asm
;---------------------------------------------------------------------
%include "stdio.inc"
%include "fblib.inc"
%include "mt19937.inc"
%assign O_RDONLY 00q
;==============================================================
section .text
global _start
_start:
;-------------------------------------------
; フレームバッファの準備
;-------------------------------------------
call fbdev_open ; /dev/fb0 をオープン
js near .fb_Error
call fb_get_screen ; スクリーン情報を取得
js .fb_Error
call fb_copy_scinfo ; scinfo_data を設定
call fb_map_screen ; メモリにマッピング
js .fb_Error
mov [fb_mem], eax ; 先頭アドレス
mov edi, [fb_mem]
mov ebx, [mmap_arg.len] ; サイズ
;-------------------------------------------
; [edi] から [edi+ebx-1] の範囲に書き込み
;-------------------------------------------
call sgenrand
shr ebx, 1
mov ecx, 10000000
.loop: call genrand
xor edx, edx
div ebx
call genrand
mov [edi+edx*2], ax ; ランダムに点を打つ
loop .loop
xor esi, esi
mov ecx, ebx
shr ecx, 1
.loop2: mov eax, [edi+esi*4] ; 画面読み出し
xor eax, 0x55555555
mov [edi+esi*4], eax ; 画面書き込み
inc esi
loop .loop2
xor esi, esi
mov ecx, ebx
shr ecx, 1
.loop3: mov eax, [edi+esi*4] ; 画面読み出し
xor eax, 0x55555555
mov [edi+esi*4], eax ; 画面書き込み
inc esi
loop .loop3
;-------------------------------------------
; フレームバッファの後始末
;-------------------------------------------
.exit call fb_unmap_screen
call fbdev_close
call Exit
;-------------------------------------------
.fb_Error:
mov eax, fberr_msg
call OutAsciiZ
jmp short .exit
fberr_msg db 'frabe buffer error!', 10, 0
;==============================================================
section .bss
fb_mem resd 1
./fbtest2 として実行すると一瞬で画面がランダムな色の点で埋めつくされます．