sys/lib/libkern/mertwist.c

   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
  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 }