［前：SplitMix64について］［ExRand目次］［次：SplitMix64とXorshiftを組み合わせる(大型)］

[ANC]SplitMix64とXorshiftを組み合わせる(小型)

SplitMix64 は周期が２⁶⁴と短くなっています。
そこで SplitMix64 と Xorshift を組み合わせた anc.h を開発しました。
これはメモリ使用量と速度をできるだけ増やさずに、周期を２²⁵⁶-２⁶⁴と長くしたものです。

多くのコンパイラでコンパイル可能です。

/* anc.h : coded by Isaku WADA(isaku@pb4.so-net.ne.jp) 2026
A PRNG combining SplitMix64 and Xorshift.

gcc, clang, bcc5.5.1, vc5.0 .. vc2026

*/
#if !defined(ANC_H__)
 #define ANC_H__
 #if defined(__BORLANDC__)
  #pragma warn -use /* 'identifier' declared but never used */
  #pragma warn -inl /* Functions containing * are not expanded inline */
 #endif
 #include <math.h>
 #include <stdlib.h>
 #if defined(_M_IX86) || defined(_MSC_VER)
typedef __int64 sqword;
typedef unsigned __int64 qword;
 #else
typedef long long sqword;
typedef unsigned long long qword;
 #endif
 #define GOLDEN_GAMMA 0x9E3779B97F4A7C15U
 #define GG_MUL_INV 0xF1DE83E19937733DU
 #define SPLITMIX_A 0xBF58476D1CE4E5B9U
 #define SPLITMIX_B 0x94D049BB133111EBU
 #define U01(A) ((sqword)((qword)(A)>>11)*(1.0/9007199254740992.0))
 #define AID 'A' 
 #define ANABC 3423902 /*35%*/
 #define ANSD (ANABC/1000000)

typedef struct Anc_t { qword xs[ANSD],cnt; } Anc_t;

static void AncNextState(qword*p) {
    enum { z=ANABC,h=100,a=z/h/h%h,b=z/h%h,c=z%h }; qword x=p[0],t;
 #if ANSD == 3
    t=p[1]; p[0]=t; t=p[2]; p[1]=t;
 #elif ANSD == 4
    t=p[1]; p[0]=t; t=p[2]; p[1]=t; t=p[3]; p[2]=t;
 #else
    int i;
    for (i=0;i<ANSD-1;i++) { t=p[i+1]; p[i]=t; }
 #endif
    t=x; x<<=a; t^=x; x=t; x>>=b; t^=x;
    x=p[ANSD-1]; t^=x; x>>=c; t^=x; p[ANSD-1]=t;
}

static void AncBackState(qword*p) {
    enum { z=ANABC,h=100,a=z/h/h%h,b=z/h%h,c=z%h }; qword x=p[ANSD-1],t;
 #if ANSD == 3
    t=p[1]; p[2]=t; t=p[0]; p[1]=t;
 #elif ANSD == 4
    t=p[2]; p[3]=t; t=p[1]; p[2]=t; t=p[0]; p[1]=t;
 #else
    int i;
    for (i=ANSD-1;i>0;i--) { t=p[i-1]; p[i]=t; }
 #endif
    t=p[ANSD-1]; x^=t; t>>=c; x^=t; t=x; x>>=b;
    t^=x; x=t; x<<=a; t^=x; p[0]=t;
}

static __inline qword AncSplitMixHash(qword z)
{ z=(z^(z>>30))*SPLITMIX_A; z=(z^(z>>27))*SPLITMIX_B; return z^(z>>31); }

static __inline qword AncRand(Anc_t*p) {
    qword z=(p->cnt+=GOLDEN_GAMMA);
    if ((short)z==0) AncNextState(p->xs);
    return AncSplitMixHash(z^p->xs[ANSD-1]);
}

static __inline qword AncBack(Anc_t*p) {
    qword z=AncSplitMixHash(p->cnt^p->xs[ANSD-1]);
    if ((short)p->cnt==0) AncBackState(p->xs);
    p->cnt-=GOLDEN_GAMMA; return z;
}

static void AncMulVecMat(qword*s,qword(*r)[ANSD]) {
    int i,j; qword t[ANSD];
    for (i=0;i<ANSD;i++) { t[i]=s[i]; s[i]=0; }
    for (i=0;i<64*ANSD;i++)
     if (t[i>>6]&((qword)1<<(i&63)))
      for (j=0;j<ANSD;j++) s[j]^=r[i][j];
}

