[MN-dev] [mndiff]: r180 - in trunk/randi: . randi.c
michael
subversion at mplayerhq.hu
Wed Apr 14 20:26:31 CEST 2010
Author: michael
Date: Wed Apr 14 20:26:30 2010
New Revision: 180
Log:
PRNG test code
Added:
trunk/randi/
trunk/randi/randi.c
Added: trunk/randi/randi.c
==============================================================================
--- /dev/null 00:00:00 1970 (empty, because file is newly added)
+++ trunk/randi/randi.c Wed Apr 14 20:26:30 2010 (r180)
@@ -0,0 +1,445 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <inttypes.h>
+#include <math.h>
+
+uint64_t add(uint64_t a, uint64_t b, uint64_t m){
+ uint64_t c= a+b;
+ if(c<a || c >= m)
+ c-=m;
+ return c;
+}
+
+// uint64_t gcd(uint64_t a, uint64_t b){
+// if(b) return gcd(b, a%b);
+// else return a;
+// }
+
+uint64_t mul(uint64_t a, uint64_t b, uint64_t m){
+ uint64_t a0= a&0xFFFFFFFF;
+ uint64_t a1= a>>32;
+ uint64_t b0= b&0xFFFFFFFF;
+ uint64_t b1= b>>32;
+ uint64_t t1= a0*b1 + a1*b0;
+ uint64_t t1a= t1<<32;
+ int i;
+
+ a0 = a0*b0 + t1a;
+ a1 = a1*b1 + (t1>>32) + (a0<t1a);
+
+ a1 %= m;
+ for(i=63; i>=0; i--){
+ t1= a1 + ((a0>>i)&1);
+ a1+= t1;
+ if(m <= a1 || a1<t1)
+ a1 -= m;
+ }
+ return a1;
+}
+
+uint64_t egcd(uint64_t a, uint64_t b, uint64_t e, uint64_t f, uint64_t g, uint64_t h){
+ if(b){
+ uint64_t q= a/b;
+
+ return egcd(b, a - b*q, g, h, e - g*q, f - h*q);
+ }else
+ return e;
+}
+
+uint64_t inverse(uint64_t a, uint64_t m){
+ uint64_t inv= egcd(a,m, 1, 0, 0, 1);
+ if(mul(a,inv,m) != 1)
+ inv+= m;
+ if(mul(a,inv,m) != 1)
+ return 0;
+
+ return inv;
+}
+
+//NR LCG
+static uint64_t get_LCG_1664525_1013904223_32(void){
+ static unsigned int seed=1;
+ seed= seed*1664525+1013904223;
+ return seed;
+}
+
+//REIMAR LCG
+static uint64_t get_LCG_2147001325_715136305_32(void){
+ static unsigned int seed=1;
+ seed= seed*2147001325+715136305;
+ return seed;
+}
+//2147001325*a+715136305
+
+static uint64_t get_LCG_921974126499971445_0_63(void){
+ static uint64_t seed=1;
+ seed= seed*9219741426499971445;
+ return seed;
+}
+
+static uint64_t get_LCG_1073217536_0_61(void){
+ static uint64_t seed=1;
+ seed= mul(seed, (1LL<<30) - (1LL<<19), ((1LL<<61)-1));
+ return seed;
+}
+
+static uint64_t get_MLFG_55_24_65_MSB64(void){
+ static uint64_t seed[128]={0,0,0,0,0,0,0,0,0,0,0,0,0,11013904223, 1664525,-1,2,4,8,3,2,2,2,3,4,5,6,2,3,7,5,3,21,4};
+ static int p=55;
+
+ uint64_t a= seed[(p-24) & 127];
+ uint64_t b= seed[(p-55) & 127];
+ seed[p & 127] = 2*a*b + a + b;
+ return seed[p++ & 127];
+}
+
+static uint64_t get_MLFG_55_24_65_MSB32(void){
+ static uint64_t seed[64]={0,0,0,0,0,0,0,0,0,0,0,0,0,11013904223, 1664525,-1,2,4,8,3,2,2,2,3,4,5,6,2,3,7,5,3,21,4};
+ static int p=55;
+
+ uint64_t a= seed[(p-24) & 63];
+ uint64_t b= seed[(p-55) & 63];
+ seed[p & 63] = 2*a*b + a + b;
+ return seed[p++ & 63]>>32;
+}
+
+static uint64_t get_ALFG_55_24_64(void){
+ static uint64_t seed[64]={0,0,0,0,0,0,0,0,0,0,0,0,0,11013904223, 1664525,-1,2,4,8,3,2,2,2,3,4,5,6,2,3,7,5,3,21,4};
+ static int p=55;
+
+ uint64_t a= seed[(p-24) & 63];
+ uint64_t b= seed[(p-55) & 63];
+ seed[p & 63] = a+b;
+ return seed[p++ & 63];
+}
+
+static uint64_t get_ALFG_55_24_64_MSB32(void){
+ static uint64_t seed[64]={0,0,0,0,0,0,0,0,0,0,0,0,0,11013904223, 1664525,-1,2,4,8,3,2,2,2,3,4,5,6,2,3,7,5,3,21,4};
+ static int p=55;
+
+ uint64_t a= seed[(p-24) & 63];
+ uint64_t b= seed[(p-55) & 63];
+ seed[p & 63] = a+b;
+ return seed[p++ & 63]>>32;
+}
+
+static uint64_t get_LFIB4(void){
+ static unsigned int seed[256]={0,0,0,0,0,0,0,0,0,0,0,0,0,11013904223, 1664525};
+ static int p=256;
+ int i;
+
+ if(p<256)
+ return seed[p++];
+
+ for(i=0 ; i<55 ; i++)
+ seed[i]+= seed[i+(256-55)] + seed[i+(256-119)] + seed[i+(256-179)];
+ for(i=55 ; i<119; i++)
+ seed[i]+= seed[i -55 ] + seed[i+(256-119)] + seed[i+(256-179)];
+ for(i=119; i<179; i++)
+ seed[i]+= seed[i -55 ] + seed[i -119 ] + seed[i+(256-179)];
+ for(i=179; i<256; i++)
+ seed[i]+= seed[i -55 ] + seed[i -119 ] + seed[i -179 ];
+
+ p=1;
+ return seed[0];
+}
+
+static uint64_t get_MRGk5_93(void){
+ static unsigned int seed[8]={11013904223, 1664525,1245452,37846321,2376483274,1253124,13487513,1236512};
+ static int p;
+
+ unsigned int a= seed[(p-1)&7];
+ unsigned int b= seed[(p-5)&7];
+
+ seed[p&7]= (107374182ULL*a + 104480ULL*b) % 0x7FFFFFFF;
+ return seed[p++&7];
+}
+
+static uint64_t get_DX47_3(void){
+ static unsigned int seed[64]={11013904223, 1664525,1245452,37846321,2376483274,1253124,13487513,1236512};
+ static int p=8;
+
+ uint64_t a= seed[(p-1)&63] + (uint64_t)seed[(p-24)&63] + seed[(p-47)&63];
+
+ seed[p&63]= (a*((1<<26) + (1<<19))) % 0x7FFFFFFF;
+ return seed[p++&63];
+}
+
+static uint64_t get_KISS99(void){
+ static uint32_t z=362436069, w=521288629, jsr=123456789, jcong=380116160;
+#define znew (z=36969*(z&65535)+(z>>16))
+#define wnew (w=18000*(w&65535)+(w>>16))
+#define MWC ((znew<<16)+wnew )
+#define SHR3 (jsr^=(jsr<<17), jsr^=(jsr>>13), jsr^=(jsr<<5))
+#define CONG (jcong=69069*jcong+1234567)
+ return (MWC^CONG)+SHR3;
+}
+
+static uint64_t get_rand(void){
+ return rand();
+}
+
+typedef struct Generator{
+ uint64_t (*func)(void);
+ const char *name;
+ uint64_t maxout;
+}Generator;
+
+Generator generator[]={
+ {get_LCG_1664525_1013904223_32 , "LCG_1664525_1013904223_32 (NR)", UINT32_MAX},
+ {get_LCG_2147001325_715136305_32, "LCG_2147001325_715136305_32" , UINT32_MAX},
+ {get_LCG_921974126499971445_0_63, "LCG_921974126499971445_0_63" , UINT64_MAX/2},
+ {get_LCG_1073217536_0_61 , "LCG_1073217536_0_61" , (1LL<<61)-2},
+ {get_MLFG_55_24_65_MSB64 , "MLFG_55_24_65_MSB64" , UINT64_MAX},
+ {get_MLFG_55_24_65_MSB32 , "MLFG_55_24_65_MSB32" , UINT32_MAX},
+ {get_ALFG_55_24_64 , "ALFG_55_24_64" , UINT64_MAX},
+ {get_ALFG_55_24_64_MSB32 , "ALFG_55_24_64_MSB32" , UINT32_MAX},
+ {get_LFIB4 , "LFIB4" , UINT32_MAX},
+ {get_DX47_3 , "DX-47-3" , 0x7FFFFFFE},
+ {get_MRGk5_93 , "MRGk5_93" , 0x7FFFFFFE},
+ {get_KISS99 , "KISS99" , UINT32_MAX},
+ {get_rand , "GNU libc rand()" , RAND_MAX},
+};
+#define GEN_COUNT 13
+
+
+#define ELEM uint64_t
+#define ELEM_BITS 64
+
+void test_bitstats(ELEM *matrix, int base_w, int bases, int maxbits){
+ int bits, i, j;
+
+ for(bits=1; bits<4; bits++){
+ int bit_pos[10]={0};
+ long bitshift[10];
+ ELEM *bitp[10];
+ int worst= 0;
+ int bits_per_chunk= (int[]){0, 8, 3, 2, 1, 1 ,1}[bits];
+ int bits_per_chunk_mask= (1<<bits_per_chunk)-1;
+ int border= (int[]){0, 0, base_w/4, base_w*63/128, base_w*63/128, base_w*127/256, base_w*127/256}[bits];
+
+ for(i=1; i<10; i++)
+ bit_pos[i] = bit_pos[i-1] + bits_per_chunk;
+
+ bit_pos[bits]= (base_w - 2*border)*maxbits;
+
+ for(;;){
+ int carry=1;
+ int remaining= 1<<(bits_per_chunk*bits);
+ uint8_t stat[256]={0};
+ for(i=0; i<bits && carry; i++){
+ bit_pos[i]++;
+ if(bit_pos[i]+bits_per_chunk-1 < bit_pos[i+1])
+ carry=0;
+ else
+ bit_pos[i]= i ? bit_pos[i-1]+bits_per_chunk : 0;
+ }
+ if(i==bits && carry)
+ break;
+
+ for(i=0; i<bits; i++){
+ bitp [i]= matrix + bit_pos[i] / maxbits + border;
+ bitshift[i]= bit_pos[i] % maxbits;
+ if(bitshift[i]+bits_per_chunk-1 >= maxbits)
+ goto next_bp;
+ if(i && bit_pos[i] < bit_pos[i-1] + bits_per_chunk - 1)
+ printf("ERROR\n");
+ }
+
+ for(j=0; j<bases; j++){
+ int pattern= (bitp[0][j*base_w] >> bitshift[0]) & bits_per_chunk_mask;
+ for(i=1; i<bits; i++)
+ pattern = (pattern<<bits_per_chunk) + ((bitp[i][j*base_w]>>bitshift[i]) & bits_per_chunk_mask);
+ if(!stat[pattern]){
+ if(!--remaining){
+ if(worst < j)
+ worst= j;
+ break;
+ }
+ stat[pattern]=1;
+ }
+ }
+ if(j>=bases){
+ worst= bases;
+ break;
+ }
+next_bp:;
+ }
+ printf("pattern found after %d(%d) tries, bits:%d\n", worst>>(bits*bits_per_chunk), worst, bits);
+ }
+}
+
+int get_rank(ELEM *matrix, int base_w, int bases, int maxbits){
+ int i, j, k, m, p;
+ int rank=0;
+ ELEM im;
+
+ for(p=i=0; i<base_w; i++){
+ for(im=1LL<<(maxbits-1); im; im>>=1){
+ for(j=p; j<bases; j++){
+ if(matrix[i + j*base_w] & im){
+ for(m=0; m<base_w; m++){
+ ELEM t= matrix[m + p*base_w];
+ matrix[m + p*base_w]= matrix[m + j*base_w];
+ matrix[m + j*base_w]= t;
+ }
+ for(k=p+1; k<bases; k++){
+ if(!(matrix[i + k*base_w] & im))
+ continue;
+ for(m=i; m<base_w; m++){
+ matrix[m + k*base_w] ^= matrix[m + p*base_w];
+ }
+ }
+ p++;
+ break;
+ }
+ }
+ }
+ fprintf(stderr, "%5d/%5d\r", i, base_w);
+ }
+ return p;
+}
+
+static ELEM ror(ELEM x, int b){
+ return (x>>(b-1)) | ((x<<1)&((1LL<<b)-1));
+}
+
+int get_rank_m(ELEM *matrix, int base_w, int bases, uint64_t mod){
+ int i, j, k, m, p;
+ int rank=0;
+
+ for(p=i=0; i<base_w; i++){
+ for(j=p; j<bases; j++){
+ if(matrix[i + j*base_w]){
+ ELEM inv= inverse(matrix[i + j*base_w], mod);
+ if(!inv)
+ continue;
+ for(m=0; m<base_w; m++){
+ ELEM t= matrix[m + p*base_w];
+ matrix[m + p*base_w]= matrix[m + j*base_w];
+ matrix[m + j*base_w]= t;
+ }
+ for(k=p+1; k<bases; k++){
+ ELEM factor;
+ if(!matrix[i + k*base_w])
+ continue;
+ factor= mul((mod-matrix[i + k*base_w]), inv, mod);
+ for(m=i; m<base_w; m++){
+ matrix[m + k*base_w] = add(matrix[m + k*base_w], mul(matrix[m + p*base_w], factor, mod), mod);
+ }
+ }
+ p++;
+ break;
+ }
+ }
+ fprintf(stderr, "%5d/%5d\r", i, base_w);
+ }
+ return p;
+}
+
+int main(){
+ int gr8[8]= {0,1,4,9,15,22,32,34};
+ int gr16[16]= {0,1,4,11,26,32,56,68,76,115,117,134,150,163,168,177};
+#define BASE_W 512
+#define BASES (BASE_W*64*2)
+#define FIXBITS 20
+#define FIXSCALER 1
+#define MASK ((1<<FIXBITS)-1)
+//#define QUAD 1
+
+#define BASE_M_W 512
+#define BASES_M (BASE_M_W*2)
+
+ ELEM *matrix = malloc(BASE_W *BASES *sizeof(ELEM));
+ ELEM *matrix_mod= malloc(BASE_M_W*BASES_M*sizeof(ELEM));
+ int gen, i;
+
+ for(gen=0; gen<GEN_COUNT; gen++){
+ uint64_t hist[BASE_M_W];
+ int i, j, rank, expected, base;
+ uint64_t (*get_r)(void)= generator[gen].func;
+ uint64_t maxout= generator[gen].maxout;
+ int maxbits= log2(maxout)+1;
+ uint64_t mod;
+
+ if(maxbits > 64) maxbits=64; //hello float rounding
+
+ mod= (maxout&1) ? (1LL<<32) : (maxout+1);
+
+ printf("\nTesting %s\n", generator[gen].name);
+
+ for(i=0; i<5000; i++)
+ get_r();
+
+ for(i=0; i<BASE_W; i++)
+ hist[i]= get_r();
+
+ base=0;
+ while(base<BASES){
+ hist[i++%BASE_W]= get_r();
+ if(i>>FIXBITS > 100*(base+1)){
+ printf("FAIL\n");
+ goto next_gen;
+ }
+#if FIXSCALER
+ if(hist[(i-BASE_W/2)%BASE_W] & MASK)
+ continue;
+
+// replace 0 by something to spare us from compensating, as a side effet this might find some weakness in a PRNG
+ hist[(i-BASE_W/2)%BASE_W] ^= i & MASK;
+#else
+ for(j=0; j<FIXBITS; j++)
+ if(hist[(i+j-BASE_W/2-FIXBITS/2)%BASE_W]&1)
+ break;
+ if(j<FIXBITS)
+ continue;
+#endif
+
+#ifdef QUAD
+ for(j=0; j<BASE_W/4; j++)
+ matrix[base*BASE_W + j]= hist[(i+j+BASE_W/2)%BASE_W] & ror(hist[(i+j+BASE_W/2)%BASE_W], maxbits);
+ for(j=BASE_W/4; j<BASE_W*3/4; j++)
+ matrix[base*BASE_W + j]= hist[(i+j)%BASE_W];
+ for(; j<BASE_W; j++)
+ matrix[base*BASE_W + j]= hist[(i+j+BASE_W/2)%BASE_W] & ror(hist[(i+j+BASE_W/2)%BASE_W], maxbits);
+#else
+ for(j=0; j<BASE_W; j++)
+ matrix[base*BASE_W + j]= hist[(i+j)%BASE_W];
+#endif
+#if FIXSCALER
+ hist[(i-BASE_W/2)%BASE_W] &= ~MASK;
+#else
+ERROR
+#endif
+ if(base < BASES_M)
+ for(j=0; j<BASE_M_W; j++)
+ matrix_mod[base*BASE_M_W + j]= matrix[base*BASE_W + j] % mod;
+
+ base++;
+ fprintf(stderr, "%5d/%5d\r", base, BASES);
+ }
+
+ test_bitstats(matrix, BASE_W, BASES, maxbits);
+
+ rank= get_rank(matrix, BASE_W, BASES, maxbits);
+
+ expected= BASE_W * maxbits;
+
+ if(rank<0 || expected < rank) printf("Internal error");
+ else if(expected > rank) printf("FAIL");
+ else printf("PASS");
+ printf(" (%d/%d)\n", rank, expected);
+
+ rank= get_rank_m(matrix_mod, BASE_M_W, BASES_M, mod);
+
+ expected= BASE_M_W;
+
+ if(rank<0 || expected < rank) printf("Internal error");
+ else if(expected > rank) printf("FAIL");
+ else printf("PASS");
+ printf(" (%d/%d)\n", rank, expected);
+
+next_gen:;
+ }
+}
More information about the Mndiff-dev
mailing list