/*
jump.c : Xorshift で戻ったり一瞬で相対ジャンプできるプログラム

                         coded by isaku@pb4.so-net.ne.jp  2008/11

void   InitXor(ULONG x[4],ULONG seed) 種で初期化
double NextUnif(ULONG x[4])           [0,1)乱数の生成
double BackUnif(ULONG x[4])           前の乱数を返す
void   JumpXor(ULONG dst[4],ULONG src[4],double interval)
                      src からジャンプしたものを dst に転送する

複数の状態ベクトルを用意して、互いに1000兆個以上間隔をあければ、
まずオーバーラップの心配なく複数系列のシミュレーションを行える。
ジャンプする距離は倍精度浮動小数点数で表せる範囲内で指定できる。

周期は、340正2823潤6692溝0938穣4634予6337京4607兆4317億6821万1455
だけあるので、同じパターンが繰り返す心配はまずない。

128 bit Xorshift は連続する4つの出力を得れば、
将来だけでなく過去の数列も分かってしまうので、暗号には向かない。

BackXor() は前のバージョンでは、最初に呼び出したときに、
一つ前の乱数を返していたが、このバージョンでは最後の乱数を返す。
そのほうが扱いやすい。

マクロ FAST_UNIF を定義してコンパイルすると、NextUnif の
精度が 53 bit から 32 bit に下がるが少し高速になる。

ある距離だけジャンプすると最初は表を作るので時間がかかるが、
2度目以降に同じ距離だけジャンプするときはキャッシュに表が残って
いるので、瞬間的にジャンプできる。キャッシュの大きさを大きくすれば
それだけ多種類の距離の表を残すことができる。キャッシュを使い切ると、
古いキャッシュを消して新しい距離の表を記録する。
表一つ当たり2Kバイトのメモリを消費する。作業用に2つの表を
使っているので、キャッシュの個数が10ならば、全体で22Kバイトの
メモリを消費している。
*/

#include <string.h>
#include <math.h>

/*┌──────────────────────────────┐
 │128 ビットXorshiftのパラメータ(11,8,19)                     │
 │論文中の他の組合せ(5,14,1),(15,4,21),(23,24,3),(5,12,29)    │
 │他の組合せ(1,3,12),(1,3,15),(2,1,21),(2,21,6),(2,25,2)      │
 │(3,2,21),(4,1,5),(6,5,17),(6,11,21),(6,11,25),(7,10,7),     │
 │(7,11,19),(7,11,20),(7,12,11),(8,11,14),(9,11,6),(9,13,17)  │
 │(9,24,1),(10,5,8),(10,11,12),(10,11,23),(11,5,12),(11,5,24) │
 │(11,5,26),(11,10,21),(11,16,1),(13,3,25),(14,3,23)          │
 │(14,13,19)),(17,7,21),(18,13,19),(19,1,2),(20,5,17)         │
 │(21,2,23),(21,9,4),(21,16,11),(22,3,12),(23,3,6),(25,3,10)  │
 │(27,5,31),(27,19,5),(29,3,30)                               │
 └──────────────────────────────┘*/
enum { A=11, B=8, C=19 };

/*┌─────────────────────────────┐
 │キャッシュの大きさ(1増えると2Kバイト必要メモリが増える)│
 └─────────────────────────────┘*/
#define JUMP_CACHE 10

/* #define FAST_UNIF /**/

#ifndef ULONG
#define ULONG unsigned long
#endif

/*┌───────────────┐
 │状態ベクトルを種sで初期化する│
 └───────────────┘*/
void InitXor(ULONG x[4],ULONG s)
{ int i; for (i=1;i<=4;i++) x[i-1]=s=1812433253L*(s^(s>>30))+i; }

/*┌────────────────┐
 │状態ベクトルを指定できるXorshift│
 └────────────────┘*/
ULONG NextXor(ULONG x[4]) {
    ULONG t=x[0],r=t; t<<=A; t^=r; r=t; r>>=B; t^=r; r=x[1]; x[0]=r;
    r=x[2]; x[1]=r; r=x[3]; x[2]=r; t^=r; r>>=C; r^=t; x[3]=r;
    return r;
}

