/* >>68 さん。的確なアドバイスありがとうございます。>>68 さんの方法で書いてみました。*/ /* 種で初期化する関数付き。XorShift の SSE2 による高速化はあきらめるしかないのか！ */ #include <stdio.h> #include <time.h> #include <emmintrin.h> unsigned long xor128(void) { static unsigned long x=123456789UL,y=362436069UL,z=521288629UL,w=88675123UL; unsigned long t;t=(x^(x<<11));x=y;y=z;z=w;return(w=(w^(w>>19))^(t^(t>>8)));}
int idx=4; union { unsigned long v[4]; __m128i m; } x={123456789,123456789,123456789,123456789},y={362436069,362436069,362436069,362436069}, z={521288629,521288629,521288629,521288629},w={ 88675123, 88675123, 88675123, 88675123};
void sxor128x4(unsigned long s) {int i; for (idx=4,i=0;i<16;i++) x.v[i]=s=1812433253*(s^(s>>30))+i; }
/* Complementary Multiply With Carry */ void initialize( unsigned long ); unsigned long cmwc( void ); enum { A = 3636507990, B = 0xffffffff, R = 1024 }; unsigned long seed[ R ], index, c;
void initialize( unsigned long n ) { unsigned long i; for (i=0; i<R; i++) { seed[ i ] = n; n = n * 1103515245 + 12345; } index = c = 0; } unsigned long cmwc( void ) { unsigned long x; unsigned long long t; t = (unsigned long long)seed[ index ] * A + c; x = B - (unsigned long)( t & B ); c = (unsigned long)( t >> 32 ); seed[ index ] = x; if ( ++index >= R ) index = 0; return x; }
' Xorshift を VBA で書く場合、Double を使うと比較的シンプルになる ' VBA の演算子 Xor と関数 Fix は 31 bit までしか扱えないので ' 53 bit まで扱える関数 uFix、uXor を用意した。 Option Explicit Public x, y, z, w As Double
Function uFix(ByVal x As Double) As Double Const Base = 2# ^ 22 Dim y As Long y = Int(x / Base): uFix = y * Base + Int(x - y * Base) End Function
Function uXor(ByVal i As Double, ByVal j As Double) As Double Const Base = 2# ^ 22 Dim u, v As Long u = Int(i / Base): v = Int(j / Base): i = i - u * Base: j = j - v * Base uXor = (u Xor v) * Base + (Int(i) Xor Int(j)) End Function
Function xor128() As Double Dim t As Double t = x * 2# ^ 11: t = t - Int(t / 2# ^ 32) * 2# ^ 32 t = uXor(x, t): x = y: y = z: z = w: t = uXor(t, uFix(t / 2# ^ 8)) w = uXor(uXor(w, uFix(w / 2# ^ 19)), t): xor128 = w End Function
Sub Test() Dim i As Long x = 123456789#: y = 362436069#: z = 521288629#: w = 88675123# For i = 1 To 2000: Cells(i, 1) = xor128: Next End Sub
116 名前：115 mailto:sage [2008/08/13(水) 18:47:55 ]
' 書き込んだ後、uFix は必要がないことに気づきました。 ' こっちの方が高速です。 Option Explicit Public x, y, z, w As Double
Function uXor(ByVal i As Double, ByVal j As Double) As Double Const Base = 2# ^ 30 Dim u, v As Long u = Int(i / Base): v = Int(j / Base): i = i - u * Base: j = j - v * Base uXor = (u Xor v) * Base + (Int(i) Xor Int(j)) End Function
Function xor128() As Double Dim t As Double t = x * 2# ^ 11: t = t - Int(t / 2# ^ 32) * 2# ^ 32 t = uXor(x, t): x = y: y = z: z = w: t = uXor(t, Int(t / 2# ^ 8)) w = uXor(uXor(w, Int(w / 2# ^ 19)), t): xor128 = w End Function
Sub Test() Dim i As Long x = 123456789#: y = 362436069#: z = 521288629#: w = 88675123# For i = 1 To 2000: Cells(i, 1) = xor128: Next End Sub
'種を使って初期化できます。116よりさらに２倍程度高速になりました Private Const p32 = 2# ^ 32, p31 = 2# ^ 31, p21 = 2# ^ 21, p11 = 2# ^ 11 Private Const m53 = 2# ^ -53, m32 = 2# ^ -32, m30 = 2# ^ -30, m19 = 2# ^ -19, m11 = 2# ^ -11, m8 = 2# ^ -8 Private x, y, z, w As Double Private f As Boolean Private Function uXor(ByVal x As Double, ByVal y As Double) As Double Dim u, v As Long If x >= p31 Then u = x - p32 Else u = x If y >= p31 Then v = y - p32 Else v = y u = u Xor v If u < 0 Then uXor = u + p32 Else uXor = u End Function Private Function XSub(ByVal x As Double, ByVal i As Long) As Double Dim s As Variant s = CDec(1812433253) * CDec(uXor(x, Int(x * m30))) + CDec(i): s = s - CDec(Int(s * m32)) * CDec(p32): XSub = s End Function Public Sub InitXor(ByVal s As Long) ' s を種にして乱数を初期化する f = True: x = XSub(s, 1): y = XSub(x, 2): z = XSub(y, 3): w = XSub(z, 4) End Sub Private Function NextXor() As Double Dim t As Double If f = 0 Then InitXor (1) t = x * p11: t = t - Int(t * m32) * p32: t = uXor(x, t): x = y: y = z: z = w: t = uXor(t, Int(t * m8)): w = uXor(uXor(w, Int(w * m19)), t): NextXor = w End Function Public Function NextUnif() As Double ' 0 以上 1 未満の乱数を返す Dim x, y As Double x = NextXor * m11: y = NextXor: NextUnif = (y * p21 + Int(x)) * m53 End Function
MTならFORTRANだけど ttp://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/VERSIONS/FORTRAN/fortran.html この中に revrand(): it generates prn in reverse order って書いてあるけど
#define main dummy #include "mt19937ar.c" #undef main
unsigned long revgrnd(void) { static unsigned long mag01[2]={0x0UL,MATRIX_A}; unsigned long y=mt[--mti],z,p,q,mt0L,mt0; int x,kk; if (mti==0) { z=mt[N-1]^mt[M-1]; x=(int)(z>>31); y=z^mag01[x]; p=(y<<1)|x; mt0L=LOWER_MASK&p; mt[0]=(mt[0]&UPPER_MASK)^mt0L; mt0=mt[0]; for (kk=N-1;kk>N-M;kk--) { z=mt[kk-1]^mt[kk-1+M-N]; x=(int)(z>>31); y=z^mag01[x]; q=(y<<1)|x; mt[kk]=(UPPER_MASK&p)|(LOWER_MASK&q); p=q; } for (;kk;kk--) { z=mt[kk-1]^mt[kk-1+M]; x=(int)(z>>31); y=z^mag01[x]; q=(y<<1)|x; mt[kk]=(UPPER_MASK&p)|(LOWER_MASK&q); p=q; } mt[0]=(UPPER_MASK&p)|mt0L; y=mt0; mti=N; } y^=(y>>11); y^=(y<<7)&0x9d2c5680UL; y^=(y<<15)&0xefc60000UL; y^=(y>>18); return y; }
int main(void) { int i; static unsigned long x[5000],y[5000]; for (i=0;i<5000;i++) x[i]=genrand_int32(); for (i=4999;i>=0;i--) y[i]=revgrnd(); for (i=0;i<5000;i++) if (x[i]!=y[i]) { printf("ERROR\n"); break; } for (i=0;i<10;i++) printf("%08lx %08lx\n",x[i],y[i]); return 0; }
