' メルセンヌツイスタ、Xorshift、SFC32 のテスト

Option Explicit

Private Type Doub1: x As Double: End Type
Private Type Shrt4: x(0 To 3) As Integer: End Type
Private Type Long2: x(0 To 1) As Long: End Type

Public Sub WriteTest()
    Dim i&, j&, k&, clo#, min#, r#, fname$
    ChDrive ThisWorkbook.Path
    ChDir ThisWorkbook.Path
    fname = "mt31-vba.txt"
    InitMt 1
    Open fname For Output As #1
    Print #1, "1000 outputs of RandMt()"
    For i = 1 To 1000
        Print #1, Format$(RandMt(), "@@@@@@@@@@ ");
        If i Mod 5 = 0 Then Print #1, ""
    Next i
    Print #1, ""
    Print #1, "1000 outputs of UnifMt()"
    For i = 1 To 1000
        Print #1, Format$(UnifMt(), "0.00000000 ");
        If i Mod 5 = 0 Then Print #1, ""
    Next i
    InitXor 1
    Print #1, ""
    Print #1, "1000 outputs of RandXor()"
    For i = 1 To 1000
        Print #1, Format$(RandXor(), "@@@@@@@@@@ ");
        If i Mod 5 = 0 Then Print #1, ""
    Next i
    Print #1, ""
    Print #1, "1000 outputs of UnifXor()"
    For i = 1 To 1000
        Print #1, Format$(UnifXor(), "0.00000000 ");
        If i Mod 5 = 0 Then Print #1, ""
    Next i
    InitSfc 1
    Print #1, ""
    Print #1, "1000 outputs of RandSfc()"
    For i = 1 To 1000
        Print #1, Format$(RandSfc(), "@@@@@@@@@@ ");
        If i Mod 5 = 0 Then Print #1, ""
    Next i
    Print #1, ""
    Print #1, "1000 outputs of UnifSfc()"
    For i = 1 To 1000
        Print #1, Format$(UnifSfc(), "0.00000000 ");
        If i Mod 5 = 0 Then Print #1, ""
    Next i
    Close #1
    MsgBox "Wrote '" & fname$ & "'."
    Range("A1").Select
End Sub

Public Sub SpeedTest()
    Dim i&, j&, k&, iter&, line&, clo#, min#, sum#
    Const MBYTE = 0.000008
    Range("A:E").Delete: DoEvents
    [A1] = "RNG Name": [B1] = "Iteration": [C1] = " Speed": sum = 0
    InitMt 1: InitXor 1: InitSfc 1
    
    line = 2: iter = 200000: Cells(line, 1) = "UnifMt":
    min = 32767: Cells(line, 2) = iter
    For i = 1 To 10
        Cells(line, 3) = "' " & CStr(i) & "0%": DoEvents: clo = Timer
        For k = 1 To iter: sum = sum + UnifMt(): Next
        clo = Timer - clo
        If (clo < min) And (clo > 0) Then min = clo
    Next i
    Cells(line, 3) = " " & Format$(MBYTE * iter / min, "0.##0") & " MB/s"
    
    line = 3: iter = 500000: Cells(line, 1) = "UnifXor":
    min = 32767: Cells(line, 2) = iter
    For i = 1 To 10
        Cells(line, 3) = "' " & CStr(i) & "0%": DoEvents: clo = Timer
        For k = 1 To iter: sum = sum + UnifXor(): Next
        clo = Timer - clo
        If (clo < min) And (clo > 0) Then min = clo
    Next i
    Cells(line, 3) = " " & Format$(MBYTE * iter / min, "0.##0") & " MB/s"
    
    line = 4: iter = 200000: Cells(line, 1) = "UnifSfc":
    min = 32767: Cells(line, 2) = iter
    For i = 1 To 10
        Cells(line, 3) = "' " & CStr(i) & "0%": DoEvents: clo = Timer
        For k = 1 To iter: sum = sum + UnifSfc(): Next
        clo = Timer - clo
        If (clo < min) And (clo > 0) Then min = clo
    Next i
    Cells(line, 3) = " " & Format$(MBYTE * iter / min, "0.##0") & " MB/s"
    
    line = 5: iter = 5000: Cells(line, 1) = "[RAND()]"
    min = 32767: Cells(line, 2) = iter
    For i = 1 To 10
        Cells(line, 3) = "' " & CStr(i) & "0%": DoEvents: clo = Timer
        For k = 1 To iter: sum = sum + [RAND()]: Next
        clo = Timer - clo
        If (clo < min) And (clo > 0) Then min = clo
    Next i
    Cells(line, 3) = " " & Format$(MBYTE * iter / min, "0.##0") & " MB/s"
    
    line = 6: iter = 5000000: Cells(line, 1) = "Rnd"
    min = 32767: Cells(line, 2) = iter
    For i = 1 To 10
        Cells(line, 3) = "' " & CStr(i) & "0%": DoEvents: clo = Timer
        For k = 1 To iter: sum = sum + Rnd(): Next
        clo = Timer - clo
        If (clo < min) And (clo > 0) Then min = clo
    Next i
    Cells(line, 3) = " " & Format$(MBYTE * iter / min, "0.0") & " MB/s"
    
    Range("A1").Select
    MsgBox "Done"
