/* xor.c : coded by isaku@pb4.so-net.ne.jp
2008/12〜2009/1,  2010/10, 2012/6, 2014/8
乱数 Xorshift のSSE2を使った高速化(整数版)

xor128 は論文に書いてあった 128 bit Xorshift である。

以前の xor130 は SSE2 を使って、128ビット演算を行うものであった。
一度に4個の整数を得るので、周期は 2^{130}-4 であった。
1ビット単位のシフト演算は SSE2 ではサポートされていないので、
他の演算を組み合わせることで128ビットシフト演算を実現した。

最初のバージョンでは

enum { a=11,b=23 };
x^=x<<a; x^=x>>b; return x;

のように、二つの Xorshift 演算を行っていた。
Diehard テストを行ったところ、良好な結果を得たが、
Dieharder テストでは必ず不合格が出るので、
次のバージョンでは、三つの Xorshift 演算を行うように改良した。

enum { a=8,b=69,c=75 };
x^=x<<a; x^=x>>b; x^=x<<c; return x;

高速に計算できるパラメータを選んだので、
古いバージョンと同じスピードで実行できた。
古い Dieharder テストで不合格がでなくなったが、
新しい Dieharder テスト(3.31.1)では不合格になってしまった。

そこで、

enum { a=101,b=99,c=8 };
t=x^x<<a; x=y; y=z; z=w; t^=t>>b; w^=w>>c^t; return w;
とした xor514 にしたところ、新しい Dieharder テストでも不合格が出なくなった。
周期は 2^{514}-4 であ

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

10億回の生成時間(単位 秒)の結果

Core i7 2600(turbo boost off)
+-----------------------------------------------+------+------+-------+
|compiler and system                            |xor128|xor514|back514|
+-----------------------------------------------+------+------+-------+
|VC++2010 64bit cl /O2 /GL /W4 on Windows7 64bit| 1.919| 1.326|  1.544|
|VC++2010 32bit cl /O2 /GL /W4 on Windows7 64bit| 2.262| 1.326|  1.248|
|gcc -O4 -Wall             on Ubuntu 14.04 64bit| 1.809| 1.009|  1.115|
|gcc -O4 -Wall -m32 -msse2 on Ubuntu 14.04 64bit| 2.333| 1.081|  1.088|
+-----------------------------------------------+------+------+-------+
*/

#include <emmintrin.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <stdint.h>
#define RND_SSE2_NOINLINE 0 /* Turn it 1 if old gcc reports some warnnings */

/*
 Xorshift RNGs. George Marsaglia.
 Journal of statistical software, Vol 8, Issue 14, Jul 2003
 http://www.jstatsoft.org/v08/i14/paper   page 5
*/

__inline uint32_t xor128(void) {
    enum { a=11,b=8,c=19 };
    static uint32_t x=123456789,y=362436069,z=521288629,w=88675123;
    uint32_t t=x^x<<a; x=y; y=z; z=w; return w^=w>>c^t^t>>b;
}

typedef union tag514 { uint32_t u[20]; __m128i m[5]; } xor514_t;

void sxor514(xor514_t*p,uint32_t s) {
   uint32_t i; p->u[0]=20; p->u[4]=s;
   for (i=1;i<16;i++) p->u[i+4]=s=1812433253*(s^(s>>30))+i;
}


#if defined(__GNUC__)

 #define XOR514HEAD 4

 #define XOR514(X,W)                                             \
    t=_mm_xor_si128(X,_mm_slli_si128(_mm_slli_epi64(X,a-64),8)); \
    t=_mm_xor_si128(t,_mm_srli_si128(_mm_srli_epi64(t,b-64),8)); \
    X=_mm_xor_si128(W,_mm_xor_si128(_mm_srli_si128(W,c/8),t))

 #define BACK514(X,W)                                            \
    t=_mm_xor_si128(_mm_xor_si128(X,W),_mm_srli_si128(W,c/8));   \
    t=_mm_xor_si128(t,_mm_srli_epi64(_mm_srli_si128(t,8),b-64)); \
    X=_mm_xor_si128(t,_mm_slli_epi64(_mm_slli_si128(t,8),a-64))

 #if RND_SSE2_NOINLINE
extern void xor514sub(xor514_t*p) __attribute__((noinline));
extern void back514sub(xor514_t*p) __attribute__((noinline));
 #endif

void xor514sub(xor514_t*p) {
    enum { a=101,b=99,c=8 }; __m128i t,x,y,z,w;
    x=_mm_load_si128(&p->m[1]);
    y=_mm_load_si128(&p->m[2]);
    z=_mm_load_si128(&p->m[3]);
    w=_mm_load_si128(&p->m[4]);
    XOR514(x,w); XOR514(y,x); XOR514(z,y); XOR514(w,z);
    _mm_store_si128(&p->m[1],x);
    _mm_store_si128(&p->m[2],y);
    _mm_store_si128(&p->m[3],z);
    _mm_store_si128(&p->m[4],w);
}

