- Main Page
- Related Pages
- Files

- File List
- Globals

# mt.c

    00001 /*
    00002  * Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
    00003  * All rights reserved.
    00004  *
    00005  * Redistribution and use in source and binary forms, with or without
    00006  * modification, are permitted provided that the following conditions
    00007  * are met:
    00008  *
    00009  *   1. Redistributions of source code must retain the above copyright
    00010  *      notice, this list of conditions and the following disclaimer.
    00011  *
    00012  *   2. Redistributions in binary form must reproduce the above copyright
    00013  *      notice, this list of conditions and the following disclaimer in the
    00014  *      documentation and/or other materials provided with the distribution.
    00015  *
    00016  *   3. The names of its contributors may not be used to endorse or promote
    00017  *      products derived from this software without specific prior written
    00018  *      permission.
    00019  *
    00020  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
    00021  * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
    00022  * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
    00023  * FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE
    00024  * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
    00025  * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
    00026  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
    00027  * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
    00028  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
    00029  * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
    00030  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
    00031  * OF THE POSSIBILITY OF SUCH DAMAGE.
    00032  */
    00033
    00034 #include <stdio.h>
    00035
    00036 /* Period parameters */
    00037 #define MT_N 624
    00038 #define MT_M 397
    00039 #define MT_MATRIX_A 0x9908b0dfUL
    00040 #define MT_UPPER_MASK 0x80000000UL
    00041 #define MT_LOWER_MASK 0x7fffffffUL
    00043 static unsigned long mt[MT_N];
    00044 static int mti = MT_N + 1;
    00050 static void init_genrand(unsigned long s)
    00051 {
    00052     mt[0] = s & 0xffffffffUL;
    00053     for (mti = 1; mti < MT_N; mti++) {
    00054         mt[mti] = (1812433253UL * (mt[mti - 1] ^ (mt[mti - 1] >> 30)) + mti);
    00055         /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
    00056         /* In the previous versions, MSBs of the seed affect   */
    00057         /* only MSBs of the array mt[].                        */
    00058         /* 2002/01/09 modified by Makoto Matsumoto             */
    00059         mt[mti] &= 0xffffffffUL;
    00060         /* for >32 bit machines */
    00061     }
    00062 }
    00063
    00067 static unsigned long genrand_int32(void)
    00068 {
    00069     unsigned long y;
    00070     static unsigned long mag01[2] = { 0x0UL, MT_MATRIX_A };
    00071     /* mag01[x] = x * MT_MATRIX_A  for x=0,1 */
    00072
    00073     if (mti >= MT_N) {          /* generate MT_N words at one time */
    00074         int kk;
    00075
    00076         if (mti == MT_N + 1)    /* if init_genrand() has not been called, */
    00077             init_genrand(5489UL);       /* a default initial seed is used */
    00078
    00079         for (kk = 0; kk < MT_N - MT_M; kk++) {
    00080             y = (mt[kk] & MT_UPPER_MASK) | (mt[kk + 1] & MT_LOWER_MASK);
    00081             mt[kk] = mt[kk + MT_M] ^ (y >> 1) ^ mag01[y & 0x1UL];
    00082         }
    00083         for (; kk < MT_N - 1; kk++) {
    00084             y = (mt[kk] & MT_UPPER_MASK) | (mt[kk + 1] & MT_LOWER_MASK);
    00085             mt[kk] = mt[kk + (MT_M - MT_N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
    00086         }
    00087         y = (mt[MT_N - 1] & MT_UPPER_MASK) | (mt[0] & MT_LOWER_MASK);
    00088         mt[MT_N - 1] = mt[MT_M - 1] ^ (y >> 1) ^ mag01[y & 0x1UL];
    00089
    00090         mti = 0;
    00091     }
    00092
    00093     y = mt[mti++];
    00094
    00095     /* Tempering */
    00096     y ^= (y >> 11);
    00097     y ^= (y << 7) & 0x9d2c5680UL;
    00098     y ^= (y << 15) & 0xefc60000UL;
    00099     y ^= (y >> 18);
    00100
    00101     return y;
    00102 }
    00103
    00107 static double genrand_res53(void)
    00108 {
    00109     unsigned long a = genrand_int32() >> 5, b = genrand_int32() >> 6;
    00110     return (1.0 * a * 67108864.0 + b) * (1.0 / 9007199254740992.0);
    00111 }
    00112
    00113 #undef MT_N
    00114 #undef MT_M
    00115 #undef MT_MATRIX_A
    00116 #undef MT_UPPER_MASK
    00117 #undef MT_LOWER_MASK
    00118
    00119 /*
    00120  * The remainder of this file is
    00121  * Copyright (C) 2010-2011 Nicolas Limare <nicolas.limare@cmla.ens-cachan.fr>
    00122  * and distributed under the same conditions as the original file.
    00123  */
    00124
    00137 #include <time.h>
    00138
    00139 /* get the high-precision timer specified by POSIX.1-2001 if available */
    00140 #if defined(unix) || defined(__unix__) || defined(__unix)
    00141 #include <unistd.h>
    00142 #if defined(_POSIX_VERSION)
    00143 #if (_POSIX_VERSION >= 200112L)
    00144 #include <sys/time.h>
    00145 #define USE_GETTIMEOFDAY
    00146 #endif
    00147 #endif
    00148 #endif
    00149
    00150 /* ensure consistency */
    00151 #include "mt.h"
    00152
    00153 /* string tag inserted into the binary */
    00154 static char mt_tag[] = "using mt " MT_VERSION;
    00163 char *mt_info(void)
    00164 {
    00165     return mt_tag;
    00166 }
    00167
    00171 void mt_init(unsigned long s)
    00172 {
    00173     init_genrand(s);
    00174     return;
    00175 }
    00176
    00180 void mt_init_auto(void)
    00181 {
    00182     unsigned long int seed;
    00183
    00184     /*
    00185      * the basic (weak) initialization
    00186      * uses the current time, in seconds
    00187      */
    00188     seed = (unsigned long int) time(NULL);
    00189
    00190 #if defined(USE_GETTIMEOFDAY)
    00191     /* gettimeofday() provides a millisecond time */
    00192     {
    00193         struct timeval tp;
    00194         (void) gettimeofday(&tp, NULL);
    00195
    00196         seed *= 1000000;
    00197         seed += tp.tv_usec;
    00198     }
    00199 #endif
    00200
    00201     mt_init(seed);
    00202     return;
    00203 }
    00204
    00208 double mt_drand53(void)
    00209 {
    00210     return genrand_res53();
    00211 }