xref: /illumos-gate/usr/src/uts/common/inet/ipf/drand48.c (revision ab073b32)
1*ab073b32Sdr /*
2*ab073b32Sdr  * CDDL HEADER START
3*ab073b32Sdr  *
4*ab073b32Sdr  * The contents of this file are subject to the terms of the
5*ab073b32Sdr  * Common Development and Distribution License (the "License").
6*ab073b32Sdr  * You may not use this file except in compliance with the License.
7*ab073b32Sdr  *
8*ab073b32Sdr  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9*ab073b32Sdr  * or http://www.opensolaris.org/os/licensing.
10*ab073b32Sdr  * See the License for the specific language governing permissions
11*ab073b32Sdr  * and limitations under the License.
12*ab073b32Sdr  *
13*ab073b32Sdr  * When distributing Covered Code, include this CDDL HEADER in each
14*ab073b32Sdr  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15*ab073b32Sdr  * If applicable, add the following below this CDDL HEADER, with the
16*ab073b32Sdr  * fields enclosed by brackets "[]" replaced with your own identifying
17*ab073b32Sdr  * information: Portions Copyright [yyyy] [name of copyright owner]
18*ab073b32Sdr  *
19*ab073b32Sdr  * CDDL HEADER END
20*ab073b32Sdr  */
21*ab073b32Sdr 
22*ab073b32Sdr /*
23*ab073b32Sdr  * Copyright 2008 Sun Microsystems, Inc.  All rights reserved.
24*ab073b32Sdr  * Use is subject to license terms.
25*ab073b32Sdr  */
26*ab073b32Sdr 
27*ab073b32Sdr /*	Copyright (c) 1988 AT&T	*/
28*ab073b32Sdr /*	  All Rights Reserved  	*/
29*ab073b32Sdr 
30*ab073b32Sdr #pragma ident	"%Z%%M%	%I%	%E% SMI"
31*ab073b32Sdr 
32*ab073b32Sdr /*
33*ab073b32Sdr  *	drand48, etc. pseudo-random number generator
34*ab073b32Sdr  *	This implementation assumes unsigned short integers of at least
35*ab073b32Sdr  *	16 bits, long integers of at least 32 bits, and ignores
36*ab073b32Sdr  *	overflows on adding or multiplying two unsigned integers.
37*ab073b32Sdr  *	Two's-complement representation is assumed in a few places.
38*ab073b32Sdr  *	Some extra masking is done if unsigneds are exactly 16 bits
39*ab073b32Sdr  *	or longs are exactly 32 bits, but so what?
40*ab073b32Sdr  *	An assembly-language implementation would run significantly faster.
41*ab073b32Sdr  */
42*ab073b32Sdr /*
43*ab073b32Sdr  *	New assumptions (supercede those stated above) for 64-bit work.
44*ab073b32Sdr  *	Longs are now 64 bits, and we are bound by standards to return
45*ab073b32Sdr  *	type long, hovever all internal calculations where long was
46*ab073b32Sdr  *	previously used (32 bit precision) are now using the int32_t
47*ab073b32Sdr  *	type (32 bit precision in both ILP32 and LP64 worlds).
48*ab073b32Sdr  */
49*ab073b32Sdr 
50*ab073b32Sdr #include <sys/mutex.h>
51*ab073b32Sdr 
52*ab073b32Sdr static kmutex_t seed_lock;
53*ab073b32Sdr static int	init48done = 0;
54*ab073b32Sdr 
55*ab073b32Sdr #define	EXPORT0(TYPE, fn, fnu)	TYPE fn() { \
56*ab073b32Sdr 	TYPE res; \
57*ab073b32Sdr 	mutex_enter(&seed_lock); \
58*ab073b32Sdr 	res = fnu(); \
59*ab073b32Sdr 	mutex_exit(&seed_lock); \
60*ab073b32Sdr 	return (res); }
61*ab073b32Sdr #define	EXPORT1(TYPE, fn, fnu)	TYPE fn(unsigned short xsubi[3]) { \
62*ab073b32Sdr 	TYPE res; \
63*ab073b32Sdr 	mutex_enter(&seed_lock); \
64*ab073b32Sdr 	res = fnu(xsubi); \
65*ab073b32Sdr 	mutex_exit(&seed_lock); \
66*ab073b32Sdr 	return (res); }
67*ab073b32Sdr 
68*ab073b32Sdr #define	N	16
69*ab073b32Sdr #define	MASK	((unsigned)(1 << (N - 1)) + (1 << (N - 1)) - 1)
70*ab073b32Sdr #define	LOW(x)	((unsigned)(x) & MASK)
71*ab073b32Sdr #define	HIGH(x)	LOW((x) >> N)
72*ab073b32Sdr #define	MUL(x, y, z)	{ int32_t l = (int32_t)(x) * (int32_t)(y); \
73*ab073b32Sdr 		(z)[0] = LOW(l); (z)[1] = HIGH(l); }
74*ab073b32Sdr #define	CARRY(x, y)	((int32_t)(x) + (int32_t)(y) > MASK)
75*ab073b32Sdr #define	ADDEQU(x, y, z)	(z = CARRY(x, (y)), x = LOW(x + (y)))
76*ab073b32Sdr #define	X0	0x330E
77*ab073b32Sdr #define	X1	0xABCD
78*ab073b32Sdr #define	X2	0x1234
79*ab073b32Sdr #define	A0	0xE66D
80*ab073b32Sdr #define	A1	0xDEEC
81*ab073b32Sdr #define	A2	0x5
82*ab073b32Sdr #define	C	0xB
83*ab073b32Sdr #define	SET3(x, x0, x1, x2)	((x)[0] = (x0), (x)[1] = (x1), (x)[2] = (x2))
84*ab073b32Sdr #define	SETLOW(x, y, n) SET3(x, LOW((y)[n]), LOW((y)[(n)+1]), LOW((y)[(n)+2]))
85*ab073b32Sdr #define	SEED(x0, x1, x2) (SET3(x, x0, x1, x2), SET3(a, A0, A1, A2), c = C)
86*ab073b32Sdr #define	REST(v)	for (i = 0; i < 3; i++) { xsubi[i] = x[i]; x[i] = temp[i]; } \
87*ab073b32Sdr 		return (v)
88*ab073b32Sdr #define	NEST(TYPE, f, F) static TYPE f(unsigned short *xsubi) { \
89*ab073b32Sdr 	int i; TYPE v; unsigned temp[3]; \
90*ab073b32Sdr 	for (i = 0; i < 3; i++) { temp[i] = x[i]; x[i] = LOW(xsubi[i]); }  \
91*ab073b32Sdr 	v = F(); REST(v); }
92*ab073b32Sdr 
93*ab073b32Sdr /* Way ugly solution to problem names, but it works */
94*ab073b32Sdr #define	x	_drand48_x
95*ab073b32Sdr #define	a	_drand48_a
96*ab073b32Sdr #define	c	_drand48_c
97*ab073b32Sdr /* End way ugly */
98*ab073b32Sdr static unsigned x[3] = { X0, X1, X2 }, a[3] = { A0, A1, A2 }, c = C;
99*ab073b32Sdr static unsigned short lastx[3];
100*ab073b32Sdr static void next(void);
101*ab073b32Sdr 
102*ab073b32Sdr static double
103*ab073b32Sdr ipf_r_drand48_u(void)
104*ab073b32Sdr {
105*ab073b32Sdr 	static double two16m = 1.0 / ((int32_t)1 << N);
106*ab073b32Sdr 
107*ab073b32Sdr 	next();
108*ab073b32Sdr 	return (two16m * (two16m * (two16m * x[0] + x[1]) + x[2]));
109*ab073b32Sdr }
110*ab073b32Sdr 
111*ab073b32Sdr NEST(double, ipf_r_erand48_u, ipf_r_drand48_u)
112*ab073b32Sdr 
113*ab073b32Sdr static long
114*ab073b32Sdr ipf_r_lrand48_u(void)
115*ab073b32Sdr {
116*ab073b32Sdr 	next();
117*ab073b32Sdr 	return ((long)((int32_t)x[2] << (N - 1)) + (x[1] >> 1));
118*ab073b32Sdr }
119*ab073b32Sdr 
120*ab073b32Sdr static void
121*ab073b32Sdr init48(void)
122*ab073b32Sdr {
123*ab073b32Sdr 	mutex_init(&seed_lock, 0L, MUTEX_DRIVER, 0L);
124*ab073b32Sdr 	init48done = 1;
125*ab073b32Sdr }
126*ab073b32Sdr 
127*ab073b32Sdr static long
128*ab073b32Sdr ipf_r_mrand48_u(void)
129*ab073b32Sdr {
130*ab073b32Sdr 	next();
131*ab073b32Sdr 	return ((long)((int32_t)x[2] << N) + x[1]);
132*ab073b32Sdr }
133*ab073b32Sdr 
134*ab073b32Sdr static void
135*ab073b32Sdr next(void)
136*ab073b32Sdr {
137*ab073b32Sdr 	unsigned p[2], q[2], r[2], carry0, carry1;
138*ab073b32Sdr 
139*ab073b32Sdr 	MUL(a[0], x[0], p);
140*ab073b32Sdr 	ADDEQU(p[0], c, carry0);
141*ab073b32Sdr 	ADDEQU(p[1], carry0, carry1);
142*ab073b32Sdr 	MUL(a[0], x[1], q);
143*ab073b32Sdr 	ADDEQU(p[1], q[0], carry0);
144*ab073b32Sdr 	MUL(a[1], x[0], r);
145*ab073b32Sdr 	x[2] = LOW(carry0 + carry1 + CARRY(p[1], r[0]) + q[1] + r[1] +
146*ab073b32Sdr 		a[0] * x[2] + a[1] * x[1] + a[2] * x[0]);
147*ab073b32Sdr 	x[1] = LOW(p[1] + r[0]);
148*ab073b32Sdr 	x[0] = LOW(p[0]);
149*ab073b32Sdr }
150*ab073b32Sdr 
151*ab073b32Sdr void
152*ab073b32Sdr ipf_r_srand48(long seedval)
153*ab073b32Sdr {
154*ab073b32Sdr 	int32_t fixseed = (int32_t)seedval;	/* limit to 32 bits */
155*ab073b32Sdr 
156*ab073b32Sdr 	if (init48done == 0)
157*ab073b32Sdr 		init48();
158*ab073b32Sdr 	mutex_enter(&seed_lock);
159*ab073b32Sdr 	SEED(X0, LOW(fixseed), HIGH(fixseed));
160*ab073b32Sdr 	mutex_exit(&seed_lock);
161*ab073b32Sdr }
162*ab073b32Sdr 
163*ab073b32Sdr unsigned short *
164*ab073b32Sdr ipf_r_seed48(unsigned short seed16v[3])
165*ab073b32Sdr {
166*ab073b32Sdr 	if (init48done == 0)
167*ab073b32Sdr 		init48();
168*ab073b32Sdr 	mutex_enter(&seed_lock);
169*ab073b32Sdr 	SETLOW(lastx, x, 0);
170*ab073b32Sdr 	SEED(LOW(seed16v[0]), LOW(seed16v[1]), LOW(seed16v[2]));
171*ab073b32Sdr 	mutex_exit(&seed_lock);
172*ab073b32Sdr 	return (lastx);
173*ab073b32Sdr }
174*ab073b32Sdr 
175*ab073b32Sdr void
176*ab073b32Sdr ipf_r_lcong48(unsigned short param[7])
177*ab073b32Sdr {
178*ab073b32Sdr 	if (init48done == 0)
179*ab073b32Sdr 		init48();
180*ab073b32Sdr 	mutex_enter(&seed_lock);
181*ab073b32Sdr 	SETLOW(x, param, 0);
182*ab073b32Sdr 	SETLOW(a, param, 3);
183*ab073b32Sdr 	c = LOW(param[6]);
184*ab073b32Sdr 	mutex_exit(&seed_lock);
185*ab073b32Sdr }
186*ab073b32Sdr 
187*ab073b32Sdr NEST(long, ipf_r_nrand48_u, ipf_r_lrand48_u)
188*ab073b32Sdr 
189*ab073b32Sdr NEST(long, ipf_r_jrand48_u, ipf_r_mrand48_u)
190*ab073b32Sdr 
191*ab073b32Sdr EXPORT0(double, ipf_r_drand48, ipf_r_drand48_u)
192*ab073b32Sdr EXPORT1(double, ipf_r_erand48, ipf_r_erand48_u)
193*ab073b32Sdr 
194*ab073b32Sdr EXPORT0(long, ipf_r_lrand48, ipf_r_lrand48_u)
195*ab073b32Sdr EXPORT1(long, ipf_r_nrand48, ipf_r_nrand48_u)
196*ab073b32Sdr 
197*ab073b32Sdr EXPORT0(long, ipf_r_mrand48, ipf_r_mrand48_u)
198*ab073b32Sdr EXPORT1(long, ipf_r_jrand48, ipf_r_jrand48_u)
199*ab073b32Sdr 
200*ab073b32Sdr #ifdef DRIVER
201*ab073b32Sdr /*
202*ab073b32Sdr 	This should print the sequences of integers in Tables 2
203*ab073b32Sdr 		and 1 of the TM:
204*ab073b32Sdr 	1623, 3442, 1447, 1829, 1305, ...
205*ab073b32Sdr 	657EB7255101, D72A0C966378, 5A743C062A23, ...
206*ab073b32Sdr  */
207*ab073b32Sdr #include <stdio.h>
208*ab073b32Sdr 
209*ab073b32Sdr main()
210*ab073b32Sdr {
211*ab073b32Sdr 	int i;
212*ab073b32Sdr 
213*ab073b32Sdr 	for (i = 0; i < 80; i++) {
214*ab073b32Sdr 		printf("%4d ", (int)(4096 * ipf_r_drand48()));
215*ab073b32Sdr 		printf("%.4X%.4X%.4X\n", x[2], x[1], x[0]);
216*ab073b32Sdr 	}
217*ab073b32Sdr }
218*ab073b32Sdr #else
219*ab073b32Sdr 
220*ab073b32Sdr #include <sys/random.h>
221*ab073b32Sdr 
222*ab073b32Sdr unsigned
223*ab073b32Sdr ipf_random()
224*ab073b32Sdr {
225*ab073b32Sdr 	static int seeded = 0;
226*ab073b32Sdr 
227*ab073b32Sdr 	if (seeded == 0) {
228*ab073b32Sdr 		long seed;
229*ab073b32Sdr 
230*ab073b32Sdr 		/*
231*ab073b32Sdr 		 * Keep reseeding until some good randomness comes from the
232*ab073b32Sdr 		 * kernel. One of two things will happen: it will succeed or
233*ab073b32Sdr 		 * it will fail (with poor randomness), thus creating NAT
234*ab073b32Sdr 		 * sessions will be "slow" until enough randomness is gained
235*ab073b32Sdr 		 * to not need to get more. It isn't necessary to initialise
236*ab073b32Sdr 		 * seed as it will just pickup whatever random garbage has
237*ab073b32Sdr 		 * been left on the heap and that's good enough until we
238*ab073b32Sdr 		 * get some good garbage.
239*ab073b32Sdr 		 */
240*ab073b32Sdr 		if (random_get_bytes((uint8_t *)&seed, sizeof (seed)) == 0)
241*ab073b32Sdr 			seeded = 1;
242*ab073b32Sdr 		ipf_r_srand48(seed);
243*ab073b32Sdr 	}
244*ab073b32Sdr 
245*ab073b32Sdr 	return (unsigned)ipf_r_lrand48();
246*ab073b32Sdr }
247*ab073b32Sdr #endif
248