don@allegra.UUCP (D. Mitchell) (04/19/84)
1,348c /* * Additive Random Number Generator (see: Knuth 3.2.2) * D. P. Mitchell 84/04/16. */ #define LENGTH 97 #define TAP 12 #define MANTISSA 56 static long y[LENGTH], *tap = &y[0], *feedback = &y[LENGTH - TAP]; static double fy[LENGTH], *ftap = &fy[0], *ffeedback = &fy[LENGTH - TAP]; static unsigned long x = 0xc0393e8d; random() { if (--tap < y) { if (x) /* effect of srandom might be delayed */ setup(); tap += LENGTH; } if (--feedback < y) feedback += LENGTH; return (*feedback += *tap) >> 1 & 0x7fffffff; } double frandom() { double sum; if (--ftap < fy) { if (x) /* fy cannot be preset portably! */ setup(); ftap += LENGTH; } if (--ffeedback < fy) ffeedback += LENGTH; sum = *ftap + *ffeedback; if (sum >= 1.0) { sum = *ftap - 0.5; /* must not overflow! */ sum += *ffeedback - 0.5; } *ffeedback = sum; return sum; } srandom(seed) long seed; { x = seed; } static setup() { long mask; double divisor; int i; for (i = 0; i < LENGTH; i++) { x = 177988*x*x + 505360173*x + 907633385; y[i] = x; } mask = 0x7fffffff >> (63 - MANTISSA); divisor = (double)mask + 1.0; for (i = 0; i < LENGTH; i++) { x = 177988*x*x + 505360173*x + 907633385; fy[i] = (double)x / 4294967296.0; x = 177988*x*x + 505360173*x + 907633385; fy[i] = ((double)(x & mask) + fy[i]) / divisor; } x = 0; } .