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