/* dxor513.c : coded by isaku@pb4.so-net.ne.jp
An Implementation of Xorshift RNG using SSE2(Double version)

       codes by Isaku WADA  Jun 2014

One billion generation time(sec):

Core i7 2600(turbo boost off)
+-----------------------------------------------+-------+-------+--------+
|compiler and system                            |dxor128|dxor513|dback513|
+-----------------------------------------------+-------+-------+--------+
|VC++2010 64bit cl /O2 /GL /W4 on Windows7 64bit|  2.121|  2.168|   1.887|
|VC++2010 32bit cl /O2 /GL /W4 on Windows7 64bit|  4.602|  1.840|   1.856|
|gcc -O4 -Wall             on Ubuntu 14.04 64bit|  2.655|  1.976|   1.863|
|gcc -O4 -Wall -m32 -msse2 on Ubuntu 14.04 64bit|  2.430|  1.847|   1.817|
+-----------------------------------------------+-------+-------+--------+

dxor513() is based on xor514.c.
xor514() does not FAIL dieharder test(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;

 This dxor513() provide 2^{513}-2 random 52bit pricision double
 with range [0,1).

*/

#if defined(_MSC_VER)
 #pragma warning(disable:4514)
 #pragma warning(disable:4710)
 #pragma warning(disable:4711)
 #pragma warning(disable:4820)
#endif

#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
*/

typedef struct dTag128 { uint32_t x,y,z,w; } dxor128_t;

double dNext128(dxor128_t*p) {
    enum { a=11,b=8,c=19 };
    uint32_t t=p->x ^ p->x<<a;  p->x=p->y; p->y=p->z; p->z=p->w;
    return (p->w ^= p->w>>c ^ t ^ t>>b)*(1.0/4294967296.0);
}

void dSrand128(dxor128_t*p,uint32_t s) {
   s^=521288629;
   p->x=s=1812433253*(s^(s>>30))+1;
   p->y=s=1812433253*(s^(s>>30))+2;
   p->z=s=1812433253*(s^(s>>30))+3;
   p->w=s=1812433253*(s^(s>>30))+4;
}

#if defined(__GNUC__)
 #define MUL513 1
#else
 #define MUL513 16 /* Expand status table */
#endif

typedef union dTag513 {
    uint32_t u[4+MUL513*32];
    double   d[2+MUL513*16];
    __m128i  m[1+MUL513*8];
    __m128d  r[1+MUL513*8];
    int      i;
} dxor513_t;

typedef union dTagM { __m128i i; __m128d r; } xmm_t;

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

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

#if RND_SSE2_NOINLINE && defined(__GNUC__)
extern void Next513Sub(dxor513_t*p) __attribute__((noinline));
extern void Back513Sub(dxor513_t*p) __attribute__((noinline));
#endif

#define CV513(X,N)                                       \
 t.i=_mm_srli_epi64(X.i,12); t.i=_mm_or_si128(t.i,e.i);  \
 t.r=_mm_sub_pd(t.r,o.r);    _mm_store_si128(&p->m[N],t.i)

void Next513Sub(dxor513_t*p) {
    const union { short  s[8]; __m128i i; } ee={{0,0,0,0x3ff0,0,0,0,0x3ff0}};
    const union { double d[2]; __m128i i; } oo={{1.0,1.0}};
#if MUL513 != 1
    int i;
#endif
    enum { a=101,b=99,c=8,m=MUL513*4-4 }; xmm_t t,x,y,z,w,e,o;
    e.i=_mm_load_si128(&ee.i); o.i=_mm_load_si128(&oo.i);
    x.i=_mm_load_si128(&p->m[m+1]); y.i=_mm_load_si128(&p->m[m+2]);
    z.i=_mm_load_si128(&p->m[m+3]); w.i=_mm_load_si128(&p->m[m+4]);
#if MUL513 == 1
    NEXT513(x,w); NEXT513(y,x);  NEXT513(z,y); NEXT513(w,z);
    _mm_store_si128(&p->m[1],x.i); _mm_store_si128(&p->m[2],y.i);
    _mm_store_si128(&p->m[3],z.i); _mm_store_si128(&p->m[4],w.i);
    CV513(x,5); CV513(y,6); CV513(z,7); CV513(w,8);
#else
    for (i=0;i<MUL513*4;i+=4) {
        NEXT513(x,w); NEXT513(y,x); NEXT513(z,y); NEXT513(w,z);
        _mm_store_si128(&p->m[1+i],x.i); _mm_store_si128(&p->m[2+i],y.i);
        _mm_store_si128(&p->m[3+i],z.i); _mm_store_si128(&p->m[4+i],w.i);
        CV513(x,m+5+i); CV513(y,m+6+i); CV513(z,m+7+i); CV513(w,m+8+i);
    }
#endif
}

