Wywaliłem diehardera. Napisałem na kolanie kilka swoich
testów, pewnie moje są gorsze, ale przynajmniej mam jasność
co do zasady działania i wiem kiedy ryzykuję przekroczenie
wartości parametrów.

Pierwszy wniosek: Do tej pory nigdy nie zaobserwowałem aby MT
oblał stabilny test. Odpaliłem kilka testów dla ranluxa,
na oko bylo widać że ranlux zalicza testy lepiej. Odpaliłem
dłuższy test, jak ktoś ma ochotę, niech odpali u siebie.

Wyniki dla std::mt19937_64
test1 24.592034%
test1 1.487763%
test1 4.233741%
test1 26.478715%
test1 4.174345%
test1 27.751500%
test1 18.710311%
test1 27.819010%
test1 39.071225%
test2 23.595284%
test2 48.625664%
test3 45.054474%
test3 0.498713%
test4 37.237054%
test4 5.056990%
test4 11.755519%
test4 45.547699%
test5 15.588600%
test5 20.128564%
test5 9.747300%
test5 41.895282%
test5 25.154676%
test5 17.381929%
test5 20.064510%
test6 2.690305%
test6 14.040318%
test6 15.422749%
test6 17.341370%
test6 7.701529%
test6 6.545549%
test6 23.770691%
test6 18.814217%
test6 8.335626%
test6 0.035693%

Ostatni wynik niby jest niski, ale paradoks dnia urodzin dla 100 dni
wymaga bardzo dużo liczb losowych, a po drugie ten kod może już
tracić stabilność numeryczną dla 99 punktów swobody. No i jeszcze
mogą być jakieś błędy w kodzie, zobaczymy jak długi test
przejdzie ranlux, gdy przejdzie źle, to bedzie wskazówka do
szukania błędu.






#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <ctime>
#include <cstring>
#include <random>

typedef double ftyp;
typedef const ftyp cftyp;

typedef __float128 dtyp;
typedef const dtyp cdtyp;


typedef int ityp;
typedef const ityp cityp;

typedef unsigned int utyp;
typedef const utyp cutyp;

typedef long long ltyp;
typedef const long long cltyp;

typedef unsigned long long ultyp;
typedef const unsigned long long cultyp;

#define CHI_STEPS (512)

//Przykłowy generator
class RLin {
private:
ultyp x;
public:
RLin() { sRand(12347512); }
RLin( cultyp seed ) { sRand(seed); }
void sRand( cultyp seed ) { x = seed; }
ultyp getUL() { return x = x * 25214903917ULL + 11ULL; }
utyp getU() { return (utyp)getUL(); }
utyp operator()() { return getU(); }
ftyp getF() { return (ftyp) ( (getUL() & ((1ull<<40)-1)) / (ftyp)(1ull<<40)); }
};


// Prawie liniowy generator ;-)
class RLin2 {
private:
ultyp a,b,c,d;
utyp n;
static cultyp m[20];
public:
RLin2( cultyp seed=12347512 ) { sRand(seed); }
void sRand( ultyp seed ) {
a = (seed = seed * 25214903917ULL + 11ULL);
b = (seed = seed * 25214903917ULL + 11ULL);
c = (seed = seed * 25214903917ULL + 11ULL);
d = (seed = seed * 25214903917ULL + 11ULL);
n = 0xFFFFFFFFu;
}
ultyp getUL() {
if( ++n >= 20 ) n = 0;
cultyp r = (a * 0xEC5ECAB872DF70B2ull + b * 0x36285E818B75145Full + c *
0xEA8DB7BF70F68277ull + d * 0x32592BEF5EC46BA8ull + 0x1D1BF584DBC7A169ull) % m[n];
d = c;
c = b;
b = a;
return a = r;
}
utyp getU() { return (utyp)getUL(); }
utyp operator()() { return getU(); }
ftyp getF() { return (ftyp) ( (getUL() & ((1ull<<40)-1)) / (ftyp)(1ull<<40)); }
};

