#include #include #include #define INITSEED 19660809UL /************************************************* C code : Chi square test in 1-dimensional Random numbers are generated by Mersenne Twister. **************************************************/ /* Period parameters */ #define P 624 #define Q 397 #define MATRIX_A 0x9908b0dfUL /* constant vector a */ #define UPPER_MASK 0x80000000UL /* most significant w-r bits */ #define LOWER_MASK 0x7fffffffUL /* least significant r bits */ #define nofseq 1000 /* No. of random numbers to generate */ #define k0 10 /* No of equally spaced intervals between 0 and 1 */ static unsigned long mt[P]; /* the array for the state vector */ static int mti = P + 1; /* mti==P+1 means mt[P] is not initialized */ /* initializes mt[P] with a seed */ void init_genrand(unsigned long s) { mt[0] = s & 0xffffffffUL; for (mti = 1; mti < P; mti++) { mt[mti] = (1664525UL * mt[mti-1] + 1UL); mt[mti] &= 0xffffffffUL; } } /* generates a random number on [0, 0xffffffff]-interval */ unsigned long genrand(void) { unsigned long y; static unsigned long mag01[2] = {0x0UL, MATRIX_A}; /* mag01[x] = x * MATRIX_A for x=0, 1 */ if (mti >= P) { /* generate P words at one time */ int kk; if (mti == P + 1) /* if init_genrand() has not been called, */ init_genrand(5489UL); /* a default initial seed is used */ for (kk = 0; kk < P - Q; kk++) { y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK); mt[kk] = mt[kk + Q] ^ (y >> 1) ^ mag01[y & 0x1UL]; } for (; kk < P - 1; kk++) { y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK); mt[kk] = mt[kk + (Q - P)] ^ (y >> 1) ^ mag01[y & 0x1UL]; } y = (mt[P - 1] & UPPER_MASK) | (mt[0] & LOWER_MASK); mt[P - 1] = mt[Q - 1] ^ (y >> 1) ^ mag01[y & 0x1UL]; mti = 0; } y = mt[mti++]; /* Tempering */ y ^= (y >> 11); y ^= (y << 7) & 0x9d2c5680UL; y ^= (y << 15) & 0xefc60000UL; y ^= (y >> 18); return y; } /* generates a random number on [0, 0x7fffffff]-interval */ long genrand_31(void) { return (long)(genrand() >> 1); } void chi2_1dim(int n,int k,float rseq[],float *chi2,int *df,int freq[]){ float interval[k0+1]; int i; int j; interval[0] = 0.0; for (i=1; ir){ rmin=r; } irseq[i]=r; rseq[i]=r; } rmax=rmax+10000; for (i= 0; i