/*┌─────────────────────┐
 │t^=t<<A の逆変換が InvLeft(t,A) で得られる│
 └─────────────────────┘*/
#define InvLeft(T,A) { /* ULONG s,r; */ s=(1UL<<A)-1; \
 do { s<<=A; r=t; r<<=A; r&=s; t^=r; } while ((s&1UL<<31)==0); }

/*┌──────────────────────┐
 │t^=t>>A の逆変換が InvRight(t,A) で得られる │
 └──────────────────────┘*/
#define InvRight(T,A) { /* ULONG s,r; */ s=(1UL<<A)-1<<32-A; \
 do { s>>=A; r=t; r>>=A; r&=s; t^=r; } while ((s&1UL)==0); }

/*┌───────────┐
 │出力順序が逆なXorshift│
 └───────────┘*/
ULONG BackXor(ULONG x[4]) {
    ULONG u=x[3],r=x[2],t=u,s; x[3]=r; t^=r; r>>=C; t^=r; r=x[1];
    x[2]=r; r=x[0]; x[1]=r; InvRight(t,B); InvLeft(t,A); x[0]=t;
    return u;
}

/*┌───────────────────────┐
 │0以上1未満の浮動小数点一様乱数(32/53bit精度)│
 └───────────────────────┘*/
double NextUnif(ULONG x[4]) {
#ifdef FAST_UNIF
    return NextXor(x)*(1.0/4294967296.0);
#else
    ULONG a=NextXor(x)>>11,b=NextXor(x);
    return(b*2097152.0+a)*(1.0/9007199254740992.0);
#endif
}

/*┌──────────────────────────┐
 │0以上1未満の浮動小数点一様乱数(32/53bit精度)(逆順)│
 └──────────────────────────┘*/
double BackUnif(ULONG x[4]) {
#ifdef FAST_UNIF
    return BackXor(x)*(1.0/4294967296.0);
#else
    ULONG b=BackXor(x),a=BackXor(x)>>11;
    return(b*2097152.0+a)*(1.0/9007199254740992.0);
#endif
}

/*┌─────────────────┐
 │相対ジャンプをするための下請け関数│
 └─────────────────┘*/
void PowerJumpXor(ULONG p[4],ULONG r[128][4]) {
    int i; ULONG q[4],*y;
    q[0]=p[0];q[1]=p[1];q[2]=p[2];q[3]=p[3]; p[3]=p[2]=p[1]=p[0]=0;
    for (i=0;i<128;i++) if (q[i>>5]&(1L<<(i&31)))
    { y=r[i]; p[0]^=y[0]; p[1]^=y[1]; p[2]^=y[2]; p[3]^=y[3]; }
}

/*┌─────────────────┐
 │高速化のためのテーブルを初期化する│
 └─────────────────┘*/
void SetJumpXor(ULONG r[128][4],double u) {
    static ULONG work[128][4],next[128][4]; int i,k; double v;
    if (u<0) { u=-u; k=1; } else k=0;
    memset(work,0,sizeof work);
    for (i=0;i<128;i++) {
        work[i][i>>5]=1L<<(i&31);
        if (k) BackUnif(work[i]); else NextUnif(work[i]);
    }
    u=floor(u); u/=2; v=floor(u);
    if (u!=v) memcpy(r,work,sizeof work);
    else {
        memset(r,0,sizeof work);
        for (i=0;i<128;i++) r[i][i>>5]=1L<<(i&31);
    }
    memcpy(next,work,sizeof work);
    for (u=v;;u=v) {
        u/=2; v=floor(u);
        for (i=0;i<128;i++) PowerJumpXor(next[i],work);
        if (u!=v) for (i=0;i<128;i++) PowerJumpXor(r[i],next);
        if (v==0.0) break;
        memcpy(work,next,sizeof work);
    }
}

/*┌──────────────────────┐
 │srcからuだけジャンプした状態をdstに転送する │
 │負の値を与えるとバックする          │
 │キャッシュの管理も行う              │
 └──────────────────────┘*/