cultyp RLin2::m[] = {
0x3EDB94BD5CE7041Eull,
0xCF64D8021AFBE487ull,
0xD7BE02752CCAD058ull,
0xC258C6DFB3B92DD8ull,
0x6DA82BC9D51F561Eull,
0x0CF26705AFC7640Eull,
0xAC25BA196CCC5A2Aull,
0x7E52DA42E37C8323ull,
0xBD72B519BBF6F08Aull,
0x68B84AF1BBA868A6ull,
0x8BC55D980AC0C053ull,
0x639366BAE2A9815Bull,
0x4E0BECBDECFFC92Dull,
0x6FCE8E25E5C82797ull,
0xDBD42E002C7B6D6Dull,
0xBDBFFA31AA65CDE1ull,
0xFE69A7DA2830FA95ull,
0xEA119C87B81B3C45ull,
0x818F544C89F8F908ull,
0x3C968AEF7ED95B40ull,
};


class RMT {
private:
std::mt19937_64 re;
public:
RMT( cultyp seed=12347512 ) { sRand(seed); }
void sRand( cultyp seed ) { re.seed(seed); }
ultyp getUL() { return re(); }
utyp getU() { return (utyp)getUL(); }
utyp operator()() { return getU(); }
ftyp getF() { return (ftyp) ( (getUL() & ((1ull<<40)-1)) / (ftyp)(1ull<<40)); }
};


class RMTSu {
private:
std::shuffle_order_engine< std::mt19937_64 , 20000 > re;
public:
RMTSu( cultyp seed=12347512 ) { sRand(seed); }
void sRand( cultyp seed ) { re.seed(seed); }
ultyp getUL() { return re(); }
utyp getU() { return (utyp)getUL(); }
utyp operator()() { return getU(); }
ftyp getF() { return (ftyp) ( (getUL() & ((1ull<<40)-1)) / (ftyp)(1ull<<40)); }
};



class RRnx {
private:
std::ranlux48 re;
public:
RRnx( cultyp seed=12347512 ) { sRand(seed); }
void sRand( cultyp seed ) { re.seed(seed); }
ultyp getUL() { return re(); }
utyp getU() { return (utyp)getUL(); }
utyp operator()() { return getU(); }
ftyp getF() { return (ftyp) ( (getUL() & ((1ull<<40)-1)) / (ftyp)(1ull<<40)); }
};


//testowany generator
typedef RRnx RndTest;

static ftyp f2( cftyp x ) { return x*x; }

// rozkład x^2
static dtyp dCh2( cdtyp x , cdtyp k ) {
return pow( 0.5 , 0.5*k ) * pow( x , 0.5*k-1 ) * exp( -0.5*x ) / tgamma( 0.5*k
);
}

// dystrybuanta x^2 numeryczna
static ftyp ddCh2( cdtyp x , cdtyp k ) {
dtyp step = x / CHI_STEPS;
if( step == 0 )
return 0;
dtyp sum=0;
dtyp old, v = 0;
for( dtyp i=0 ; i<=x ; i+=step ) {
old = v;
v = dCh2( i , k );
if( isinf(v) ) v = 0;
sum += (v + old) * 0.5 * step; // trapezami
}
return (ftyp)sum;
// return (ftyp)(sum < 1-sum ? sum : 1-sum) ; // przedział lewy lub prawy
}


// N-kubełków wypełniamy równomiernie
static ftyp test1( RndTest &rnd , cultyp SIZE, cutyp BUCKET, cutyp SKIP=0 ) {
ultyp *const freq = new ultyp[BUCKET];
memset( freq , 0 , sizeof(ultyp) * BUCKET );

for( ultyp i=0 ; i<SIZE ; i++ ) {
freq[ rnd() % BUCKET ] ++;
for( utyp j=0 ; j<SKIP ; j++ ) rnd();
}
ftyp sum = 0;
for( utyp i=0 ; i<BUCKET ; i++ ) {
// printf("%llu,",freq[i]);
sum += f2(freq[i] - (ftyp)SIZE/BUCKET) / ((ftyp)SIZE/BUCKET);
}
// printf("\n");
delete[] freq;
return ddCh2( sum , BUCKET-1 );
}


