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