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 | |
64 | static void mtprng_refresh(struct mtprng_state *); |
65 | |
66 | /* |
67 | * Initialize the generator from a seed |
68 | */ |
69 | void |
70 | mtprng_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 | |
87 | void |
88 | mtprng_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 | */ |
139 | void |
140 | mtprng_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 | */ |
166 | uint32_t |
167 | mtprng_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 | */ |
213 | uint32_t |
214 | mtprng_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 | |