// całka wielomianu piątego stopnia
static ftyp test2( RndTest &rnd , cultyp SIZE, cutyp SKIP1=0 , cutyp SKIP2=0 ) {
ultyp ok = 0;
for( ultyp i=0 ; i<SIZE ; i++ ) {
cftyp x = rnd.getF() * 2.0 - 1.0;
for( utyp j=0 ; j<SKIP1 ; j++ ) rnd();
cftyp y = ((((4.03 * x + 0.35) * x - 5.83) * x - 0.63) * x + 1.81) * x +
0.55; // horner
ok += y >= rnd.getF();
for( utyp j=0 ; j<SKIP2 ; j++ ) rnd();
}
// printf("%.14lf\n" , ok / (double)SIZE * 2 - 0.82 );
ftyp sum = 0;
sum += f2(ok - SIZE*0.41) / (SIZE*0.41);
ok = SIZE - ok;
sum += f2(ok - SIZE*0.59) / (SIZE*0.59);
return ddCh2( sum , 1 );
}


// pięć 7-wymiarowych sfer
static ftyp test3( RndTest &rnd , cultyp SIZE ) {
cftyp sf[5][8] = {
{0.6909275387413800,0.2643721119500700,0.35555111472
49490,0.5734918751753870,0.2892100607510660,0.320642
9267302160,0.6386669614352290,0.1846687629586090*0.1
846687629586090},
{0.7423938774969430,0.2995040355715900,0.41840102863
49800,0.6758922111708670,0.3378314392175530,0.636259
3002617360,0.3229378267191350,0.1131459328345950*0.1
131459328345950},
{0.3506312669720500,0.6557646378409120,0.35579511499
96370,0.7002436001319440,0.3976298568770290,0.602149
9079186470,0.2164524640422310,0.1645823939819820*0.1
645823939819820},
{0.4523639138322320,0.3345844920724630,0.35285726310
68530,0.7479594776872550,0.6326152258086950,0.231005
7792346920,0.3267629486043010,0.1733424318023030*0.1
733424318023030},
{0.3796008176170290,0.4209442320279780,0.68776309522
80020,0.3223769583739340,0.4397911571431910,0.792334
8305746910,0.2916006030049170,0.1399564120918510*0.1
399564120918510}
};
ultyp freq[6] = {0,0,0,0,0,0};
for( ultyp i=0 ; i<SIZE ; i++ ) {
utyp k=0;
while( k<5 && sf[k][7] < f2(rnd.getF()-sf[k][0]) + f2(rnd.getF()-sf[k][1]) +
f2(rnd.getF()-sf[k][2]) + f2(rnd.getF()-sf[k][3]) + f2(rnd.getF()-sf[k][4]) +
f2(rnd.getF()-sf[k][5]) + f2(rnd.getF()-sf[k][6]))
k++;
freq[k] ++ ;
}
ftyp sum = 0;
ftyp sump = 0;
ftyp p;
for( utyp i=0 ; i<5 ; i++ ) {
cftyp r = sqrt(sf[i][7]);
p = 16 * pow(M_PI,3) / 105 * pow(r,7);
sum += f2( freq[i] - p * SIZE ) / (p * SIZE);
sump += p;
}
p = (1.0 - sump);
sum += f2( freq[5] - p * SIZE ) / (p * SIZE);
return ddCh2( sum , 5 );
}

//przedziały
static ftyp test4( RndTest &rnd , cultyp SIZE , cutyp SKIP=0) {
cftyp rngs[] = {0,0.1249265665,0.1729261032,0.2464769642,0.34760909
33,0.3792258092,0.4183281037,0.432585527,0.561137729
9,0.6213443817,0.6453121416,0.680396048,0.7939877372
,0.8844306844,0.9670833906,1.00000000001};
cutyp BUCKETS = sizeof(rngs) / sizeof(rngs[0]) - 1;
ultyp freq[BUCKETS];
memset( freq , 0 , sizeof(freq) );
for( ultyp i=0 ; i<SIZE ; i++ ) {
cftyp v = rnd.getF();
for( utyp j=0 ; j<SKIP ; j++ ) rnd();
utyp j=1;
while( v > rngs[j] )
j++;
freq[j-1] ++ ;
}
ftyp sum = 0;
for( utyp i=0 ; i<BUCKETS ; i++ )
sum += f2( freq[i] - (rngs[i+1] - rngs[i]) * SIZE ) / ((rngs[i+1] - rngs[i])
* SIZE);
return ddCh2( sum , BUCKETS-1 );
}


