vpi/mt19937int.c

   1
   2 /*
   3  * NOTE: This code as been slightly modified to interface with the
   4  * PLI implementations of $random. The copyright and license
   5  * information are given in the comment block below.
   6  *
   7  * The Modifications include:
   8  *
   9  *    Remove the "main" function, as this is being used for its functions.
  10  *
  11  *    Change the function prototypes to use ANSI/ISO C syntax.
  12  */
  13 #ifdef HAVE_CVS_IDENT
  14 #ident "$Id: mt19937int.c,v 1.5 2004/10/04 01:10:58 steve Exp $"
  15 #endif
  16
  17 /* A C-program for MT19937: Integer version (1998/4/6)            */
  18 /*  genrand() generates one pseudorandom unsigned integer (32bit) */
  19 /* which is uniformly distributed among 0 to 2^32-1  for each     */
  20 /* call. sgenrand(seed) set initial values to the working area    */
  21 /* of 624 words. Before genrand(), sgenrand(seed) must be         */
  22 /* called once. (seed is any 32-bit integer except for 0).        */
  23 /*   Coded by Takuji Nishimura, considering the suggestions by    */
  24 /* Topher Cooper and Marc Rieffel in July-Aug. 1997.              */
  25
  26 /* This library is free software; you can redistribute it and/or   */
  27 /* modify it under the terms of the GNU Library General Public     */
  28 /* License as published by the Free Software Foundation; either    */
  29 /* version 2 of the License, or (at your option) any later         */
  30 /* version.                                                        */
  31 /* This library is distributed in the hope that it will be useful, */
  32 /* but WITHOUT ANY WARRANTY; without even the implied warranty of  */
  33 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.            */
  34 /* See the GNU Library General Public License for more details.    */
  35 /* You should have received a copy of the GNU Library General      */
  36 /* Public License along with this library; if not, write to the    */
  37 /* Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA   */
  38 /* 02111-1307  USA                                                 */
  39
  40 /* Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura.       */
  41 /* When you use this, send an email to: matumoto@math.keio.ac.jp   */
  42 /* with an appropriate reference to your work.                     */
  43
  44 /* REFERENCE                                                       */
  45 /* M. Matsumoto and T. Nishimura,                                  */
  46 /* "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform  */
  47 /* Pseudo-Random Number Generator",                                */
  48 /* ACM Transactions on Modeling and Computer Simulation,           */
  49 /* Vol. 8, No. 1, January 1998, pp 3--30.                          */
  50
  51 #include "sys_priv.h"
  52
  53 /* Period parameters */
  54 #define N 624
  55 #define M 397
  56 #define MATRIX_A 0x9908b0df   /* constant vector a */
  57 #define UPPER_MASK 0x80000000 /* most significant w-r bits */
  58 #define LOWER_MASK 0x7fffffff /* least significant r bits */
  59
  60 /* Tempering parameters */
  61 #define TEMPERING_MASK_B 0x9d2c5680
  62 #define TEMPERING_MASK_C 0xefc60000
  63 #define TEMPERING_SHIFT_U(y)  (y >> 11)
  64 #define TEMPERING_SHIFT_S(y)  (y << 7)
  65 #define TEMPERING_SHIFT_T(y)  (y << 15)
  66 #define TEMPERING_SHIFT_L(y)  (y >> 18)
  67
  68 /* initializing the array with a NONZERO seed */
  69 void
  70 sgenrand(struct context_s *context, unsigned long seed)
  71 {
  72     unsigned long *mt = context->mt;
  73     int mti;
  74
  75     /* setting initial seeds to mt[N] using         */
  76     /* the generator Line 25 of Table 1 in          */
  77     /* [KNUTH 1981, The Art of Computer Programming */
  78     /*    Vol. 2 (2nd Ed.), pp102]                  */
  79     mt[0]= seed & 0xffffffff;
  80     for (mti=1; mti<N; mti++)
  81         mt[mti] = (69069 * mt[mti-1]) & 0xffffffff;
  82
  83     context->mti = mti;
  84 }
  85
  86 unsigned long
  87 genrand(struct context_s *context)
  88 {
  89     unsigned long y;
  90     static unsigned long mag01[2]={0x0, MATRIX_A};
  91     /* mag01[x] = x * MATRIX_A  for x=0,1 */
  92
  93     unsigned long *mt = context->mt;
  94     int mti = context->mti;
  95
  96     if (mti >= N) { /* generate N words at one time */
  97         int kk;
  98
  99         if (mti == N+1)   /* if sgenrand() has not been called, */
 100             sgenrand(context, 4357); /* a default initial seed is used   */
 101
 102         for (kk=0;kk<N-M;kk++) {
 103             y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
 104             mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1];
 105         }
 106         for (;kk<N-1;kk++) {
 107             y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
 108             mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1];
 109         }
 110         y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
 111         mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1];
 112
 113         mti = 0;
 114     }
 115
 116     y = mt[mti++];
 117     y ^= TEMPERING_SHIFT_U(y);
 118     y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B;
 119     y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C;
 120     y ^= TEMPERING_SHIFT_L(y);
 121
 122     context->mti = mti;
 123
 124     return y;
 125 }
 126
 127 /*
 128  * $Log: mt19937int.c,v $
 129  * Revision 1.5  2004/10/04 01:10:58  steve
 130  *  Clean up spurious trailing white space.
 131  *
 132  * Revision 1.4  2003/05/14 04:18:16  steve
 133  *  Use seed to store random number context.
 134  *
 135  * Revision 1.3  2002/08/12 01:35:04  steve
 136  *  conditional ident string using autoconfig.
 137  *
 138  * Revision 1.2  2002/08/11 23:47:04  steve
 139  *  Add missing Log and Ident strings.
 140  *
 141  */