*** FIXES/rand.cc~ 2004-07-15 11:17:11.000000000 +0200 --- FIXES/rand.cc 2004-07-28 11:04:46.000000000 +0200 *************** *** 206,212 **** void fill_randu(int n, double *p) { ! for (int i=0; i < n; i++) p[i] = randu(); } /* --- 206,212 ---- void fill_randu(int n, double *p) { ! for (int i=0; i < n; i++) p[i] = rand53(); } /* *** FIXES/randmtzig.c~ 2004-07-13 23:29:25.000000000 +0200 --- FIXES/randmtzig.c 2004-07-28 11:37:47.000000000 +0200 *************** *** 106,116 **** --- 106,118 ---- #define ZIGGURAT_NOR_R 3.6541528853610088 #define ZIGGURAT_NOR_INV_R 0.27366123732975828 #define TWO_TO_POWER_31 2147483648.0 + #define TWO_TO_POWER_52 4503599627370496.0 #define NOR_SECTION_AREA 0.00492867323399 #define ZIGGURAT_EXP_R 7.69711747013104972 #define ZIGGURAT_EXP_INV_R 0.129918765548341586 #define TWO_TO_POWER_32 4294967296.0 + #define TWO_TO_POWER_53 9007199254740992.0 #define EXP_SECTION_AREA 0.0039496598225815571993 static unsigned long state[MT_N]; /* the array for the state vector */ *************** *** 118,126 **** static int initf = 0; static int initt = 1; static unsigned long *next; ! static unsigned long ki[ZIGGURAT_TABLE_SIZE]; static double wi[ZIGGURAT_TABLE_SIZE], fi[ZIGGURAT_TABLE_SIZE]; ! static unsigned long ke[ZIGGURAT_TABLE_SIZE]; static double we[ZIGGURAT_TABLE_SIZE], fe[ZIGGURAT_TABLE_SIZE]; /* initializes state[MT_N] with a seed */ --- 120,128 ---- static int initf = 0; static int initt = 1; static unsigned long *next; ! static unsigned long long ki[ZIGGURAT_TABLE_SIZE]; static double wi[ZIGGURAT_TABLE_SIZE], fi[ZIGGURAT_TABLE_SIZE]; ! static unsigned long long ke[ZIGGURAT_TABLE_SIZE]; static double we[ZIGGURAT_TABLE_SIZE], fe[ZIGGURAT_TABLE_SIZE]; /* initializes state[MT_N] with a seed */ *************** *** 289,294 **** --- 291,297 ---- generator makes this cycle extremely long. 3) The license on the original code was unclear, thus rewriting the code from the article means we are free of copyright issues. + 4) 53-bit algorithms are used to increase the randomness It should be stated that the authors made my life easier, by the fact that the algorithm developed in the text of the article is for a 256 level *************** *** 296,304 **** One modification to the algorithm developed in the article, is that it is assumed that 0 <= x < Inf, and "unsigned long"s are used, thus resulting in ! terms like 2^32 in the code. As the normal distribution is defined between ! -Inf < x < Inf, we effectively only have 31 bit integers plus a sign. Thus ! in Marsaglia and Tsang, terms like 2^32 become 2^31. It appears that I'm slightly slower than the code in the article, this is partially due to a better generator of random integers than they --- 299,307 ---- One modification to the algorithm developed in the article, is that it is assumed that 0 <= x < Inf, and "unsigned long"s are used, thus resulting in ! terms like 2^53 in the code. As the normal distribution is defined between ! -Inf < x < Inf, we effectively only have 52 bit integers plus a sign. Thus ! in Marsaglia and Tsang, terms like 2^32 become 2^52. It appears that I'm slightly slower than the code in the article, this is partially due to a better generator of random integers than they *************** *** 315,331 **** // Ziggurat tables for the normal distribution x1 = ZIGGURAT_NOR_R; ! wi[255] = x1 / TWO_TO_POWER_31; fi[255] = exp (-0.5 * x1 * x1); /* Index zero is special for tail strip, where Marsaglia and Tsang * defines this as ! * k_0 = 2^31 * r * f(r) / v, w_0 = 0.5^31 * v / f(r), f_0 = 1, * where v is the area of each strip of the ziggurat. */ ! ki[0] = (unsigned long int) (x1 * fi[255] / NOR_SECTION_AREA * ! TWO_TO_POWER_31); ! wi[0] = NOR_SECTION_AREA / fi[255] / TWO_TO_POWER_31; fi[0] = 1.; for (i = 254; i > 0; i--) --- 318,334 ---- // Ziggurat tables for the normal distribution x1 = ZIGGURAT_NOR_R; ! wi[255] = x1 / TWO_TO_POWER_52; fi[255] = exp (-0.5 * x1 * x1); /* Index zero is special for tail strip, where Marsaglia and Tsang * defines this as ! * k_0 = 2^52 * r * f(r) / v, w_0 = 0.5^52 * v / f(r), f_0 = 1, * where v is the area of each strip of the ziggurat. */ ! ki[0] = (unsigned long long) (x1 * fi[255] / NOR_SECTION_AREA * ! TWO_TO_POWER_52); ! wi[0] = NOR_SECTION_AREA / fi[255] / TWO_TO_POWER_52; fi[0] = 1.; for (i = 254; i > 0; i--) *************** *** 334,341 **** * need inverse operator of y = exp(-0.5*x*x) -> x = sqrt(-2*ln(y)) */ x = sqrt(-2. * log(NOR_SECTION_AREA / x1 + fi[i+1])); ! ki[i+1] = (unsigned long int)(x / x1 * TWO_TO_POWER_31); ! wi[i] = x / TWO_TO_POWER_31; fi[i] = exp (-0.5 * x * x); x1 = x; } --- 337,344 ---- * need inverse operator of y = exp(-0.5*x*x) -> x = sqrt(-2*ln(y)) */ x = sqrt(-2. * log(NOR_SECTION_AREA / x1 + fi[i+1])); ! ki[i+1] = (unsigned long long)(x / x1 * TWO_TO_POWER_52); ! wi[i] = x / TWO_TO_POWER_52; fi[i] = exp (-0.5 * x * x); x1 = x; } *************** *** 344,360 **** // Zigurrat tables for the exponential distribution x1 = ZIGGURAT_EXP_R; ! we[255] = x1 / TWO_TO_POWER_32; fe[255] = exp (-x1); /* Index zero is special for tail strip, where Marsaglia and Tsang * defines this as ! * k_0 = 2^32 * r * f(r) / v, w_0 = 0.5^32 * v / f(r), f_0 = 1, * where v is the area of each strip of the ziggurat. */ ! ke[0] = (unsigned long int) (x1 * fe[255] / EXP_SECTION_AREA * ! TWO_TO_POWER_32); ! we[0] = EXP_SECTION_AREA / fe[255] / TWO_TO_POWER_32; fe[0] = 1.; for (i = 254; i > 0; i--) --- 347,363 ---- // Zigurrat tables for the exponential distribution x1 = ZIGGURAT_EXP_R; ! we[255] = x1 / TWO_TO_POWER_53; fe[255] = exp (-x1); /* Index zero is special for tail strip, where Marsaglia and Tsang * defines this as ! * k_0 = 2^53 * r * f(r) / v, w_0 = 0.5^53 * v / f(r), f_0 = 1, * where v is the area of each strip of the ziggurat. */ ! ke[0] = (unsigned long long) (x1 * fe[255] / EXP_SECTION_AREA * ! TWO_TO_POWER_53); ! we[0] = EXP_SECTION_AREA / fe[255] / TWO_TO_POWER_53; fe[0] = 1.; for (i = 254; i > 0; i--) *************** *** 363,370 **** * need inverse operator of y = exp(-x) -> x = -ln(y) */ x = - log(EXP_SECTION_AREA / x1 + fe[i+1]); ! ke[i+1] = (unsigned long int)(x / x1 * TWO_TO_POWER_32); ! we[i] = x / TWO_TO_POWER_32; fe[i] = exp (-x); x1 = x; } --- 366,373 ---- * need inverse operator of y = exp(-x) -> x = -ln(y) */ x = - log(EXP_SECTION_AREA / x1 + fe[i+1]); ! ke[i+1] = (unsigned long long)(x / x1 * TWO_TO_POWER_53); ! we[i] = x / TWO_TO_POWER_53; fe[i] = exp (-x); x1 = x; } *************** *** 377,383 **** * Here is the guts of the algorithm. As Marsaglia and Tsang state the * algorithm in their paper * ! * 1) Calculate a 32-bit random signed integer j and let i be the index * provided by the rightmost 8-bits of j * 2) Set x = j * w_i. If j < k_i return x * 3) If i = 0, then return x from the tail --- 380,386 ---- * Here is the guts of the algorithm. As Marsaglia and Tsang state the * algorithm in their paper * ! * 1) Calculate a 53-bit random integer j and let i be the index * provided by the rightmost 8-bits of j * 2) Set x = j * w_i. If j < k_i return x * 3) If i = 0, then return x from the tail *************** *** 392,404 **** if (initt) create_ziggurat_tables(); while (1) { ! long ri = (signed long) randi (); ! #if SIZEOF_LONG > 4 ! if (ri > 2147483647) ri-=4294967296; ! #endif ! const int idx = ri & 0xFF; ! const double x = ri * wi[idx]; ! if ((ri&LMASK) < ki[idx]) /* LMASK is all but the top bit */ return x; /* 99.3% of the time we return here 1st try */ else if (idx == 0) { --- 395,407 ---- if (initt) create_ziggurat_tables(); while (1) { ! unsigned long r1 = randi(); ! long long ri = (((long long) (r1 & 0x03FFFFFFUL)) << 26) | ! ((long long) (randi() & 0x03FFFFFFUL)); ! int si = ((r1 & 0x80000000UL) ? 1 : 0); ! const int idx = (int)(ri & 0xFFULL); ! const double x = (si ? - ri * wi[idx] : ri * wi[idx]); ! if (ri < ki[idx]) /* All but the top bit */ return x; /* 99.3% of the time we return here 1st try */ else if (idx == 0) { *************** *** 413,425 **** double xx, yy; do { ! xx = - ZIGGURAT_NOR_INV_R * log (randu()); ! yy = - log (randu()); } while ( yy+yy <= xx*xx); ! return (ri < 0 ? -ZIGGURAT_NOR_R-xx : ZIGGURAT_NOR_R+xx); } ! else if ((fi[idx-1] - fi[idx]) * randu() + fi[idx] < exp(-0.5*x*x)) return x; } --- 416,428 ---- double xx, yy; do { ! xx = - ZIGGURAT_NOR_INV_R * log (rand53()); ! yy = - log (rand53()); } while ( yy+yy <= xx*xx); ! return (si ? -ZIGGURAT_NOR_R-xx : ZIGGURAT_NOR_R+xx); } ! else if ((fi[idx-1] - fi[idx]) * rand53() + fi[idx] < exp(-0.5*x*x)) return x; } *************** *** 430,437 **** if (initt) create_ziggurat_tables(); while (1) { ! unsigned long ri = randi (); ! const int idx = ri & 0xFF; const double x = ri * we[idx]; if (ri < ke[idx]) return x; // 98.9% of the time we return here 1st try --- 433,442 ---- if (initt) create_ziggurat_tables(); while (1) { ! unsigned long long ri = ! (((unsigned long long) (randi () & 0x07FFFFFFUL)) << 26) | ! ((unsigned long long) (randi () & 0x03FFFFFFUL)); ! const int idx = (int)(ri & 0xFFULL); const double x = ri * we[idx]; if (ri < ke[idx]) return x; // 98.9% of the time we return here 1st try *************** *** 442,450 **** * For the exponential tail, the method of Marsaglia[5] provides: * x = r - ln(U); */ ! return ZIGGURAT_NOR_R - log(randu()); } ! else if ((fe[idx-1] - fe[idx]) * randu() + fe[idx] < exp(-x)) return x; } } --- 447,455 ---- * For the exponential tail, the method of Marsaglia[5] provides: * x = r - ln(U); */ ! return ZIGGURAT_NOR_R - log(rand53()); } ! else if ((fe[idx-1] - fe[idx]) * rand53() + fe[idx] < exp(-x)) return x; } }