//cztery bity
static ftyp test5( RndTest &rnd , cultyp SIZE , utyp B1=0, utyp B2=4, utyp B3=16,
utyp B4=24, cutyp SKIP=0 ) {
cutyp BUCKETS = 16;
ultyp freq[BUCKETS];
memset( freq , 0 , sizeof(freq) );
cutyp bits1[4] = {B1,B2,B3,B4};
cutyp bits2[4] = {1u<<B1,1u<<B2,1u<<B3,1u<<B4};
for( ultyp i=0 ; i<SIZE ; i++ ) {
utyp v=0;
for( utyp j=0 ; j<4 ; j++ ) {
cutyp r = rnd();
v = (v << 1) | (( r & bits2[j] ) >> bits1[j]);
for( utyp k=0 ; k<SKIP ; k++ )
rnd();
}
freq[v]++;
}
ftyp sum = 0;
for( utyp i=0 ; i<BUCKETS ; i++ ) {
// printf("%llu,",freq[i]);
sum += f2( freq[i] - (ftyp)SIZE/BUCKETS ) / ((ftyp)SIZE/BUCKETS);
}
return ddCh2( sum , BUCKETS-1 );

}


//BIRTHDAY PROBLEM
static ftyp test6( RndTest &rnd , cultyp SIZE , cutyp BUCKETS, cutyp SKIP=0 ) {
ultyp *const freq = new ultyp[BUCKETS];
memset( freq , 0 , sizeof(ultyp) * BUCKETS );
ftyp *const p = new ftyp[BUCKETS]; // rozkład dla paradoksu dnia urodziń
char *const days = new char[BUCKETS];
for( utyp i=0 ; i<BUCKETS ; i++ ) {
p[i] = (ftyp)(i+1) / BUCKETS;
for( utyp j=0 ; j<i ; j++ ) {
p[i] *= 1.0 - ( (ftyp)(j+1) / BUCKETS );
}
// printf("%.12lf\n",(double)p[i]);
}

ftyp sum = 0;
for( utyp i=0 ; i<=BUCKETS ; i++ )
sum += p[i];
if( fabs( sum - 1.0 ) > 0.000000001 )
abort();


for( ultyp i=0 ; i<SIZE ; i++ ) {
memset( days , 0 , sizeof(char)*BUCKETS );
utyp v=0;
while(true) {
cutyp r = rnd() % BUCKETS;
if( days[r] )
break;
days[r] = 1;
v++;
}
if( v-1 >= BUCKETS ) abort();
freq[v-1] ++;
}
// for( utyp i=0 ; i<BUCKETS ; i++ )
// printf("%llu," , freq[i] );

sum = 0;
for( utyp i=0 ; i<BUCKETS ; i++ )
sum += f2( freq[i] - SIZE * p[i] ) / (SIZE*p[i]);
// printf("\nsum %lf\n" , (double)sum );

delete[] p;
delete[] days;
delete[] freq;
return ddCh2( sum , BUCKETS-1 );
}



ftyp rl( cftyp chi ) {
return (chi < 1-chi ? chi : 1-chi) * 100;
}

