' ＶＢＡによるメルセンヌツイスタ、Xorshift、SFC32

Option Explicit

' メルセンヌツイスタのパラメータ(ダイナミッククリエーター(dcmt0.6.7)の結果)
Private Const MTN = 644, MTM = 322, MTA = 12, MTB = 7, MTC = 15, MTD = 18
Private Const MXA = &H4C745D7E, UMK = &H78000000, LMK = &H7FFFFFF
Private Const MKB = &H5B5B9B80, MKC = &H7EF20000

' Xorshift(124) のパラメータ
Private Const XA = 28, XB = 17, XC = 3

' 補助的な定数の宣言
Private Const LngMax = &H7FFFFFFF, LngMin = &H80000000
Private Const MTL = MTN - MTM, MTK = MTN - 1, MTJ = MTL - 1, MTP = MTN - 2
Private Const PWA = 2& ^ MTA, PWB = 2& ^ MTB, PWC = 2& ^ MTC, PWD = 2& ^ MTD
Private Const KB = MKB \ PWB, KC = MKC \ PWC
Private Const P22 = 2& ^ 22, P9 = 512
Private Const P32 = 2# ^ 32, P31 = 2# ^ 31, Pc31 = 2@ ^ 31
Private Const M53 = 2# ^ -53, M31 = 2# ^ -31, M30 = 2# ^ -30
Private Const PA = 2& ^ XA, PB = 2& ^ XB, PC = 2& ^ XC
Private Const XMA = (2& ^ (31 - XA)) - 1
Private Const P28 = 2& ^ 28, P21 = 2& ^ 21, P20 = 2& ^ 20
Private Const P11 = 2048, P10 = 1024, P3 = 8
#If Win64 Then
Private Const Pq32 = 2^ ^ 32, Pq31 = 2^ ^ 31
#End If

' 乱数の状態
Private mt(0 To MTK), mti As Long
Private xx, yy, zz, ww As Long, okxor As Boolean
Private aa, bb, cc, cnt As Long, oksfc As Boolean

' 初期化の補助関数
Private Function Ri(ByRef r As Long, ByVal i As Long) As Long
    Dim c As Currency
    r = r Xor Int(r * M30): c = CDec(r) * 181243.3253:
    c = c - Fix(c * M31) * Pc31: c = c * 10000 + i
    c = c - Fix(c * M31) * Pc31: r = CLng(c): Ri = r
End Function

' s を種にして乱数を初期化する
Public Sub InitMt(ByVal s As Long)
    s = s And LngMax: mt(0) = s
    For mti = 1 To MTK: mt(mti) = Ri(s, mti): Next mti
    mti = MTN
End Sub

' 31 ビットの整数乱数
Public Function RandMt() As Long
    Dim y, k As Long
    If mti = 0 Then InitMt 1
    If mti = MTN Then
        mti = 0
        For k = 0 To MTJ
            y = (mt(k) And UMK) Or (mt(k + 1) And LMK)
            mt(k) = mt(k + MTM) Xor (y \ 2) Xor (-(y And 1) And MXA)
        Next k
        For k = MTL To MTP
            y = (mt(k) And UMK) Or (mt(k + 1) And LMK)
            mt(k) = mt(k - MTL) Xor (y \ 2) Xor (-(y And 1) And MXA)
        Next k
        y = (mt(MTK) And UMK) Or (mt(0) And LMK)
        mt(MTK) = mt(MTM - 1) Xor (y \ 2) Xor (-(y And 1) And MXA)
    End If
    y = mt(mti): mti = mti + 1
    y = y Xor (y \ PWA): y = y Xor ((y And KB) * PWB)
    y = y Xor ((y And KC) * PWC): y = y Xor (y \ PWD): RandMt = y
End Function

' 0 以上 1 未満の乱数を返す
Public Function UnifMt() As Double
    Dim x As Long
    x = RandMt \ P9: UnifMt = (RandMt * P22 + x) * M53
End Function

' 時間を種にして乱数を初期化する
Public Sub RandomizeMt()
    InitMt (CLng(Timer() * 1000))
