-
Data: 2016-10-09 15:26:41
Temat: kontynuacja generatory: mersen vs ranlux
Od: "M.M." <m...@g...com> szukaj wiadomości tego autora
[ pokaż wszystkie nagłówki ]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;
}
Następne wpisy z tego wątku
- 09.10.16 15:55 bartekltg
- 09.10.16 23:12 M.M.
- 10.10.16 09:31 M.M.
- 11.10.16 17:44 Borneq
- 11.10.16 20:20 M.M.
- 11.10.16 20:53 M.M.
- 11.10.16 21:47 bartekltg
- 11.10.16 22:00 bartekltg
- 12.10.16 04:14 M.M.
Najnowsze wątki z tej grupy
- Na grupie comp.os.linux.advocacy CrudeSausage twierdzi, że Micro$lop używa SI do szyfrowania formatu dok. XML
- Błąd w Sofcie Powodem Wymiany 3 Duńskich Fregat Typu Iver Huitfeldt
- Grok zaczął nadużywać wulgaryzmów i wprost obrażać niektóre znane osoby
- Can you activate BMW 48V 10Ah Li-Ion battery, connecting to CAN-USB laptop interface ?
- We Wrocławiu ruszyła Odra 5, pierwszy w Polsce komputer kwantowy z nadprzewodzącymi kubitami
- Ada-Europe - AEiC 2025 early registration deadline imminent
- John Carmack twierdzi, że gdyby gry były optymalizowane, to wystarczyły by stare kompy
- Ada-Europe Int.Conf. Reliable Software Technologies, AEiC 2025
- Linuks od wer. 6.15 przestanie wspierać procesory 486 i będzie wymagać min. Pentium
- ,,Polski przemysł jest w stanie agonalnym" - podkreślił dobitnie, wskazując na brak zamówień.
- Rewolucja w debugowaniu!!! SI analizuje zrzuty pamięci systemu M$ Windows!!!
- Brednie w wiki - hasło Dehomag
- Perfidne ataki krakerów z KRLD na skrypciarzy JS i Pajton
- Instytut IDEAS może zacząć działać: "Ma to być unikalny w europejskiej skali ośrodek badań nad sztuczną inteligencją."
- Instytut IDEAS może zacząć działać: "Ma to być unikalny w europejskiej skali ośrodek badań nad sztuczną inteligencją."
Najnowsze wątki
- 2025-08-06 Gdynia => Konsultant wdrożeniowy (systemy controlingowe) <=
- 2025-08-06 Białystok => Inżynier oprogramowania .Net <=
- 2025-08-06 "[...] sejmowe wystąpienie posłanki Klaudii Jachiry, która zakończyła je słowami ,,Sława Ukrainie"."
- 2025-08-05 "Chiny przekraczają w wydobyciu 4 mld ton węgla, Indie i USA ponad 1 mld, a Rosja 500 mln ton [...]"
- 2025-08-05 Panuje się 181 159,42 zł./mies. na posła w 2026r.
- 2025-08-05 "Chiny przekraczają w wydobyciu 4 mld ton węgla, Indie i USA ponad 1 mld, a Rosja 500 mln ton [...]"
- 2025-08-05 Czy cos fi przechodzi przez trafo separujące?
- 2025-08-05 kajaki i promile
- 2025-08-05 Re: Tesla jest bezpieczna, wczoraj spaliła się doszczętnie na Ursynowie i nikomu się nic nie stało
- 2025-08-05 Gdynia => Przedstawiciel handlowy / KAM (branża TSL) <=
- 2025-08-05 Re: Atak na lekarza w Oławie. Policja zatrzymała sprawcę na lotnisku Polska Agencja Prasowa 4 sierpnia 2025, 12:16 FACEBOOK X E-MAIL KOPIUJ LINK W szpitalu w Oławie 37-letni pacjent zaatakował lekarza, po tym, jak ten odmówił mu wypisania długoterminowego
- 2025-08-05 B2B i książka przychodów i rozchodów
- 2025-08-04 Re: Atak na lekarza w Oławie. Policja zatrzymała sprawcę na lotnisku Polska Agencja Prasowa 4 sierpnia 2025, 12:16 FACEBOOK X E-MAIL KOPIUJ LINK W szpitalu w Oławie 37-letni pacjent zaatakował lekarza, po tym, jak ten odmówił mu wypisania długoterminowego
- 2025-08-04 Na grupie comp.os.linux.advocacy CrudeSausage twierdzi, że Micro$lop używa SI do szyfrowania formatu dok. XML
- 2025-08-04 Na grupie comp.os.linux.advocacy CrudeSausage twierdzi, że Micro$lop używa SI do szyfrowania formatu dok. XML