| 1 | #include <stdlib.h> | 
|---|
| 2 |  | 
|---|
| 3 | static randbuf rand48buf; | 
|---|
| 4 | #define A_0  0xE66D | 
|---|
| 5 | #define A_1  0xDEEC | 
|---|
| 6 | #define A_2  0x5 | 
|---|
| 7 | #define C 0xB | 
|---|
| 8 | static randbuf a = { A_0, A_1, A_2 }; | 
|---|
| 9 | static unsigned short c = C; | 
|---|
| 10 |  | 
|---|
| 11 | static void calc_next(randbuf buf) { | 
|---|
| 12 |         randbuf tmp; | 
|---|
| 13 |         long t; | 
|---|
| 14 |         t = buf[0] * a[0] + c; | 
|---|
| 15 |         tmp[0] = t & 0xffff; | 
|---|
| 16 |         tmp[1] = (t >> 16) & 0xffff; | 
|---|
| 17 |         t = buf[1] * a[0] + buf[0] * a[1] + tmp[1]; | 
|---|
| 18 |         tmp[1] = t & 0xffff; | 
|---|
| 19 |         tmp[2] = (t >> 16) & 0xffff; | 
|---|
| 20 |         t = buf[2] * a[0] + buf[1] * a[1] + buf[0] * a[2] + tmp[2]; | 
|---|
| 21 |         tmp[2] = t & 0xffff; | 
|---|
| 22 |         buf[0] = tmp[0]; | 
|---|
| 23 |         buf[1] = tmp[1]; | 
|---|
| 24 |         buf[2] = tmp[2]; | 
|---|
| 25 | } | 
|---|
| 26 |  | 
|---|
| 27 |  | 
|---|
| 28 | long lrand48(void) { | 
|---|
| 29 |         return nrand48(rand48buf); | 
|---|
| 30 | } | 
|---|
| 31 |  | 
|---|
| 32 | long mrand48(void) { | 
|---|
| 33 |         return jrand48(rand48buf); | 
|---|
| 34 | } | 
|---|
| 35 |  | 
|---|
| 36 | double drand48(void) { | 
|---|
| 37 |         return frand48(rand48buf); | 
|---|
| 38 | } | 
|---|
| 39 |  | 
|---|
| 40 | void srand48(long seed) { | 
|---|
| 41 |         rand48buf[1] = (seed >> 16) & 0xffff; | 
|---|
| 42 |         rand48buf[2] = seed & 0xffff; | 
|---|
| 43 |         rand48buf[0] = 0x330e; | 
|---|
| 44 |         a[0] = A_0; | 
|---|
| 45 |         a[1] = A_1; | 
|---|
| 46 |         a[2] = A_2; | 
|---|
| 47 |         c = C; | 
|---|
| 48 | } | 
|---|
| 49 |  | 
|---|
| 50 |  | 
|---|
| 51 | double frand48(randbuf buf) { | 
|---|
| 52 |         double ret; | 
|---|
| 53 |         static double two16m = 1.0 / (1L << 16); | 
|---|
| 54 |         ret = (two16m * (two16m * (two16m * buf[0] + buf[1]) + buf[2])); | 
|---|
| 55 |         calc_next(buf); | 
|---|
| 56 |         return ret; | 
|---|
| 57 | } | 
|---|
| 58 |  | 
|---|
| 59 |  | 
|---|
| 60 |  | 
|---|
| 61 | unsigned short *seed48(randbuf buf) { | 
|---|
| 62 |         static randbuf oldx; | 
|---|
| 63 |         int i; | 
|---|
| 64 |         for (i = 0; i < 3; i++) { | 
|---|
| 65 |                 oldx[i] = rand48buf[i]; | 
|---|
| 66 |                 rand48buf[i] = buf[i]; | 
|---|
| 67 |         } | 
|---|
| 68 |         a[0] = A_0; | 
|---|
| 69 |         a[1] = A_1; | 
|---|
| 70 |         a[2] = A_2; | 
|---|
| 71 |         c = C; | 
|---|
| 72 |         return (unsigned short *)&oldx; | 
|---|
| 73 | } | 
|---|
| 74 |  | 
|---|
| 75 | void lcong48(unsigned short param[7]) { | 
|---|
| 76 |         int i; | 
|---|
| 77 |         for (i = 0; i < 3; i++) { | 
|---|
| 78 |                 rand48buf[i] = param[i]; | 
|---|
| 79 |                 a[i] = param[i + 3]; | 
|---|
| 80 |         } | 
|---|
| 81 |         c = param[6]; | 
|---|
| 82 | } | 
|---|
| 83 |  | 
|---|
| 84 | long jrand48(randbuf buf) { | 
|---|
| 85 |         long ret; | 
|---|
| 86 |         ret = buf[2] << 16 | buf[1]; | 
|---|
| 87 |         calc_next(buf); | 
|---|
| 88 |         return ret; | 
|---|
| 89 | } | 
|---|
| 90 |  | 
|---|
| 91 | long nrand48(randbuf buf) { | 
|---|
| 92 |         return jrand48(buf) & 0x7FFFFFFFL; | 
|---|
| 93 | } | 
|---|