End Sub

Public Sub InitXor(ByVal s As Long)
    okxor = True: s = s And LngMax: xx = s
    yy = Ri(s, 1): zz = Ri(s, 2): ww = Ri(s, 3)
End Sub

Public Function RandXor() As Long
    Dim t As Long
    If Not okxor Then InitXor 1
    t = xx: xx = xx And XMA: t = t Xor (xx * PA)
    xx = yy: yy = zz: zz = ww: t = t Xor (t \ PB)
    ww = ww Xor (ww \ PC): ww = ww Xor t: RandXor = ww
End Function

Public Function UnifXor() As Double
    Dim x As Long
    x = RandXor \ P9: UnifXor = (RandXor * P22 + x) * M53
End Function

Public Sub RandomizeXor()
    InitXor (CLng(Timer() * 1000))
End Sub

Public Function RandSfc()
#If Win64 Then
    Dim tmp, s, t As Long, q As LongLong
    If Not oksfc Then InitSfc 1
    q = CLngLng(aa) + CLngLng(bb) + CLngLng(cnt)
    If q >= Pq31 Then q = q - Pq32
    If q < -Pq31 Then q = q + Pq32
    tmp = CLng(q): t = bb
    If cnt = LngMax Then cnt = LngMin Else cnt = cnt + 1
    If t >= 0 Then t = t \ P9 Else t = ((t And LngMax) \ P9) Or P22
    aa = bb Xor t: t = cc: t = (t And (P28 - 1)) * P3
    If (cc And P28) = P28 Then t = t Or LngMin
    q = CLngLng(cc) + CLngLng(t)
    If q >= Pq31 Then q = q - Pq32
    If q < -Pq31 Then q = q + Pq32
    bb = CLng(q): t = cc
    If t >= 0 Then t = t \ P11 Else t = ((t And LngMax) \ P11) Or P20
    s = cc: s = (s And (P10 - 1)) * P21
    If (cc And P10) = P10 Then s = s Or LngMin
    s = s Or t: q = CLngLng(s) + CLngLng(tmp)
    If q >= Pq31 Then q = q - Pq32
    If q < -Pq31 Then q = q + Pq32
    cc = CLng(q): RandSfc = tmp And LngMax
#Else
    Dim tmp, s, t As Long, r As Double
    If Not oksfc Then InitSfc 1
    r = CDbl(aa) + CDbl(bb) + CDbl(cnt): t = bb
    If r >= P31 Then r = r - P32
    If r < -P31 Then r = r + P32
    tmp = CLng(r)
    If cnt = LngMax Then cnt = LngMin Else cnt = cnt + 1
    If t >= 0 Then t = t \ P9 Else t = ((t And LngMax) \ P9) Or P22
    aa = bb Xor t: t = cc: t = (t And (P28 - 1)) * P3
    If (cc And P28) = P28 Then t = t Or LngMin
    r = CDbl(cc) + CDbl(t)
    If r >= P31 Then r = r - P32
    If r < -P31 Then r = r + P32
    bb = CLng(r): t = cc
    If t >= 0 Then t = t \ P11 Else t = ((t And LngMax) \ P11) Or P20
    s = cc: s = (s And (P10 - 1)) * P21
    If (cc And P10) = P10 Then s = s Or LngMin
    r = CDbl(s Or t) + CDbl(tmp)
    If r >= P31 Then r = r - P32
    If r < -P31 Then r = r + P32
    cc = CLng(r): RandSfc = tmp And LngMax
#End If
End Function

Public Sub InitSfc(ByVal s As Long)
    Dim i As Long
    oksfc = True: aa = 0: cc = 0: bb = s: cnt = 1
    For i = 1 To 12: Call RandSfc: Next
End Sub

Public Function UnifSfc() As Double
    Dim x As Long
    x = RandSfc \ P9: UnifSfc = (RandSfc * P22 + x) * M53
End Function

Public Sub RandomizeSfc()
    InitSfc (CLng(Timer() * 1000))
End Sub