static const char*AncJump(Anc_t*p,double lag) {
    enum { n=ANSD,m=64*n }; typedef struct mat { qword c[m][n]; } mat;
    typedef struct mat2 { mat a,b; } mat2; mat2*x=0; qword a=GOLDEN_GAMMA;
    int i,back=0,mod=(unsigned short)(p->cnt*GG_MUL_INV);
    if (lag<0) { lag=-lag; back=1; a=0-a; }
    i=(int)(modf(lag/65536,&lag)*65536+0.25);
    if (i) {
        if (back==0) { if (i>65536-mod) AncNextState(p->xs); }
        else         { if (i>mod) AncBackState(p->xs); }
        p->cnt+=a*i;
    }
    if (lag) {
 #if defined(__cplusplus)
        try { x=new mat2; } catch(...) { x=0; };
        if (x==0) return"new threw a exception";
 #else
        if ((x=(mat2*)malloc(sizeof*x))==0) return"malloc returned null";
 #endif
        for (i=0;i<m;i++) {
            int k; qword stat[n];
            for (k=0;k<n;k++) stat[k]=0;
            k=(i>>6)%n; stat[k]=(qword)1<<(i&63);
            if (back) AncBackState(stat); else AncNextState(stat);
            for (k=0;k<n;k++) x->a.c[i][k]=stat[k];
        }
    }
    a<<=16;
    for (;;) {
        if (modf(lag/2.0,&lag)!=0)
        { AncMulVecMat(p->xs,x->a.c); p->cnt+=a; }
        if (lag==0.0) break;
        x->b=x->a; a<<=1;
        for (i=0;i<m;i++) AncMulVecMat(x->a.c[i],x->b.c);
    }
 #if defined(__cplusplus)
    delete x;
 #else
    free(x);
 #endif
    return 0;
}

static void AncInitHelp(Anc_t*p);

static const char*AncInit(Anc_t*p,qword seed) {
    int i;
    for (i=0;i<ANSD;i++) p->xs[i]=AncSplitMixHash(seed+=GOLDEN_GAMMA);
    p->cnt=0; AncInitHelp(p); return 0;
}

static qword AncStrHash(qword r,const char*key) {
    if (key) for (;*key;key++) r=AncSplitMixHash(r^(unsigned char)*key);
    return r;
}

static const char*AncInitByStr(Anc_t*p,const char*key) {
    int i; qword seed=1;
    for (i=0;i<ANSD;i++) p->xs[i]=AncStrHash(seed+=GOLDEN_GAMMA,key);
    p->cnt=AncStrHash(seed+=GOLDEN_GAMMA,key); AncInitHelp(p); return 0;
}

static void AncInitHelp(Anc_t*p) {
 #if defined(__GNUC__)
    (void)AncInit; (void)AncInitByStr; (void)AncJump;
 #endif
    p->cnt&=~(qword)65535; p->cnt-=GOLDEN_GAMMA; AncRand(p);
}

#endif

#if !defined(ExRand)
 #define ExRand      AncRand
 #define ExBack      AncBack
 #define ExName     "AncRand"
 #define ExBName    "AncBack"
 #define ExInit      AncInit
 #define ExInitByStr AncInitByStr
 #define ExJump      AncJump
 #define ExSExp      0
 #define ExID        AID
 #define ExNsd       ANSD
 #define Ex_t        Anc_t
 #define ExAlign     8
#endif

AncRandをtestU01で検定する

QnxRand は AncRand をより一般化した疑似乱数です。
3A は QnxRand を AncRand と同じにするパラメータです。
まず testU01 で検定した結果の要約を示します。