End Sub

Private Sub SHR2(ByRef g As Long2, ByVal k As Byte)
    If k = 0 Then Exit Sub
    If k >= 64 Then g.x(0) = 0: g.x(1) = 0: Exit Sub
    If k > 32 Then g.x(0) = LogicalShiftRight(g.x(1), k - 32): g.x(1) = 0: Exit Sub
    If k = 32 Then g.x(0) = g.x(1): g.x(1) = 0: Exit Sub
    g.x(0) = LogicalShiftRight(g.x(0), k) + LogicalShiftLeft(g.x(1), 32 - k): g.x(1) = LogicalShiftRight(g.x(1), k)
End Sub

Private Sub SHL2(ByRef g As Long2, ByVal k As Byte)
    If k = 0 Then Exit Sub
    If k >= 64 Then g.x(0) = 0: g.x(1) = 0: Exit Sub
    If k > 32 Then g.x(1) = LogicalShiftLeft(g.x(0), k - 32): g.x(0) = 0: Exit Sub
    If k = 32 Then g.x(1) = g.x(0): g.x(0) = 0: Exit Sub
    g.x(1) = LogicalShiftLeft(g.x(1), k) + LogicalShiftRight(g.x(0), 32 - k): g.x(0) = LogicalShiftLeft(g.x(0), k)
End Sub

Private Sub DisplayResult(a() As Double, ByVal n As Long, name As String)
    Const mega = 2& ^ 20: Static cnt(0 To 63) As Long
    Dim d As Doub1, s As Shrt4, g As Long2, i, j, k As Long
    Application.ScreenUpdating = True: Range("A:J").Delete
    [B1] = name: [C1] = "' n=" & CStr(n): DoEvents
    For i = 0 To 63: cnt(i) = 0: Next
    For i = 0 To n - 1
        d.x = a(i):  LSet g = d
        If d.x <> 0 Then
            k = 1012 - g.x(1) \ mega
            g.x(1) = (g.x(1) And (mega - 1)) Or mega
            If (k > 0) Then SHR2 g, k
            If (k < 0) Then SHL2 g, -k
        End If
        LSet s = g
        For k = 0 To 63
            If (s.x(k \ 16) And (2& ^ (k And 15))) <> 0 Then: cnt(k) = cnt(k) + 1
        Next k
    Next i
    For i = 0 To 63
        Cells(i \ 4 + 2, i Mod 4 + 1) = "' " & i - 63 & "=" & cnt(i)
    Next
    MsgBox "Done " & name
End Sub

Public Sub BitTest()
    Const n = 10000: Static a(0 To n - 1) As Double: Dim i As Long
    InitMt 1: InitXor 1: InitSfc 1
    For i = 0 To n - 1: a(i) = UnifMt: Next
    DisplayResult a, n, "UnifMt"
    For i = 0 To n - 1: a(i) = UnifXor: Next
    DisplayResult a, n, "UnifXor"
    For i = 0 To n - 1: a(i) = UnifSfc: Next
    DisplayResult a, n, "UnifSfc"
    For i = 0 To n - 1: a(i) = Rnd: Next
    DisplayResult a, n, "Rnd"
    For i = 0 To n - 1: a(i) = [RAND()]: Next
    DisplayResult a, n, "[RAND()]"
End Sub

#If False Then
#include <stdio.h>
#define N 644
#define M 322
#define MATRIX_A 0x4c745d7eL
#define UPPER_MASK 0x78000000L
#define LOWER_MASK 0x07FFFFFFL
static long mt[N];
static int mti=N+1;
static long xx,yy,zz,ww;
static unsigned long aa,bb,cc,cnt;

