/* dxor.c :          coded by isaku@pb4.so-net.ne.jp 2008/12〜2009/1
乱数 Xorshift のSSE2を使った高速化(浮動小数点数版)

dxor128 関数は速度の比較のために、論文の Xorshift から、単純に
double に変換したもの。周期は 2^{128}-1 で解像度は 32 bit である。

dxor156 は 52 bit Xorshift 演算
 t=x^x<<19; x=y; y=z; z^=t^t>>33;
を行い、3つの xmm 型 x,y,z に状態を保持する。
xmm 型の上位 64 bit は使用していない。
周期は 2^{156}-1 で、解像度は 52 bit である。
計算にかかる時間は Visual C++ 2008、Core 2 Quad で
dxor128 の 35% 以下になった。
Celeron D では dxor128 の 53% の計算時間になった(速くなった)。

また、Diehard テストを行ったが、良好な結果を得た。

高速化ではないが、さらに高解像度な乱数も作成した。
dxor189 は 63 bit Xorshift 演算
 t=x^x<<6; x=y; y=z; z^=t^t>>19;
を行う。戻り値の型は long double である。
周期は 2^{189}-1 で、解像度は 63 bit である。
Visual C++ では long double は 8 byte なので無効である。
一方 gcc では 12 byte or 16 byte である。
しかし、いずれの場合でも 10 byte のみしか使っておらず、
内部表現は全く同じである。

以下のプログラムは dxor128 , dxor156 , dxor189 のそれぞれについて、
1億個 × 10の乱数を足し合わせて、最速の実行時間を比較する。

How to make:
(VC++2008) set VS=C:\Program Files\Microsoft Visual Studio 9.0
           set SDK=C:\Program Files\Microsoft SDKs\Windows\v6.0A
           set PATH=%VS%\Common7\IDE;%VS%\VC\BIN;%SDK%\bin
           set INCLUDE=%VS%\VC\INCLUDE;%SDK%\include
           set LIB=%VS%\VC\LIB;%SDK%\lib
           cl /O2 /GL /W4 dxor.c
(gcc-4 on cygwin)   gcc-4 -O3 -Wall -msse3 dxor.c
(gcc 4 on linux)    gcc   -O3 -Wall -msse3 dxor.c

1億回の生成時間(単位 秒)
+-----------------+------------------+-------+-------+-------+
| machine         |compiler & system |dxor128|dxor156|dxor189|
+-----------------+------------------+-------+-------+-------+
|Core 2 Quad Q9450|VC++2008 on XP    | 0.734 | 0.234 | ----- |
|Core 2 Quad Q9450|gcc-4 on cygwin   | 1.046 | 0.328 | 0.859 |
|Core 2 Quad Q9450|gcc 4 on linux    | 1.280 | 0.430 | 0.940 |
|Core 2 Quad Q9450|gcc 4 on linux x64| 0.490 | 0.470 | 0.960 |
|Celeron D 2.66GHz|VC++2008 on XP    | 1.140 | 0.593 | ----- |
|Celeron D 2.66GHz|gcc-4 on cygwin   | 3.360 | 0.907 | 3.047 |
|Celeron D 2.66GHz|gcc 4 on linux    | 3.420 | 1.490 | 4.480 |
|Celeron D 2.66GHz|gcc 4 on linux x64| 1.590 | 1.650 | 3.450 |
+-----------------+------------------+-------+-------+-------+
*/

#include <stdio.h>
#include <time.h>

__inline double dxor128(void) {
    static unsigned x=123456789,y=362436069,z=521288629,w=88675123;
    unsigned t=x^x<<11; x=y; y=z; z=w; w^=w>>19^t^t>>8;
    return w*(1.0/4294967296.0);
}

#if _MSC_VER >= 1300 || defined(__GNUC__)
#include <emmintrin.h>

#define AVALABLE_D156

__inline double dxor156(void) {
    static union { int i[4]; double d[2]; __m128i m; }
    x = {{123456789,(362436069&0xFFFFF)|0x3FF00000,0,0}},
    y = {{521288629,( 88675123&0xFFFFF)|0x3FF00000,0,0}},
    z = {{  5783321,(  6615241&0xFFFFF)|0x3FF00000,0,0}};
    const union { int i[4]; __m128i m; } mask={{-1,0xFFFFF,0,0}};
    __m128i t=x.m, s;
    s=_mm_slli_epi64(t,19);     t=_mm_xor_si128(t,s);
    t=_mm_and_si128(t,mask.m);  s=_mm_srli_epi64(t,33);
    s=_mm_xor_si128(s,t); t=y.m; x.m=t; t=z.m; y.m=t;
    t=_mm_xor_si128(t,s); z.m=t; return z.d[0]-1.0;
}

#ifdef __GNUC__

#define AVALABLE_D189

__inline long double dxor189(void) {
    static union { int i[4]; long double d; __m128i m; }
    x = {{123456789,362436069|0x80000000,0x3FFF,0}},
    y = {{521288629, 88675123|0x80000000,0x3FFF,0}},
    z = {{  5783321,  6615241|0x80000000,0x3FFF,0}};
    const union { int i[4]; __m128i m; } mask={{-1,0x7FFFFFFF,0,0}};
    __m128i t=x.m, s;
    s=_mm_slli_epi64(t,6);      t=_mm_xor_si128(t,s);
    t=_mm_and_si128(t,mask.m);  s=_mm_srli_epi64(t,19);
    s=_mm_xor_si128(s,t); t=y.m; x.m=t; t=z.m; y.m=t;
    t=_mm_xor_si128(t,s); z.m=t; return z.d-1.0;
}

#endif
#endif

#define N 100000000

int main(void) {
    clock_t lap,min1,min2; unsigned i,k; long double sum;

    min1=999999999; sum=0;
    for (k=0;k<10;k++) {
        lap=clock();
        for (i=0;i<N;i++) sum+=dxor128();
        lap=clock()-lap; printf(" %3lu",lap);
#ifdef __GNUC__
        fflush(stdout);
#endif
        if (lap<min1) min1=lap;
    }
    printf("(%lu)%.0Lf\n",min1,sum);

#ifdef AVALABLE_D156

    min2=999999999; sum=0;
    for (k=0;k<10;k++) {
        lap=clock();
        for (i=0;i<N;i++) sum+=dxor156();
        lap=clock()-lap; printf(" %3lu",lap);
#ifdef __GNUC__
        fflush(stdout);
#endif
        if (lap<min2) min2=lap;
    }
    printf("(%lu)%.0Lf ",min2,sum);
    if (min2<=min1)
         printf("%.2f%%fast\n",100.0*(min1-min2)/min1);
    else printf("%.2f%%slow\n",100.0*(min2-min1)/min1);
#endif

#ifdef AVALABLE_D189
    min2=999999999; sum=0;
    for (k=0;k<10;k++) {
        lap=clock();
        for (i=0;i<N;i++) sum+=dxor189();
        lap=clock()-lap; printf(" %3lu",lap);
#ifdef __GNUC__
        fflush(stdout);
#endif
        if (lap<min2) min2=lap;
    }
    printf("(%lu)%.0Lf ",min2,sum);
    if (min2<=min1)
         printf("%.2f%%fast\n",100.0*(min1-min2)/min1);
    else printf("%.2f%%slow\n",100.0*(min2-min1)/min1);
#endif
    return 0;
}