[431c9e0] | 1 | /* drand.c |
---|
| 2 | * |
---|
| 3 | * Pseudorandom number generator |
---|
| 4 | * |
---|
| 5 | * |
---|
| 6 | * |
---|
| 7 | * SYNOPSIS: |
---|
| 8 | * |
---|
| 9 | * double y, drand(); |
---|
| 10 | * |
---|
| 11 | * drand( &y ); |
---|
| 12 | * |
---|
| 13 | * |
---|
| 14 | * |
---|
| 15 | * DESCRIPTION: |
---|
| 16 | * |
---|
| 17 | * Yields a random number 1.0 <= y < 2.0. |
---|
| 18 | * |
---|
| 19 | * The three-generator congruential algorithm by Brian |
---|
| 20 | * Wichmann and David Hill (BYTE magazine, March, 1987, |
---|
| 21 | * pp 127-8) is used. The period, given by them, is |
---|
| 22 | * 6953607871644. |
---|
| 23 | * |
---|
| 24 | * Versions invoked by the different arithmetic compile |
---|
| 25 | * time options DEC, IBMPC, and MIEEE, produce |
---|
| 26 | * approximately the same sequences, differing only in the |
---|
| 27 | * least significant bits of the numbers. The UNK option |
---|
| 28 | * implements the algorithm as recommended in the BYTE |
---|
| 29 | * article. It may be used on all computers. However, |
---|
| 30 | * the low order bits of a double precision number may |
---|
| 31 | * not be adequately random, and may vary due to arithmetic |
---|
| 32 | * implementation details on different computers. |
---|
| 33 | * |
---|
| 34 | * The other compile options generate an additional random |
---|
| 35 | * integer that overwrites the low order bits of the double |
---|
| 36 | * precision number. This reduces the period by a factor of |
---|
| 37 | * two but tends to overcome the problems mentioned. |
---|
| 38 | * |
---|
| 39 | */ |
---|
| 40 | |
---|
| 41 | |
---|
| 42 | /* Three-generator random number algorithm |
---|
| 43 | * of Brian Wichmann and David Hill |
---|
| 44 | * BYTE magazine, March, 1987 pp 127-8 |
---|
| 45 | * |
---|
| 46 | * The period, given by them, is (p-1)(q-1)(r-1)/4 = 6.95e12. |
---|
| 47 | */ |
---|
| 48 | |
---|
| 49 | #include "mconf.h" |
---|
| 50 | #ifdef ANSIPROT |
---|
| 51 | static int ranwh ( void ); |
---|
| 52 | #else |
---|
| 53 | static int ranwh(); |
---|
| 54 | #endif |
---|
| 55 | |
---|
| 56 | static int sx = 1; |
---|
| 57 | static int sy = 10000; |
---|
| 58 | static int sz = 3000; |
---|
| 59 | |
---|
| 60 | static union { |
---|
| 61 | double d; |
---|
| 62 | unsigned short s[4]; |
---|
| 63 | } unkans; |
---|
| 64 | |
---|
| 65 | /* This function implements the three |
---|
| 66 | * congruential generators. |
---|
| 67 | */ |
---|
| 68 | |
---|
| 69 | static int ranwh() |
---|
| 70 | { |
---|
| 71 | int r, s; |
---|
| 72 | |
---|
| 73 | /* sx = sx * 171 mod 30269 */ |
---|
| 74 | r = sx/177; |
---|
| 75 | s = sx - 177 * r; |
---|
| 76 | sx = 171 * s - 2 * r; |
---|
| 77 | if( sx < 0 ) |
---|
| 78 | sx += 30269; |
---|
| 79 | |
---|
| 80 | |
---|
| 81 | /* sy = sy * 172 mod 30307 */ |
---|
| 82 | r = sy/176; |
---|
| 83 | s = sy - 176 * r; |
---|
| 84 | sy = 172 * s - 35 * r; |
---|
| 85 | if( sy < 0 ) |
---|
| 86 | sy += 30307; |
---|
| 87 | |
---|
| 88 | /* sz = 170 * sz mod 30323 */ |
---|
| 89 | r = sz/178; |
---|
| 90 | s = sz - 178 * r; |
---|
| 91 | sz = 170 * s - 63 * r; |
---|
| 92 | if( sz < 0 ) |
---|
| 93 | sz += 30323; |
---|
| 94 | /* The results are in static sx, sy, sz. */ |
---|
| 95 | return 0; |
---|
| 96 | } |
---|
| 97 | |
---|
| 98 | /* drand.c |
---|
| 99 | * |
---|
| 100 | * Random double precision floating point number between 1 and 2. |
---|
| 101 | * |
---|
| 102 | * C callable: |
---|
| 103 | * drand( &x ); |
---|
| 104 | */ |
---|
| 105 | |
---|
| 106 | int drand( a ) |
---|
| 107 | double *a; |
---|
| 108 | { |
---|
| 109 | unsigned short r; |
---|
| 110 | #ifdef DEC |
---|
| 111 | unsigned short s, t; |
---|
| 112 | #endif |
---|
| 113 | |
---|
| 114 | /* This algorithm of Wichmann and Hill computes a floating point |
---|
| 115 | * result: |
---|
| 116 | */ |
---|
| 117 | ranwh(); |
---|
| 118 | unkans.d = sx/30269.0 + sy/30307.0 + sz/30323.0; |
---|
| 119 | r = unkans.d; |
---|
| 120 | unkans.d -= r; |
---|
| 121 | unkans.d += 1.0; |
---|
| 122 | |
---|
| 123 | /* if UNK option, do nothing further. |
---|
| 124 | * Otherwise, make a random 16 bit integer |
---|
| 125 | * to overwrite the least significant word |
---|
| 126 | * of unkans. |
---|
| 127 | */ |
---|
| 128 | #ifdef UNK |
---|
| 129 | /* do nothing */ |
---|
| 130 | #else |
---|
| 131 | ranwh(); |
---|
| 132 | r = sx * sy + sz; |
---|
| 133 | #endif |
---|
| 134 | |
---|
| 135 | #ifdef DEC |
---|
| 136 | /* To make the numbers as similar as possible |
---|
| 137 | * in all arithmetics, the random integer has |
---|
| 138 | * to be inserted 3 bits higher up in a DEC number. |
---|
| 139 | * An alternative would be put it 3 bits lower down |
---|
| 140 | * in all the other number types. |
---|
| 141 | */ |
---|
| 142 | s = unkans.s[2]; |
---|
| 143 | t = s & 07; /* save these bits to put in at the bottom */ |
---|
| 144 | s &= 0177770; |
---|
| 145 | s |= (r >> 13) & 07; |
---|
| 146 | unkans.s[2] = s; |
---|
| 147 | t |= r << 3; |
---|
| 148 | unkans.s[3] = t; |
---|
| 149 | #endif |
---|
| 150 | |
---|
| 151 | #ifdef IBMPC |
---|
| 152 | unkans.s[0] = r; |
---|
| 153 | #endif |
---|
| 154 | |
---|
| 155 | #ifdef MIEEE |
---|
| 156 | unkans.s[3] = r; |
---|
| 157 | #endif |
---|
| 158 | |
---|
| 159 | *a = unkans.d; |
---|
| 160 | return 0; |
---|
| 161 | } |
---|