void InitMt(long s) {
    mt[0]= s & 0x7fffffffL;
    for (mti=1; mti<N; mti++) {
        mt[mti] = (1812433253L * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
        mt[mti] &= 0x7fffffffL;
    }
}

long RandMt(void) {
    long y;
    if (mti >= N) {
        int kk;
        if (mti == N+1) InitMt(1);
        for (kk=0;kk<N-M;kk++) {
            y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
            mt[kk] = mt[kk+M] ^ (y >> 1) ^ (y&1?MATRIX_A:0);
        }
        for (;kk<N-1;kk++) {
            y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
            mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ (y&1?MATRIX_A:0);
        }
        y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
        mt[N-1] = mt[M-1] ^ (y >> 1) ^ (y&1?MATRIX_A:0);
        mti = 0;
    }
    y = mt[mti++];
    y ^= (y >> 12);
    y ^= (y <<  7) & 0x5b5b9b80L;
    y ^= (y << 15) & 0x7ef20000L;
    y ^= (y >> 18);
    return y;
}

double UnifMt(void) {
    long x=RandMt()>>9;
    return(RandMt()*4194304.0+x)*(1.0/9007199254740992.0);
}

void InitXor(long s) {
    xx=s&0x7fffffffL;
    yy=(1812433253L*(xx^(xx>>30))+1)&0x7fffffffL;
    zz=(1812433253L*(yy^(yy>>30))+2)&0x7fffffffL;
    ww=(1812433253L*(zz^(zz>>30))+3)&0x7fffffffL;
}

long RandXor(void) {
    enum { a=28, b=17, c=3 }; long t=xx;
    xx&=(1L<<(31-a))-1; t^=xx<<a; xx=yy; yy=zz; zz=ww;
    t^=t>>b; ww^=ww>>c; ww^=t; return ww;
}

double UnifXor(void) {
    long x=RandXor()>>9;
    return(RandXor()*4194304.0+x)*(1.0/9007199254740992.0);
}

long RandSfc(void) {
    unsigned long tmp=(aa+bb+cnt++)&0xffffffffUL;
    aa=bb^(bb>>9);
    bb=(cc+(cc<<3))&0xffffffffUL;
    cc=(((cc<<21)|(cc>>(32-21)))+tmp)&0xffffffffUL;
    return tmp&0x7fffffffL;
}

void InitSfc(long s) {
    int i;
    aa=cc=0; bb=s; cnt=1;
    for (i=0;i<12;i++) RandSfc();
}

double UnifSfc(void) {
    long x=RandSfc()>>9;
    return(RandSfc()*4194304.0+x)*(1.0/9007199254740992.0);
}

int main(void) {
    int i;
    InitMt(1);
    printf("1000 outputs of RandMt()\n");
    for (i=0; i<1000; i++) {
      printf("%10ld ",RandMt());
      if (i%5==4) printf("\n");
    }
    printf("\n1000 outputs of UnifMt()\n");
    for (i=0; i<1000; i++) {
      printf("%10.8f ", UnifMt());
      if (i%5==4) printf("\n");
    }
    InitXor(1);
    printf("\n1000 outputs of RandXor()\n");
    for (i=0; i<1000; i++) {
      printf("%10ld ",RandXor());
      if (i%5==4) printf("\n");
    }
    printf("\n1000 outputs of UnifXor()\n");
    for (i=0; i<1000; i++) {
      printf("%10.8f ", UnifXor());
      if (i%5==4) printf("\n");
    }
    InitSfc(1);
    printf("\n1000 outputs of RandSfc()\n");
    for (i=0; i<1000; i++) {
      printf("%10ld ",RandSfc());
      if (i%5==4) printf("\n");
    }
    printf("\n1000 outputs of UnifSfc()\n");
    for (i=0; i<1000; i++) {
      printf("%10.8f ", UnifSfc());
      if (i%5==4) printf("\n");
    }
    return 0;
}
#End If
#If False Then
/* メルセンヌツイスタ・ダイナミッククリエーター(dcmt0.6.7) */
#include <stdio.h>
#include "../include/dc.h"
int main(void) {
    mt_struct *mts;
    init_dc(4172);
    mts = get_mt_parameter(31,19937);
    if (mts == NULL) { printf ("error\n"); return 1; }
    printf("aaa     =0x%08x\n",mts->aaa);
    printf("mm      =%d\n",mts->mm);
    printf("nn      =%d\n",mts->nn);
    printf("rr      =%d\n",mts->rr);
    printf("ww      =%d\n",mts->ww);
    printf("wmask   =0x%08x\n",mts->wmask);
    printf("umask   =0x%08x\n",mts->umask);
    printf("lmask   =0x%08x\n",mts->lmask);
    printf("shift0  =%d\n",mts->shift0);
    printf("shift1  =%d\n",mts->shift1);
    printf("shiftB  =%d\n",mts->shiftB);
    printf("shiftC  =%d\n",mts->shiftC);
    printf("maskB   =0x%08x\n",mts->maskB);
    printf("maskC   =0x%08x\n",mts->maskC);
    free_mt_struct(mts);
    return 0;
}
#End If