1 /* -*-C-*- dran.c */ 2 3 #include "elefunt.h" 4 5 /*********************************************************************** 6 Random number generator - based on Algorithm 266 by Pike and Hill 7 (modified by Hansson), Communications of the ACM, Vol. 8, No. 10, 8 October 1965. 9 10 This subprogram is intended for use on computers with fixed point 11 wordlength of at least 29 bits. It is best if the floating point 12 significand has at most 29 bits. 13 ***********************************************************************/ 14 15 #define DEFAULT_INITIAL_SEED 100001L 16 17 static long iy = DEFAULT_INITIAL_SEED; 18 19 #if STDC 20 long 21 (ranset)(long seed) 22 #else 23 long 24 (ranset)(seed) 25 long seed; 26 #endif 27 { 28 long old_seed; 29 30 old_seed = iy; 31 if (seed < 0) 32 seed = -seed; 33 if (seed < 0) /* then had seed = LONG_MIN */ 34 seed = LONG_MAX; 35 iy = (seed == 0) ? DEFAULT_INITIAL_SEED : seed; 36 37 /* Guard against bad DEFAULT_INITIAL_SEED */ 38 if (iy < 0) 39 iy = -iy; 40 if (iy < 0) 41 iy = LONG_MAX; 42 43 /* Ensure that (iy * 125) is positive (i.e., does not overflow) */ 44 if (iy > (LONG_MAX / 125)) /* see below for magic constant 125 */ 45 iy = (LONG_MAX / 125); 46 47 return (old_seed); 48 } 49 50 dp_t(ran)51dp_t 52 (ran)(VOID_ARG) 53 { 54 dp_t result; 55 56 iy = iy * 125; 57 iy = iy - (iy / 2796203L) * 2796203L; 58 result = ((dp_t) (iy)) / 2796203.0e0; 59 60 iy = iy * 125; 61 iy = iy - (iy / 2796203L) * 2796203L; 62 result += (((dp_t) (iy)) / 2796203.0e0) / 536870912.0e+00; 63 64 return (result); 65 } 66