1/*	$NetBSD: mertwist.c,v 1.7 2008/02/17 22:49:11 matt Exp $	*/
2/*-
3 * Copyright (c) 2008 The NetBSD Foundation, Inc.
4 * All rights reserved.
5 *
6 * This code is derived from software contributed to The NetBSD Foundation
7 * by Matt Thomas <matt@3am-software.com>.
8 *
9 * Redistribution and use in source and binary forms, with or without
10 * modification, are permitted provided that the following conditions
11 * are met:
12 * 1. Redistributions of source code must retain the above copyright
13 *    notice, this list of conditions and the following disclaimer.
14 * 2. Redistributions in binary form must reproduce the above copyright
15 *    notice, this list of conditions and the following disclaimer in the
16 *    documentation and/or other materials provided with the distribution.
17 *
18 * THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
19 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
20 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
21 * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
22 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
23 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
24 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
26 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
27 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
28 * POSSIBILITY OF SUCH DAMAGE.
29 */
30
31#if defined(_KERNEL) || defined(_STANDALONE)
32#include <sys/param.h>
33#include <sys/types.h>
34#include <sys/systm.h>
35
36#include <lib/libkern/libkern.h>
37#else
38#include <stdlib.h>
39#include <string.h>
40#include <inttypes.h>
41#include <assert.h>
42#define	KASSERT(x)	assert(x)
43#endif
44
45/*
46 * Mersenne Twister.  Produces identical output compared to mt19937ar.c
47 * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
48 */
49
50#define	MATRIX_A(a)		(((a) & 1) ? 0x9908b0df : 0)
51#define	TEMPERING_MASK_B	0x9d2c5680
52#define	TEMPERING_MASK_C	0xefc60000
53#define	UPPER_MASK		0x80000000
54#define	LOWER_MASK		0x7fffffff
55#define	MIX(u,l)		(((u) & UPPER_MASK) | ((l) & LOWER_MASK))
56
57#define	KNUTH_MULTIPLIER	0x6c078965
58
59#ifndef MTPRNG_RLEN
60#define	MTPRNG_RLEN		624
61#endif
62#define	MTPRNG_POS1		397
63
64static void mtprng_refresh(struct mtprng_state *);
65
66/*
67 * Initialize the generator from a seed
68 */
69void
70mtprng_init32(struct mtprng_state *mt, uint32_t seed)
71{
72	size_t i;
73
74	/*
75	 * Use Knuth's algorithm for expanding this seed over its
76	 * portion of the key space.
77	 */
78	mt->mt_elem[0] = seed;
79	for (i = 1; i < MTPRNG_RLEN; i++) {
80		mt->mt_elem[i] = KNUTH_MULTIPLIER
81		    * (mt->mt_elem[i-1] ^ (mt->mt_elem[i-1] >> 30)) + i;
82	}
83
84	mtprng_refresh(mt);
85}
86
87void
88mtprng_initarray(struct mtprng_state *mt, const uint32_t *key, size_t keylen)
89{
90	uint32_t *mp;
91	size_t i, j, k;
92
93	/*
94	 * Use Knuth's algorithm for expanding this seed over its
95	 * portion of the key space.
96	 */
97	mt->mt_elem[0] = 19650218UL;
98	for (i = 1; i < MTPRNG_RLEN; i++) {
99		mt->mt_elem[i] = KNUTH_MULTIPLIER
100		    * (mt->mt_elem[i-1] ^ (mt->mt_elem[i-1] >> 30)) + i;
101	}
102
103	KASSERT(keylen > 0);
104
105	i = 1;
106	j = 0;
107	k = (keylen < MTPRNG_RLEN ? MTPRNG_RLEN : keylen);
108
109	mp = &mt->mt_elem[1];
110	for (; k-- > 0; mp++) {
111		mp[0] ^= (mp[-1] ^ (mp[-1] >> 30)) * 1664525UL;
112		mp[0] += key[j] + j;
113		if (++i == MTPRNG_RLEN) {
114			KASSERT(mp == mt->mt_elem + MTPRNG_RLEN - 1);
115			mt->mt_elem[0] = mp[0];
116			i = 1;
117			mp = mt->mt_elem;
118		}
119		if (++j == keylen)
120			j = 0;
121	}
122	for (j = MTPRNG_RLEN; --j > 0; mp++) {
123		mp[0] ^= (mp[-1] ^ (mp[-1] >> 30)) * 1566083941UL;
124		mp[0] -= i;
125		if (++i == MTPRNG_RLEN) {
126			KASSERT(mp == mt->mt_elem + MTPRNG_RLEN - 1);
127			mt->mt_elem[0] = mp[0];
128			i = 1;
129			mp = mt->mt_elem;
130		}
131	}
132	mt->mt_elem[0] = 0x80000000;
133	mtprng_refresh(mt);
134}
135
136/*
137 * Generate an array of 624 untempered numbers
138 */
139void
140mtprng_refresh(struct mtprng_state *mt)
141{
142	uint32_t y;
143	size_t i, j;
144	/*
145	 * The following has been refactored to avoid the need for 'mod 624'
146	*/
147	for (i = 0, j = MTPRNG_POS1; j < MTPRNG_RLEN; i++, j++) {
148		y = MIX(mt->mt_elem[i], mt->mt_elem[i+1]);
149		mt->mt_elem[i] = mt->mt_elem[j] ^ (y >> 1) ^ MATRIX_A(y);
150	}
151	for (j = 0; i < MTPRNG_RLEN - 1; i++, j++) {
152		y = MIX(mt->mt_elem[i], mt->mt_elem[i+1]);
153		mt->mt_elem[i] = mt->mt_elem[j] ^ (y >> 1) ^ MATRIX_A(y);
154	}
155	y = MIX(mt->mt_elem[MTPRNG_RLEN - 1], mt->mt_elem[0]);
156	mt->mt_elem[MTPRNG_RLEN - 1] =
157	    mt->mt_elem[MTPRNG_POS1 - 1] ^ (y >> 1) ^ MATRIX_A(y);
158}
159
160/*
161 * Extract a tempered PRN based on the current index.  Then recompute a
162 * new value for the index.  This avoids having to regenerate the array
163 * every 624 iterations. We can do this since recomputing only the next
164 * element and the [(i + 397) % 624] one.
165 */
166uint32_t
167mtprng_rawrandom(struct mtprng_state *mt)
168{
169	uint32_t x, y;
170	const size_t i = mt->mt_idx;
171	size_t j;
172
173	/*
174	 * First generate the random value for the current position.
175	 */
176	x = mt->mt_elem[i];
177	x ^= x >> 11;
178	x ^= (x << 7) & TEMPERING_MASK_B;
179	x ^= (x << 15) & TEMPERING_MASK_C;
180	x ^= x >> 18;
181
182	/*
183	 * Next recalculate the next sequence for the current position.
184	 */
185	y = mt->mt_elem[i];
186	if (__predict_true(i < MTPRNG_RLEN - 1)) {
187		/*
188		 * Avoid doing % since it can be expensive.
189		 * j = (i + MTPRNG_POS1) % MTPRNG_RLEN;
190		 */
191		j = i + MTPRNG_POS1;
192		if (j >= MTPRNG_RLEN)
193			j -= MTPRNG_RLEN;
194		mt->mt_idx++;
195	} else {
196		j = MTPRNG_POS1 - 1;
197		mt->mt_idx = 0;
198	}
199	y = MIX(y, mt->mt_elem[mt->mt_idx]);
200	mt->mt_elem[i] = mt->mt_elem[j] ^ (y >> 1) ^ MATRIX_A(y);
201
202	/*
203	 * Return the value calculated in the first step.
204	 */
205	return x;
206}
207
208/*
209 * This is a non-standard routine which attempts to return a cryptographically
210 * strong random number by collapsing 2 32bit values outputed by the twister
211 * into one 32bit value.
212 */
213uint32_t
214mtprng_random(struct mtprng_state *mt)
215{
216	uint32_t a;
217
218	mt->mt_count = (mt->mt_count + 13) & 31;
219	a = mtprng_rawrandom(mt);
220	a = (a << mt->mt_count) | (a >> (32 - mt->mt_count));
221	return a + mtprng_rawrandom(mt);
222}
223