void back514sub(xor514_t*p) {
    enum { a=101,b=99,c=8 }; __m128i t,x,y,z,w;
    x=_mm_load_si128(&p->m[1]);
    y=_mm_load_si128(&p->m[2]);
    z=_mm_load_si128(&p->m[3]);
    w=_mm_load_si128(&p->m[4]);
    BACK514(w,z); BACK514(z,y); BACK514(y,x); BACK514(x,w);
    _mm_store_si128(&p->m[1],x);
    _mm_store_si128(&p->m[2],y);
    _mm_store_si128(&p->m[3],z);
    _mm_store_si128(&p->m[4],w);
}

#else

 #define XOR514HEAD 16

void xor514sub(xor514_t*p) {
    enum { a=101,b=99,c=8 }; __m128i t=p->m[1],s=p->m[4];
    p->m[1]=p->m[2]; p->m[2]=p->m[3]; p->m[3]=s; 
    t=_mm_xor_si128(t,_mm_slli_si128(_mm_slli_epi64(t,a-64),8));
    t=_mm_xor_si128(t,_mm_srli_si128(_mm_srli_epi64(t,b-64),8));
    p->m[4]=_mm_xor_si128(s,_mm_xor_si128(_mm_srli_si128(s,c/8),t));
}

void back514sub(xor514_t*p) {
    enum { a=101,b=99,c=8 }; __m128i t=p->m[4],s=p->m[3];
    p->m[4]=s; p->m[3]=p->m[2]; p->m[2]=p->m[1];
    t=_mm_xor_si128(_mm_xor_si128(t,s),_mm_srli_si128(s,c/8));
    t=_mm_xor_si128(t,_mm_srli_epi64(_mm_srli_si128(t,8),b-64));
    p->m[1]=_mm_xor_si128(t,_mm_slli_epi64(_mm_slli_si128(t,8),a-64));
}

#endif

uint32_t xor514(xor514_t*p) {
    if (p->u[0]==20) { p->u[0]=XOR514HEAD; xor514sub(p); }
    return p->u[p->u[0]++];
}

uint32_t back514(xor514_t*p) {
    int i=p->u[0]-1; uint32_t r=p->u[i];
    if (i==XOR514HEAD) { i=20; back514sub(p); }
    p->u[0]=i; return r;
}


#define N 1000000000

int main(int argc,char**argv) {
    clock_t lap,min1,min2; uint32_t i,k,sum; xor514_t arg[1];

    if (argc==2) {
        uint32_t buf[1024];
        sxor514(arg,atoi(argv[1]));
        for (;;) {
            for (i=0;i<1024;i++) buf[i]=xor514(arg);
            fwrite(buf,4,1024,stdout);
        }
    }

    printf("sizeof(size_t)=%d\n",(int)sizeof(size_t));
    sxor514(arg,1);

    min1=0x7fffffff; sum=0; printf("xor128 ");
    for (k=0;k<8;k++) {
        lap=clock();
        for (i=0;i<N;i++) sum+=xor128();
        lap=(clock()-lap)*1000/CLOCKS_PER_SEC;
        printf(" %3lu",lap);
#if defined(__GNUC__)
        fflush(stdout);
#endif
        if (lap<min1) min1=lap;
    }
    printf("(%lu)%08x\n",min1,sum);

    min2=0x7fffffff; sum=0; printf("xor514 ");
    for (k=0;k<8;k++) {
        lap=clock();
        for (i=0;i<N;i++) sum+=xor514(arg);
        lap=(clock()-lap)*1000/CLOCKS_PER_SEC;
        printf(" %3lu",lap);
#if defined(__GNUC__)
        fflush(stdout);
#endif
        if (lap<min2) min2=lap;
    }
    printf("(%lu)%08x ",min2,sum);
    if (min2<min1) printf("%.2f%%fast\n",100.0*(min1-min2)/min1);
    else           printf("%.2f%%slow\n",100.0*(min2-min1)/min1);

    min2=0x7fffffff; printf("back514");
    for (k=0;k<8;k++) {
        lap=clock();
        for (i=0;i<N;i++) sum-=back514(arg);
        lap=(clock()-lap)*1000/CLOCKS_PER_SEC;
        printf(" %3lu",lap);
#if defined(__GNUC__)
        fflush(stdout);
#endif
        if (lap<min2) min2=lap;
    }
    printf("(%lu)%08x ",min2,sum);
    if (min2<min1) printf("%.2f%%fast\n",100.0*(min1-min2)/min1);
    else           printf("%.2f%%slow\n",100.0*(min2-min1)/min1);
    return 0;
}