| 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 | } | 
|---|