LuaでXorShiftかいてみたよ

関連：Luaでメルセンヌ・ツイスタ（MT乱数）かいてみたよ

唐突ですがLuaスクリプトはじめました。Luaスクリプトに対する個人的な感想はさておき、乱数生成の精度というか挙動がちょっと気に入らなかったので、知る人ぞ知るXorShiftという乱数生成アルゴリズムをLuaに移植してみました。学習初日から作りたいものも決まってないのになに書いてるんでしょうね (・д・)

移植元のCソースに忠実な動作をするよう頑張りましたが、穴があるかもしれません。ビット演算周りにLuaBitOpモジュールを使っています（Lua for Windowsの場合5.1.4.22より付属）。

2009-02-16 20時: random() の発生区間が [0,1) ではなく [0,1] になっている気がしたので修正。
2009-04-13 7時: 恥ずかしながら math.random() が引数を受け取ることを最近知りました。よって処理を追加。
xorshift.lua (pastebin link):

-- XorShift: A random number generator
-- ported to Lua by gocha

module("xorshift", package.seeall)

require "bit"

local seed128 = {}

function randomseed(s)
    for i = 0, 3 do
        -- s = 1812433253 * (bit.bxor(s, bit.rshift(s, 30))) + i
        s = bit.bxor(s, bit.rshift(s, 30))
        local s_lo = bit.band(s, 0xffff)
        local s_hi = bit.rshift(s, 16)
        local s_lo2 = bit.band(1812433253 * s_lo, 0xffffffff)
        local s_hi2 = bit.band(1812433253 * s_hi, 0xffff)
        s = bit.bor(bit.lshift(bit.rshift(s_lo2, 16) + s_hi2, 16),
            bit.band(s_lo2, 0xffff))
        -- s = bit.band(s + i, 0xffffffff)
        local s_lim = -bit.tobit(s)
        -- assumes i<2^31
        if (s_lim > 0 and s_lim <= i) then
            s = i - s_lim
        else
            s = s + i
        end
        seed128[i+1] = s
    end
end

function random_int32()
    local t = bit.bxor(seed128[1], bit.lshift(seed128[1], 11))
    seed128[1], seed128[2], seed128[3] = seed128[2], seed128[3], seed128[4]
    seed128[4] = bit.bxor(bit.bxor(seed128[4], bit.rshift(seed128[4], 19)), bit.bxor(t, bit.rshift(t, 8)))
    return seed128[4]
end

function random(...)
    -- local r = xorshift.random_int32() * (1.0/4294967296.0)
    local rtemp = xorshift.random_int32()
    local r = (bit.band(rtemp, 0x7fffffff) * (1.0/4294967296.0)) + (bit.tobit(rtemp) < 0 and 0.5 or 0)
    local arg = {...}
    if #arg == 0 then
        return r
    elseif #arg == 1 then
        local u = math.floor(arg[1])
        if 1 <= u then
            return math.floor(r*u)+1
        else
            error("bad argument #1 to 'random' (internal is empty)")
        end
    elseif #arg == 2 then
        local l, u = math.floor(arg[1]), math.floor(arg[2])
        if l <= u then
            return math.floor((r*(u-l+1))+l)
        else
            error("bad argument #2 to 'random' (internal is empty)")
        end
    else
        error("wrong number of arguments")
    end
end

math.randomseed(), math.random() と同じ感覚で xorshift.randomseed(), xorshift.random() が利用できます。xorshift.random_int32() は内部用に作ったつもりなので、外部からの利用はあんまりおすすめしません。でも、好きにしてください。とりあえず以下に適当なサンプルもおいておきます。

require "xorshift"

local randomseed = xorshift.randomseed
local random = xorshift.random

randomseed(os.time())
for i = 1, 10 do
    local rf = random()
    local ri = math.floor(rf * 10)
    print(ri, rf)
end

オリジナルよりもちょっとややこしい書き方で実装してありますね。random() では signed int32 で整数演算がおこなわれる環境でも unsigned な算術演算を再現するために、最上位ビットを特別扱いして演算しています。randomseed() の式がオリジナルよりも長いのは、符号の問題もありますが、それよりも32bit同士の乗算で演算誤差が出てしまったために*1、上位16bitと下位16bitに分割しながら乗算をおこなうようにしたという理由が大きいです。最後の +i はさすがにまず誤差が出ることはないように思います。でもこの際なので開き直って、ここも32bitの範囲を超えることがないようにしてみました。

これほど遠回りの多い書き方になると、これ絶対「高速」という利点を殺してるよとか、ライブラリで用意した方がいろんな意味で早いよとか、「まあそんなことより（自主規制）と（禁則事項）したい」とか言われそうですが、Lua学習の副産物程度の位置づけのものなので、そのあたりは深く考えないことにします。どうせ4さいあたまだし

なにかの役に立つことがあったら、コメントで教えてもらえるときっとうれしくなれるのです。

*1:今回のかけ算は演算結果の下の位（下位32bit）が重要でそれより上の位はどうでもいいのですが、算術的には当然上の位の方が重要視されるわけで、うまく収まらなかった下の位はすぱっと切られちゃうんですよね。かなしいっ！