Mndiff-dev

Download

mndiff-dev@mplayerhq.hu

April 2010

1 participants
1 discussions

[mndiff]: r180 - in trunk/randi: . randi.c
by michael 14 Apr '10

14 Apr '10

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:; + } +}

1 0