#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
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
#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;
}