jardin | eb5d44e | 2003-12-23 08:09:43 +0000 | [diff] [blame] | 1 | /*********************************************************************/ |
| 2 | /* */ |
| 3 | /* current processor time in seconds */ |
| 4 | /* difference between two calls is processor time spent by your code */ |
| 5 | /* needs: <sys/types.h>, <sys/times.h> */ |
| 6 | /* depends on compiler and OS */ |
| 7 | /* */ |
| 8 | /*********************************************************************/ |
| 9 | |
| 10 | #include <sys/types.h> |
| 11 | #include <sys/times.h> |
| 12 | |
Paul Jakma | 41b36e9 | 2006-12-08 01:09:50 +0000 | [diff] [blame^] | 13 | /* |
| 14 | * Prototypes. |
| 15 | */ |
| 16 | unsigned long timer(void); |
| 17 | void init_rand(long); |
| 18 | double rand01(void); |
| 19 | double randg01(void); |
| 20 | long nrand(long); |
| 21 | void free_arc(void *); |
| 22 | |
| 23 | unsigned long timer () |
jardin | eb5d44e | 2003-12-23 08:09:43 +0000 | [diff] [blame] | 24 | { struct tms hold; |
| 25 | |
| 26 | times(&hold); |
hasso | 29e50b2 | 2005-09-01 18:18:47 +0000 | [diff] [blame] | 27 | return (unsigned long) ((float) (hold.tms_utime) / 60.0); |
jardin | eb5d44e | 2003-12-23 08:09:43 +0000 | [diff] [blame] | 28 | } |
| 29 | |
| 30 | |
| 31 | /*********************************************************************/ |
| 32 | /* */ |
| 33 | /* Family of random number generators */ |
| 34 | /* */ |
| 35 | /* Initialisation: */ |
| 36 | /* void init_rand ( seed ); */ |
| 37 | /* long seed - any positive number */ |
| 38 | /* if seed<=0 init_rand takes time */ |
| 39 | /* from timer instead of seed */ |
| 40 | /* */ |
| 41 | /* Whole number uniformly distributed on [0,n): */ |
| 42 | /* long nrand (n); */ |
| 43 | /* long n */ |
| 44 | /* */ |
| 45 | /* Real number uniformly distributed on [0,1] */ |
| 46 | /* double rand01(); */ |
| 47 | /* */ |
| 48 | /* Real number with Gauss(0,1) disitribution: */ |
| 49 | /* double randg01(); */ |
| 50 | /* */ |
| 51 | /* Algorithm: */ |
| 52 | /* x(n+1) = (x(n) * 5^13) mod 2^31 */ |
| 53 | /* */ |
| 54 | /*********************************************************************/ |
| 55 | |
| 56 | unsigned long internal_seed; |
| 57 | |
| 58 | void init_rand ( init_seed ) |
| 59 | |
| 60 | long init_seed; |
| 61 | |
| 62 | { internal_seed = ( init_seed > 0 ) |
| 63 | ? (unsigned long) init_seed |
| 64 | : (unsigned long) timer(); |
| 65 | |
| 66 | |
| 67 | /* only odd numbers are acceptable */ |
| 68 | if ( internal_seed % 2 == 0 ) internal_seed --; |
| 69 | } |
| 70 | |
| 71 | /*********************************************************************/ |
| 72 | /* */ |
| 73 | /* Internal function irand may depend on OS and compiler */ |
| 74 | /* */ |
| 75 | /* irand assumption: */ |
| 76 | /* unsigned long i,j; */ |
| 77 | /* if i*j > max(unsigned long) */ |
| 78 | /* 1. No overflow interruption */ |
| 79 | /* 2. i*j = i*j mod max(unsigned long) */ |
| 80 | /* */ |
| 81 | /* This assumption is true for a lot of computers. */ |
| 82 | /* If your computer fails: */ |
| 83 | /* rename: irand <---> xrand */ |
| 84 | /* */ |
| 85 | /*********************************************************************/ |
| 86 | |
| 87 | #define A 1220703125 |
| 88 | #define B 2147483647 |
| 89 | #define BF 2147483647. |
| 90 | |
| 91 | static long irand () |
| 92 | |
| 93 | { internal_seed = ( internal_seed * A ) & B; |
| 94 | return (long) internal_seed ; |
| 95 | } |
| 96 | |
hasso | 29e50b2 | 2005-09-01 18:18:47 +0000 | [diff] [blame] | 97 | #if 0 /* Not used. */ |
jardin | eb5d44e | 2003-12-23 08:09:43 +0000 | [diff] [blame] | 98 | /*********************************************************************/ |
| 99 | /* */ |
| 100 | /* computer independent variant of irand */ |
| 101 | /* */ |
| 102 | /*********************************************************************/ |
| 103 | |
| 104 | |
| 105 | #define T15 32768 |
| 106 | #define T16 65536 |
| 107 | #define A1 37252 |
| 108 | #define A2 29589 |
| 109 | |
| 110 | static long xrand() |
| 111 | |
| 112 | { unsigned long is1, is2; |
| 113 | |
| 114 | is1 = internal_seed / T15; |
| 115 | is2 = internal_seed % T15; |
| 116 | |
| 117 | internal_seed = ( (((is2 * A1) + (is1 * A2))% T16 )* T15 + (is2 * A2) ) & B; |
| 118 | return (long) ( internal_seed ) ; |
| 119 | } |
hasso | 29e50b2 | 2005-09-01 18:18:47 +0000 | [diff] [blame] | 120 | #endif |
jardin | eb5d44e | 2003-12-23 08:09:43 +0000 | [diff] [blame] | 121 | |
| 122 | /*********************************************************************/ |
| 123 | |
| 124 | |
| 125 | double rand01() |
| 126 | |
hasso | 29e50b2 | 2005-09-01 18:18:47 +0000 | [diff] [blame] | 127 | { return (double) (irand() / BF) ; |
jardin | eb5d44e | 2003-12-23 08:09:43 +0000 | [diff] [blame] | 128 | } |
| 129 | |
| 130 | /*********************************************************************/ |
| 131 | |
| 132 | #define NK 12 |
| 133 | |
| 134 | double randg01() |
| 135 | |
| 136 | { int i; |
| 137 | double sum = 0; |
| 138 | |
| 139 | for ( i = 0; i < NK; i++ ) sum += rand01(); |
| 140 | return sum - 6.; |
| 141 | |
| 142 | /* if NK != 12 then you must return (12/NK)*sum - (NK/2) */ |
| 143 | } |
| 144 | |
| 145 | #undef NK |
| 146 | |
| 147 | |
| 148 | /*********************************************************************/ |
| 149 | |
| 150 | long nrand ( n ) |
| 151 | |
| 152 | long n; |
| 153 | |
| 154 | { return (long) ( rand01() * (double) n ); |
| 155 | } |
| 156 | |
| 157 | /*********************************************************************/ |
| 158 | |
| 159 | #undef A |
| 160 | #undef A1 |
| 161 | #undef A2 |
| 162 | #undef B |
| 163 | #undef BF |
| 164 | #undef T15 |
| 165 | #undef T16 |