void Back513Sub(dxor513_t*p) {
    const union { short  s[8]; __m128i i; } ee={{0,0,0,0x3ff0,0,0,0,0x3ff0}};
    const union { double d[2]; __m128i i; } oo={{1.0,1.0}};
#if MUL513 != 1
    int i;
#endif
    enum { a=101,b=99,c=8,m=MUL513*4-4 }; xmm_t t,x,y,z,w,e,o;
    e.i=_mm_load_si128(&ee.i); o.i=_mm_load_si128(&oo.i);
    x.i=_mm_load_si128(&p->m[1]); y.i=_mm_load_si128(&p->m[2]);
    z.i=_mm_load_si128(&p->m[3]); w.i=_mm_load_si128(&p->m[4]);
#if MUL513 == 1
    BACK513(w,z); BACK513(z,y); BACK513(y,x); BACK513(x,w);
    _mm_store_si128(&p->m[4],w.i); _mm_store_si128(&p->m[3],z.i);
    _mm_store_si128(&p->m[2],y.i); _mm_store_si128(&p->m[1],x.i);
    CV513(x,5); CV513(y,6); CV513(z,7); CV513(w,8);
#else
    for (i=MUL513*4;i>0;) {
        i-=4; BACK513(w,z); BACK513(z,y); BACK513(y,x); BACK513(x,w);
        _mm_store_si128(&p->m[4+i],w.i); _mm_store_si128(&p->m[3+i],z.i);
        _mm_store_si128(&p->m[2+i],y.i); _mm_store_si128(&p->m[1+i],x.i);
        CV513(w,m+8+i); CV513(z,m+7+i); CV513(y,m+6+i); CV513(x,m+5+i);
    }
#endif
}

double dNext513(dxor513_t*p) {
    if (p->i==2+MUL513*16) { p->i=2+MUL513*8; Next513Sub(p); }
    return p->d[p->i++];
}

double dBack513(dxor513_t*p) {
    double r=p->d[--p->i];
    if (p->i==2+MUL513*8) { p->i=2+MUL513*16; Back513Sub(p); }
    return r;
}

void dSrand513(dxor513_t*p,uint32_t s) {
   uint32_t i; p->i=2+MUL513*16; s^=521288629;
   for (i=0;i<16;i++) p->u[i+MUL513*16-12]=s=1812433253*(s^(s>>30))+i;
   dNext513(p); dBack513(p);
}

#define N 1000000000

int main(void) {
    clock_t lap,min1,min2; uint32_t i,k; double sum;
#if defined(__GNUC__)
    dxor128_t tbl128[1]={{123456789,362436069,521288629,88675123}};
    dxor513_t tbl513[1];
#else
    static dxor128_t tbl128[1]={{123456789,362436069,521288629,88675123}};
    static dxor513_t tbl513[1];
#endif

    printf("sizeof(size_t)=%d\n",(int)sizeof(size_t));
    dSrand513(tbl513,1);

    {
        double n[100],b[100];
        for (i=0;i<100;i++) n[i]=dNext513(tbl513);
        while (i--) b[i]=dBack513(tbl513);
        for (i=0;i<100;i++) if (n[i]!=b[i]) printf("*");
    }

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

    min2=0x7fffffff; sum=0; printf("dNext513");
    for (k=0;k<6;k++) {
        lap=clock();
        for (i=0;i<N;i++) sum+=dNext513(tbl513)-0.5;
        lap=(clock()-lap)*1000/CLOCKS_PER_SEC;
        printf(" %3lu",lap);
#if defined(__GNUC__)
        fflush(stdout);
#endif
        if (lap<min2) min2=lap;
    }
    printf("(%lu)%g ",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("dBack513");
    for (k=0;k<6;k++) {
        lap=clock();
        for (i=0;i<N;i++) sum-=dBack513(tbl513)-0.5;
        lap=(clock()-lap)*1000/CLOCKS_PER_SEC;
        printf(" %3lu",lap);
#if defined(__GNUC__)
        fflush(stdout);
#endif
        if (lap<min2) min2=lap;
    }
    printf("(%lu)%g ",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;
}