1/* $NetBSD: mertwist.c,v 1.8 2008/04/28 20:24:06 martin 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