type                    , c, test                      , p-value
QnxRand 3A rawhi Seed(1), B, Passed                    , -      
QnxRand 3A rawhi Seed(2), B, HammingIndep,L=30,r=27    , 5.0e-5 
QnxRand 3A rawhi Seed(3), B, Passed                    , -      
QnxRand 3A rawhi Seed(4), B, Passed                    , -      
QnxRand 3A rawhi Seed(5), B, Passed                    , -      
QnxRand 3A rawhi Seed(6), B, Passed                    , -      
QnxRand 3A rawlo Seed(1), B, Passed                    , -      
QnxRand 3A rawlo Seed(2), B, Passed                    , -      
QnxRand 3A rawlo Seed(3), B, Passed                    , -      
QnxRand 3A rawlo Seed(4), B, Passed                    , -      
QnxRand 3A rawlo Seed(5), B, LinearComp,r=29           , 6.7e-5 
QnxRand 3A rawlo Seed(6), B, Passed                    , -      
QnxRand 3A revhi Seed(1), B, ClosePairs NJumps,t=3     , 0.9991 
QnxRand 3A revhi Seed(2), B, Passed                    , -      
QnxRand 3A revhi Seed(3), B, Passed                    , -      
QnxRand 3A revhi Seed(4), B, Passed                    , -      
QnxRand 3A revhi Seed(5), B, Passed                    , -      
QnxRand 3A revhi Seed(6), B, Passed                    , -      
QnxRand 3A revlo Seed(1), B, Passed                    , -      
QnxRand 3A revlo Seed(2), B, CollisionOver,t=14        , 0.9994 
-                       , B, RandomWalk1 J(L=10000,r=0), 0.9992 
QnxRand 3A revlo Seed(3), B, MaxOft AD,t=16            , 7.9e-4 
QnxRand 3A revlo Seed(4), B, Passed                    , -      
QnxRand 3A revlo Seed(5), B, Passed                    , -      
QnxRand 3A revlo Seed(6), B, Passed                    , -      
QnxRand 3A swphi Seed(1), B, Passed                    , -      
QnxRand 3A swphi Seed(2), B, Passed                    , -      
QnxRand 3A swphi Seed(3), B, MatrixRank,L=30,r=26      , 0.9997 
QnxRand 3A swphi Seed(4), B, Passed                    , -      
QnxRand 3A swphi Seed(5), B, Passed                    , -      
QnxRand 3A swphi Seed(6), B, Passed                    , -      
QnxRand 3A swplo Seed(1), B, Passed                    , -      
QnxRand 3A swplo Seed(2), B, Passed                    , -      
QnxRand 3A swplo Seed(3), B, Passed                    , -      
QnxRand 3A swplo Seed(4), B, Passed                    , -      
QnxRand 3A swplo Seed(5), B, Passed                    , -      
QnxRand 3A swplo Seed(6), B, Passed                    , -

このように rdrand とほぼ同じ品質であることがわかります。

AncRandをPractRand ver0.96で検定する

次に PractRand ver0.96 での検定です。
以下の６つの検定をしました。

(1)RNG_test stdin64 -tlmin 1GB -multithreaded
(2)RNG_test stdin64 -tlmin 1GB -tf 2 -multithreaded
(3)RNG_test stdin64 -tlmin 1GB -te 1 -tlmax 8TB -multithreaded
(4)RNG_test stdin64 -tlmin 1GB -te 1 -tf 2 -tlmax 8TB -multithreaded
(5)RNG_test stdin64 -tlmin 1GB -te 10 -tlmax 8TB -multithreaded
(6)RNG_test stdin64 -tlmin 1GB -te 10 -tf 2 -tlmax 8TB -multithreaded

本来は -te 1 や -te 10 でも既定の 32TB まで検定したかったのですが、
-te 0 に比べて時間がかかるので 8TB までにしました。
(1),(3),(4),(5),(6)は、種１で異常はありませんでした。
(2)RNG_test stdin64 -tlmin 1GB -tf 2 -multithreadedは、種１で、

RNG_test stdin64 -tlmin 1GB -tf 2 -multithreaded
QnxRand te0 tf2 3A Seed(1)
RNG_test using PractRand version 0.96
RNG = RNG_stdin64, seed = unknown
test set = core, folding = extra
...
length= 32 terabytes (2^45 bytes), time= 143463 seconds
  Test Name                         Raw       Processed     Evaluation       
  [Low4/16]BCFN(2+1,13-0U)          R=  +9.3  p =  1.5e-4   unusual

となりましたが、種２で、1GB から 32TB まで異常はありませんでした。

PractRand ver0.96で QncRand を種を1～40にして 1GB, 2GB, 4GB の3つの長さで、異常のあったものを表示します。
この120回分を1つのテストとします。
-te 0 を te0、-te 1 を te1、-tf 1 を tf1、-tf 2 を te2 と表します。
１つのテストを(te0,te1)と(tf1,tf2)の4個の組み合わせをテストして合計480回調べました。

