' ＶＢＡによるXorshift

Option Explicit

' 定数の宣言
Private Const P32 = 2# ^ 32, P31 = 2 ^ 31, P22 = 2# ^ 22, P18 = 2 ^ 18
Private Const P13 = 2 ^ 13, P9 = 2 ^ 9
Private Const M53 = 2# ^ -53, M32 = 2# ^ -32, M30 = 2# ^ -30
Private Const K18 = (2 ^ 18) - 1


' 乱数の状態
Private x, y, z, w As Long
Private flag As Boolean


' 初期化の補助関数
Private Function Ri(ByRef r As Double, ByVal i As Long) As Long
    Dim s As Variant
    Dim shft As Double
    Dim a As Long
    If r >= P31 Then a = r - P32 Else a = r
    a = a Xor Int(r * M30)
    If a < 0 Then r = a + P32 Else r = a
    s = 1812433253 * CDec(r) + i: r = s - CDec(Int(s * M32)) * P32
    If r >= P31 Then Ri = r - P31 Else Ri = r
End Function


' s を種にして乱数を初期化する
Public Sub InitXor(ByVal s As Long)
    Dim r As Double
    If s < 0 Then r = P32 + s Else r = s
    flag = True: x = Ri(r, 1): y = Ri(r, 2): z = Ri(r, 3): w = Ri(r, 4)
End Sub


' 31 ビットの整数乱数
Public Function NextXor() As Long
    Dim t As Long
    If flag = 0 Then InitXor (1)
    t = x: x = y: y = z: z = w: t = t Xor ((t And K18) * P13)
    t = t Xor (t \ P18): w = w Xor t: NextXor = w
End Function


' 0 以上 1 未満の浮動小数点数乱数を返す
Public Function NextUnifXor() As Double
    Dim x As Long
    x = NextXor \ P9: NextUnifXor = (NextXor * P22 + x) * M53
End Function


' 時間を種にして乱数を初期化する
Public Sub RandomizeXor()
    InitXor (Timer() * 200)
End Sub