1*0Sstevel@tonic-gate /*
2*0Sstevel@tonic-gate * CDDL HEADER START
3*0Sstevel@tonic-gate *
4*0Sstevel@tonic-gate * The contents of this file are subject to the terms of the
5*0Sstevel@tonic-gate * Common Development and Distribution License, Version 1.0 only
6*0Sstevel@tonic-gate * (the "License"). You may not use this file except in compliance
7*0Sstevel@tonic-gate * with the License.
8*0Sstevel@tonic-gate *
9*0Sstevel@tonic-gate * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
10*0Sstevel@tonic-gate * or http://www.opensolaris.org/os/licensing.
11*0Sstevel@tonic-gate * See the License for the specific language governing permissions
12*0Sstevel@tonic-gate * and limitations under the License.
13*0Sstevel@tonic-gate *
14*0Sstevel@tonic-gate * When distributing Covered Code, include this CDDL HEADER in each
15*0Sstevel@tonic-gate * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
16*0Sstevel@tonic-gate * If applicable, add the following below this CDDL HEADER, with the
17*0Sstevel@tonic-gate * fields enclosed by brackets "[]" replaced with your own identifying
18*0Sstevel@tonic-gate * information: Portions Copyright [yyyy] [name of copyright owner]
19*0Sstevel@tonic-gate *
20*0Sstevel@tonic-gate * CDDL HEADER END
21*0Sstevel@tonic-gate */
22*0Sstevel@tonic-gate /*
23*0Sstevel@tonic-gate * Copyright 1996 Sun Microsystems, Inc. All rights reserved.
24*0Sstevel@tonic-gate * Use is subject to license terms.
25*0Sstevel@tonic-gate */
26*0Sstevel@tonic-gate
27*0Sstevel@tonic-gate /* Copyright (c) 1984 AT&T */
28*0Sstevel@tonic-gate /* All Rights Reserved */
29*0Sstevel@tonic-gate
30*0Sstevel@tonic-gate #pragma ident "%Z%%M% %I% %E% SMI"
31*0Sstevel@tonic-gate
32*0Sstevel@tonic-gate /*LINTLIBRARY*/
33*0Sstevel@tonic-gate /*
34*0Sstevel@tonic-gate * drand48, etc. pseudo-random number generator
35*0Sstevel@tonic-gate * This implementation assumes unsigned short integers of at least
36*0Sstevel@tonic-gate * 16 bits, long integers of at least 32 bits, and ignores
37*0Sstevel@tonic-gate * overflows on adding or multiplying two unsigned integers.
38*0Sstevel@tonic-gate * Two's-complement representation is assumed in a few places.
39*0Sstevel@tonic-gate * Some extra masking is done if unsigneds are exactly 16 bits
40*0Sstevel@tonic-gate * or longs are exactly 32 bits, but so what?
41*0Sstevel@tonic-gate * An assembly-language implementation would run significantly faster.
42*0Sstevel@tonic-gate */
43*0Sstevel@tonic-gate #define N 16
44*0Sstevel@tonic-gate #define MASK ((unsigned)(1 << (N - 1)) + (1 << (N - 1)) - 1)
45*0Sstevel@tonic-gate #define LOW(x) ((unsigned)(x) & MASK)
46*0Sstevel@tonic-gate #define HIGH(x) LOW((x) >> N)
47*0Sstevel@tonic-gate #define MUL(x, y, z) { long l = (long)(x) * (long)(y); \
48*0Sstevel@tonic-gate (z)[0] = LOW(l); (z)[1] = HIGH(l); }
49*0Sstevel@tonic-gate #define CARRY(x, y) ((long)(x) + (long)(y) > MASK)
50*0Sstevel@tonic-gate #define ADDEQU(x, y, z) (z = CARRY(x, (y)), x = LOW(x + (y)))
51*0Sstevel@tonic-gate #define X0 0x330E
52*0Sstevel@tonic-gate #define X1 0xABCD
53*0Sstevel@tonic-gate #define X2 0x1234
54*0Sstevel@tonic-gate #define A0 0xE66D
55*0Sstevel@tonic-gate #define A1 0xDEEC
56*0Sstevel@tonic-gate #define A2 0x5
57*0Sstevel@tonic-gate #define C 0xB
58*0Sstevel@tonic-gate #define SET3(x, x0, x1, x2) ((x)[0] = (x0), (x)[1] = (x1), (x)[2] = (x2))
59*0Sstevel@tonic-gate #define SETLOW(x, y, n) SET3(x, LOW((y)[n]), LOW((y)[(n)+1]), LOW((y)[(n)+2]))
60*0Sstevel@tonic-gate #define SEED(x0, x1, x2) (SET3(x, x0, x1, x2), SET3(a, A0, A1, A2), c = C)
61*0Sstevel@tonic-gate #define REST(v) \
62*0Sstevel@tonic-gate for (i = 0; i < 3; i++) { \
63*0Sstevel@tonic-gate xsubi[i] = x[i]; \
64*0Sstevel@tonic-gate x[i] = temp[i]; \
65*0Sstevel@tonic-gate } \
66*0Sstevel@tonic-gate return (v)
67*0Sstevel@tonic-gate
68*0Sstevel@tonic-gate #define NEST(TYPE, f, F) \
69*0Sstevel@tonic-gate TYPE f(xsubi) \
70*0Sstevel@tonic-gate register unsigned short *xsubi; \
71*0Sstevel@tonic-gate { \
72*0Sstevel@tonic-gate register int i; \
73*0Sstevel@tonic-gate register TYPE v; \
74*0Sstevel@tonic-gate unsigned temp[3]; \
75*0Sstevel@tonic-gate \
76*0Sstevel@tonic-gate for (i = 0; i < 3; i++) { \
77*0Sstevel@tonic-gate temp[i] = x[i]; \
78*0Sstevel@tonic-gate x[i] = LOW(xsubi[i]); \
79*0Sstevel@tonic-gate } \
80*0Sstevel@tonic-gate v = F(); \
81*0Sstevel@tonic-gate REST(v); \
82*0Sstevel@tonic-gate }
83*0Sstevel@tonic-gate
84*0Sstevel@tonic-gate static unsigned x[3] = { X0, X1, X2 }, a[3] = { A0, A1, A2 }, c = C;
85*0Sstevel@tonic-gate static unsigned short lastx[3];
86*0Sstevel@tonic-gate static void next();
87*0Sstevel@tonic-gate
88*0Sstevel@tonic-gate double
drand48()89*0Sstevel@tonic-gate drand48()
90*0Sstevel@tonic-gate {
91*0Sstevel@tonic-gate static double two16m = 1.0 / (1L << N);
92*0Sstevel@tonic-gate
93*0Sstevel@tonic-gate next();
94*0Sstevel@tonic-gate return (two16m * (two16m * (two16m * x[0] + x[1]) + x[2]));
95*0Sstevel@tonic-gate }
96*0Sstevel@tonic-gate
NEST(double,erand48,drand48)97*0Sstevel@tonic-gate NEST(double, erand48, drand48)
98*0Sstevel@tonic-gate
99*0Sstevel@tonic-gate long
100*0Sstevel@tonic-gate lrand48()
101*0Sstevel@tonic-gate {
102*0Sstevel@tonic-gate next();
103*0Sstevel@tonic-gate return (((long)x[2] << (N - 1)) + (x[1] >> 1));
104*0Sstevel@tonic-gate }
105*0Sstevel@tonic-gate
106*0Sstevel@tonic-gate long
mrand48()107*0Sstevel@tonic-gate mrand48()
108*0Sstevel@tonic-gate {
109*0Sstevel@tonic-gate next();
110*0Sstevel@tonic-gate return (((long)x[2] << N) + x[1]);
111*0Sstevel@tonic-gate }
112*0Sstevel@tonic-gate
113*0Sstevel@tonic-gate static void
next()114*0Sstevel@tonic-gate next()
115*0Sstevel@tonic-gate {
116*0Sstevel@tonic-gate unsigned p[2], q[2], r[2], carry0, carry1;
117*0Sstevel@tonic-gate
118*0Sstevel@tonic-gate MUL(a[0], x[0], p);
119*0Sstevel@tonic-gate ADDEQU(p[0], c, carry0);
120*0Sstevel@tonic-gate ADDEQU(p[1], carry0, carry1);
121*0Sstevel@tonic-gate MUL(a[0], x[1], q);
122*0Sstevel@tonic-gate ADDEQU(p[1], q[0], carry0);
123*0Sstevel@tonic-gate MUL(a[1], x[0], r);
124*0Sstevel@tonic-gate x[2] = LOW(carry0 + carry1 + CARRY(p[1], r[0]) + q[1] + r[1] +
125*0Sstevel@tonic-gate a[0] * x[2] + a[1] * x[1] + a[2] * x[0]);
126*0Sstevel@tonic-gate x[1] = LOW(p[1] + r[0]);
127*0Sstevel@tonic-gate x[0] = LOW(p[0]);
128*0Sstevel@tonic-gate }
129*0Sstevel@tonic-gate
130*0Sstevel@tonic-gate void
srand48(seedval)131*0Sstevel@tonic-gate srand48(seedval)
132*0Sstevel@tonic-gate long seedval;
133*0Sstevel@tonic-gate {
134*0Sstevel@tonic-gate SEED(X0, LOW(seedval), HIGH(seedval));
135*0Sstevel@tonic-gate }
136*0Sstevel@tonic-gate
137*0Sstevel@tonic-gate unsigned short *
seed48(seed16v)138*0Sstevel@tonic-gate seed48(seed16v)
139*0Sstevel@tonic-gate unsigned short seed16v[3];
140*0Sstevel@tonic-gate {
141*0Sstevel@tonic-gate SETLOW(lastx, x, 0);
142*0Sstevel@tonic-gate SEED(LOW(seed16v[0]), LOW(seed16v[1]), LOW(seed16v[2]));
143*0Sstevel@tonic-gate return (lastx);
144*0Sstevel@tonic-gate }
145*0Sstevel@tonic-gate
146*0Sstevel@tonic-gate void
lcong48(param)147*0Sstevel@tonic-gate lcong48(param)
148*0Sstevel@tonic-gate unsigned short param[7];
149*0Sstevel@tonic-gate {
150*0Sstevel@tonic-gate SETLOW(x, param, 0);
151*0Sstevel@tonic-gate SETLOW(a, param, 3);
152*0Sstevel@tonic-gate c = LOW(param[6]);
153*0Sstevel@tonic-gate }
154*0Sstevel@tonic-gate
NEST(long,nrand48,lrand48)155*0Sstevel@tonic-gate NEST(long, nrand48, lrand48)
156*0Sstevel@tonic-gate
157*0Sstevel@tonic-gate NEST(long, jrand48, mrand48)
158*0Sstevel@tonic-gate
159*0Sstevel@tonic-gate #ifdef DRIVER
160*0Sstevel@tonic-gate /*
161*0Sstevel@tonic-gate * This should print the sequences of integers in Tables 2
162*0Sstevel@tonic-gate * and 1 of the TM:
163*0Sstevel@tonic-gate * 1623, 3442, 1447, 1829, 1305, ...
164*0Sstevel@tonic-gate * 657EB7255101, D72A0C966378, 5A743C062A23, ...
165*0Sstevel@tonic-gate */
166*0Sstevel@tonic-gate #include <stdio.h>
167*0Sstevel@tonic-gate
168*0Sstevel@tonic-gate main()
169*0Sstevel@tonic-gate {
170*0Sstevel@tonic-gate int i;
171*0Sstevel@tonic-gate
172*0Sstevel@tonic-gate for (i = 0; i < 80; i++) {
173*0Sstevel@tonic-gate printf("%4d ", (int)(4096 * drand48()));
174*0Sstevel@tonic-gate printf("%.4X%.4X%.4X\n", x[2], x[1], x[0]);
175*0Sstevel@tonic-gate }
176*0Sstevel@tonic-gate }
177*0Sstevel@tonic-gate #endif
178