AncRand 3A te0 tf1 Seed(7) 4GB  [Low16/64]BCFN(2+3,13-1U) p=2.1e-4 unusual 
AncRand 3A te0 tf1 Seed(8) 4GB  DC6-9x1Bytes-1 p=1-1.5e-3 unusual 
AncRand 3A te0 tf1 Seed(17) 4GB  [Low1/64]DC6-9x1Bytes-1 p=1-2.8e-3 unusual 
AncRand 3A te0 tf1 Seed(23) 4GB  [Low4/64]BCFN(2+0,13-2U) p=1-8.7e-4 unusual 
AncRand 3A te0 tf1 Seed(28) 2GB  FPF/16:all p=1-9.2e-5 unusual 
AncRand 3A te0 tf1 Seed(30) 1GB  [Low4/64]DC6-9x1Bytes-1 p=1-3.6e-4 mildly suspicious
AncRand 3A te0 tf1 Seed(31) 2GB  DC6-9x1Bytes-1 p=1.3e-3 unusual 
AncRand 3A te0 tf1 Seed(31) 4GB  DC6-9x1Bytes-1 p=2.0e-4 mildly suspicious
AncRand 3A te0 tf1 Seed(35) 2GB  DC6-9x1Bytes-1 p=1.8e-3 unusual 
AncRand 3A te0 tf1 Seed(37) 1GB  [Low16/64]DC6-9x1Bytes-1 p=1.6e-3 unusual

AncRand 3A te0 tf2 Seed(6) 2GB  [Low4/16]Gap-16:A p=1-3.0e-4 unusual 
AncRand 3A te0 tf2 Seed(8) 4GB  DC6-9x1Bytes-1 p=1-1.5e-3 unusual 
AncRand 3A te0 tf2 Seed(11) 4GB  [Low1/16]DC6-9x1Bytes-1 p=1.2e-3 unusual 
AncRand 3A te0 tf2 Seed(16) 1GB  [Low8/64]FPF/16:all p=1.8e-4 unusual 
AncRand 3A te0 tf2 Seed(28) 2GB  FPF/16:all p=1-9.2e-5 unusual 
AncRand 3A te0 tf2 Seed(30) 1GB  [Low4/64]DC6-9x1Bytes-1 p=1-3.6e-4 unusual 
AncRand 3A te0 tf2 Seed(30) 2GB  [Low1/16]DC6-9x1Bytes-1 p=1.0e-3 unusual 
AncRand 3A te0 tf2 Seed(31) 1GB  [Low4/32]DC6-9x1Bytes-1 p=1-1.3e-4 mildly suspicious
AncRand 3A te0 tf2 Seed(31) 2GB  DC6-9x1Bytes-1 p=1.3e-3 unusual 
AncRand 3A te0 tf2 Seed(31) 4GB  DC6-9x1Bytes-1 p=2.0e-4 unusual

AncRand 3A te1 tf1 Seed(12) 4GB  [Low1/64][C8]BCFN(2+5,13-8U) p=6.1e-6 unusual 
AncRand 3A te1 tf1 Seed(14) 1GB  [Low1/64]BCFN(0+1,13-4U) p=9.5e-5 unusual 
AncRand 3A te1 tf1 Seed(18) 1GB  [Low16/64]FPF/64:all p=3.2e-5 unusual 
AncRand 3A te1 tf1 Seed(26) 2GB  NS3[3:hw:both] p=8.2e-5 unusual 
AncRand 3A te1 tf1 Seed(27) 4GB  DC7-6x2Bytes-1:indep p=1-1e-4 unusual 
AncRand 3A te1 tf1 Seed(29) 4GB  NS3[6:hw:both] p=5.1e-5 unusual 
AncRand 3A te1 tf1 Seed(30) 1GB  [Low4/64]DC7-9x1Bytes-1:both p=1-3.6e-4 unusual 
AncRand 3A te1 tf1 Seed(31) 4GB  DC7-9x1Bytes-1:both p=2.0e-4 unusual 
AncRand 3A te1 tf1 Seed(40) 1GB  [Low1/64]DC7-6x2Bytes-1:indep p=1-1e-4 unusual

AncRand 3A te1 tf2 Seed(10) 2GB  [Low1/32]NS3[4:hw:both] p=5.2e-6 unusual 
AncRand 3A te1 tf2 Seed(31) 1GB  [Low4/32]DC7-9x1Bytes-1:both p=1-1.3e-4 unusual 
AncRand 3A te1 tf2 Seed(31) 4GB  DC7-9x1Bytes-1:both p=2.0e-4 unusual

te0 で mildly suspicious が出現しますが、p-value を見ると極端ではありません。
また、rdrandをtestU01とPractRand ver0.96で検定するの結論より問題ありません。
従って、rdrand とほぼ同じ品質であることがわかります。

ANC の SSE や AVX のプログラムは、周期の長さが変わらないのに多くのメモリを必要とするので用意しませんでした。

［前：SplitMix64について］［ExRand目次］［次：SplitMix64とXorshiftを組み合わせる(大型)］