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 http://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 として実行すると一瞬で画面がランダムな色の点で埋めつくされます.