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 TWOP29 536870912.0e+00L /* 2^29 */ 16 #define TWOP58 288230376151711744.0e+00L /* 2^(2*29) */ 17 #define TWOP87 154742504910672534362390528.0e+00L /* 2^(3*29) */ 18 19 #define DEFAULT_INITIAL_SEED 100001L 20 21 static long iy = DEFAULT_INITIAL_SEED; 22 23 #if STDC 24 long 25 (qranset)(long seed) 26 #else 27 long 28 (qranset)(seed) 29 long seed; 30 #endif 31 { 32 long old_seed; 33 34 old_seed = iy; 35 if (seed < 0) 36 seed = -seed; 37 if (seed < 0) /* then had seed = LONG_MIN */ 38 seed = LONG_MAX; 39 iy = (seed == 0) ? DEFAULT_INITIAL_SEED : seed; 40 41 /* Guard against bad DEFAULT_INITIAL_SEED */ 42 if (iy < 0) 43 iy = -iy; 44 if (iy < 0) 45 iy = LONG_MAX; 46 47 /* Ensure that (iy * 125) is positive (i.e., does not overflow) */ 48 if (iy > (LONG_MAX / 125)) /* see below for magic constant 125 */ 49 iy = (LONG_MAX / 125); 50 51 return (old_seed); 52 } 53 54 qp_t(qran)55qp_t 56 (qran)(VOID_ARG) 57 { 58 qp_t result; 59 60 iy = iy * 125; 61 iy = iy - (iy / 2796203L) * 2796203L; 62 result = ((qp_t) (iy)) / 2796203.0e+00L; 63 64 iy = iy * 125; 65 iy = iy - (iy / 2796203L) * 2796203L; 66 result += (((qp_t) (iy)) / 2796203.0e+00L) / TWOP29; 67 68 iy = iy * 125; 69 iy = iy - (iy / 2796203L) * 2796203L; 70 result += (((qp_t) (iy)) / 2796203.0e+00L) / TWOP58; 71 72 iy = iy * 125; 73 iy = iy - (iy / 2796203L) * 2796203L; 74 result += (((qp_t) (iy)) / 2796203.0e+00L) / TWOP87; 75 76 return (result); 77 } 78