int main( int argc , char *argv[] ) {
RndTest rnd( time(NULL) );

cultyp m=1000;

printf("test1 %9.6lf%%\n" , rl(test1( rnd , m*1000000 , 16 ) )); fflush(stdout);
printf("test1 %9.6lf%%\n" , rl(test1( rnd , m*1000000 , 10 ) )); fflush(stdout);
printf("test1 %9.6lf%%\n" , rl(test1( rnd , m*1000000 , 16 , 1 ) ));
fflush(stdout);
printf("test1 %9.6lf%%\n" , rl(test1( rnd , m*1000000 , 16 , 2 ) ));
fflush(stdout);
printf("test1 %9.6lf%%\n" , rl(test1( rnd , m*1000000 , 16 , 3 ) ));
fflush(stdout);
printf("test1 %9.6lf%%\n" , rl(test1( rnd , m*1000000 , 20 , 4 ) ));
fflush(stdout);
printf("test1 %9.6lf%%\n" , rl(test1( rnd , m*1000000 , 32 , 5 ) ));
fflush(stdout);
printf("test1 %9.6lf%%\n" , rl(test1( rnd , m*300000 , 16 , 20 ) ));
fflush(stdout);
printf("test1 %9.6lf%%\n" , rl(test1( rnd , m*100000 , 16 , 80 ) ));
fflush(stdout);

printf("test2 %9.6lf%%\n" , rl(test2( rnd , m*100000 ) ));fflush(stdout);
printf("test2 %9.6lf%%\n" , rl(test2( rnd , m*100000 ) ));fflush(stdout);

printf("test3 %9.6lf%%\n" , rl(test3( rnd , m*100000 ) ));fflush(stdout);
printf("test3 %9.6lf%%\n" , rl(test3( rnd , m*100000 ) ));fflush(stdout);

printf("test4 %9.6lf%%\n" , rl(test4( rnd , m*100000 ) ));fflush(stdout);
printf("test4 %9.6lf%%\n" , rl(test4( rnd , m*100000 , 1 ) ));fflush(stdout);
printf("test4 %9.6lf%%\n" , rl(test4( rnd , m*100000 , 2 ) ));fflush(stdout);
printf("test4 %9.6lf%%\n" , rl(test4( rnd , m*100000 , 3 ) ));fflush(stdout);

printf("test5 %9.6lf%%\n" , rl(test5( rnd , m*100000 ) ));fflush(stdout);
printf("test5 %9.6lf%%\n" , rl(test5( rnd , m*100000 , 0 , 1 , 2 , 3 )
));fflush(stdout);
printf("test5 %9.6lf%%\n" , rl(test5( rnd , m*100000 , 0 , 1 , 2 , 3 , 20 )
));fflush(stdout);
printf("test5 %9.6lf%%\n" , rl(test5( rnd , m*100000 , 0 , 1 , 2 , 3 , 80 )
));fflush(stdout);
printf("test5 %9.6lf%%\n" , rl(test5( rnd , m*100000 , 2 , 4 , 6 , 8 )
));fflush(stdout);
printf("test5 %9.6lf%%\n" , rl(test5( rnd , m*100000 , 2 , 4 , 6 , 8 , 20 )
));fflush(stdout);
printf("test5 %9.6lf%%\n" , rl(test5( rnd , m*100000 , 2 , 4 , 6 , 8 , 80 )
));fflush(stdout);

printf("test6 %9.6lf%%\n" , rl(test6( rnd , m*10000000 , 6 ) ));
fflush(stdout);
printf("test6 %9.6lf%%\n" , rl(test6( rnd , m*10000000 , 7 ) ));
fflush(stdout);
printf("test6 %9.6lf%%\n" , rl(test6( rnd , m*10000000 , 10 ) ));
fflush(stdout);
printf("test6 %9.6lf%%\n" , rl(test6( rnd , m*10000000 , 11 ) ));
fflush(stdout);
printf("test6 %9.6lf%%\n" , rl(test6( rnd , m*10000000 , 14 ) ));
fflush(stdout);
printf("test6 %9.6lf%%\n" , rl(test6( rnd , m*10000000 , 15 ) ));
fflush(stdout);
printf("test6 %9.6lf%%\n" , rl(test6( rnd , m*10000000 , 16 ) ));
fflush(stdout);
printf("test6 %9.6lf%%\n" , rl(test6( rnd , m*10000000 , 17 ) ));
fflush(stdout);
printf("test6 %9.6lf%%\n" , rl(test6( rnd , m*10000000 , 31 ) ));
fflush(stdout);
printf("test6 %9.6lf%%\n" , rl(test6( rnd , m*10000000 , 100 ) ));
fflush(stdout);

return 0;
}