void JumpXor(ULONG dst[4],ULONG src[4],double u) {
    static double b[JUMP_CACHE]; static int c;
    static ULONG tbl[JUMP_CACHE][128][4]; int i;
    if (u==0.0) return;
    for (i=0;i<JUMP_CACHE&&b[i]!=0.0;i++) if (u==b[i]) break;
    if (i==JUMP_CACHE) {
        i=c++; b[i]=u; SetJumpXor(tbl[i],u);
        if (c==JUMP_CACHE) c=0;
    }
    if (b[i]==0.0) { b[i]=u; SetJumpXor(tbl[i],u); }
    if (dst!=src) memcpy(dst,src,sizeof*dst*4);
    PowerJumpXor(dst,tbl[i]);
}

/*────────────────────────────────*/
/* 以上、乱数のプログラム終了。以下動作確認のためのプログラム。 */
/*────────────────────────────────*/

#if 1

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

/*┌────────────────────┐
 │MS-DOS コンパイラ用のシンプルな clock() │
 └────────────────────┘*/
#if CLOCKS_PER_SEC == 1
#undef CLOCKS_PER_SEC
#endif
#ifndef CLOCKS_PER_SEC
#define CLOCKS_PER_SEC 1000
#include <dos.h>
static long clock(void) {
    union REGS t; t.h.ah=0x2c; intdos(&t,&t);
    return t.h.dl*10L+t.h.dh*1000L+t.h.cl*60000L+t.h.ch*3600000L;
}
#endif

/*┌─────────┐
 │Xorshiftの初期状態│
 └─────────┘*/
#define DEF_X 1812433254UL
#define DEF_Y 3713160357UL
#define DEF_Z 3109174145UL
#define DEF_W   64984499UL

/*┌────────────────┐
 │検証のためのオリジナルのXorshift│
 └────────────────┘*/
ULONG xor128(void) {
    static ULONG x=DEF_X,y=DEF_Y,z=DEF_Z,w=DEF_W;
    ULONG t=x^(x<<A); x=y; y=z; z=w; return w^=(w>>C)^t^(t>>B);
}

/*┌───────────────────────┐
 │0以上1未満の浮動小数点一様乱数(32/53bit精度)│
 └───────────────────────┘*/
double xor_unif(void) {
#ifdef FAST_UNIF
    return xor128()*(1.0/4294967296.0);
#else
    ULONG a=xor128()>>11,b=xor128();
    return(b*2097152.0+a)*(1.0/9007199254740992.0);
#endif
}

/*┌─────────────┐
 │実験で生成させる乱数の個数│
 └─────────────┘*/
#define N 16777216L

/*┌───────┐
 │メインルーチン│
 └───────┘*/
int main(void) {
    ULONG i,lap; static ULONG x[4]={DEF_X,DEF_Y,DEF_Z,DEF_W};
    double period=2.0;

    printf("xor128()とNextXor()の最初の4個の出力\n");
    for (i=0;i<4;i++)
        printf("xor()%.17f next()%.17f\n",xor_unif(),NextUnif(x));

    lap=clock();
    for (i=0;i<N;i++) xor_unif();
    printf("xor_unif()で%luの生成時間:%lumsec\n",N,clock()-lap);

    lap=clock(); JumpXor(x,x,(double)N);
    printf("JumpXor()で%luの移動時間:%lumsec\n",N,clock()-lap);

    printf("ジャンプできたかどうか確認する\n");
    for (i=0;i<4;i++)
        printf("xor()%.17f next()%.17f\n",xor_unif(),NextUnif(x));

    lap=clock(); JumpXor(x,x,-N-4.0);
    printf("JumpXor()で%luの逆移動時間:%lumsec\n",N+4,clock()-lap);

    printf("逆順に出力して元に戻ったかどうかの確認\n");
    for (i=0;i<4;i++) printf("back()%.17f\n",BackUnif(x));

    lap=clock();
    for (i=0;i<7;i++) period*=period;
    JumpXor(x,x,period); BackUnif(x);
    printf("JumpXor()で%eの移動時間:%lumsec\n",period,clock()-lap);
    printf("1周期ジャンプして同じものか繰り返されるか\n");
    for (i=0;i<4;i++) printf("next()%.17f\n",NextUnif(x));
    return 0;
}

#endif