ThreeB 1.1

- Main Page
- Related Pages
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

mersenne.cpp
Go to the documentation of this file.

    00001 /*
    00002     This file is part of B-cubed.
    00003
    00004     Copyright (C) 2009, 2010, 2011, Edward Rosten and Susan Cox
    00005
    00006     B-cubed is free software; you can redistribute it and/or
    00007     modify it under the terms of the GNU Lesser General Public
    00008     License as published by the Free Software Foundation; either
    00009     version 3.0 of the License, or (at your option) any later version.
    00010
    00011     This library is distributed in the hope that it will be useful,
    00012     but WITHOUT ANY WARRANTY; without even the implied warranty of
    00013     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    00014     Lesser General Public License for more details.
    00015
    00016     You should have received a copy of the GNU General Public License
    00017     along with this program.  If not, see <http://www.gnu.org/licenses/>
    00018 */
    00019
    00020 /**************************   mersenne.cpp   **********************************
    00021 * Author:        Agner Fog
    00022 * Date created:  2001
    00023 * Last modified: 2008-11-16
    00024 * Project:       randomc.h
    00025 * Platform:      Any C++
    00026 * Description:
    00027 * Random Number generator of type 'Mersenne Twister'
    00028 *
    00029 * This random number generator is described in the article by
    00030 * M. Matsumoto & T. Nishimura, in:
    00031 * ACM Transactions on Modeling and Computer Simulation,
    00032 * vol. 8, no. 1, 1998, pp. 3-30.
    00033 * Details on the initialization scheme can be found at
    00034 * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
    00035 *
    00036 * Further documentation:
    00037 * The file ran-instructions.pdf contains further documentation and
    00038 * instructions.
    00039 *
    00040 * Copyright 2001-2008 by Agner Fog.
    00041 * GNU General Public License http://www.gnu.org/licenses/gpl.html
    00042 *******************************************************************************/
    00043
    00044 #include "randomc.h"
    00045
    00046 void CRandomMersenne::Init0(int seed) {
    00047    // Seed generator
    00048    const uint32_t factor = 1812433253UL;
    00049    mt[0]= seed;
    00050    for (mti=1; mti < MERS_N; mti++) {
    00051       mt[mti] = (factor * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
    00052    }
    00053 }
    00054
    00055 void CRandomMersenne::RandomInit(int seed) {
    00056    // Initialize and seed
    00057    Init0(seed);
    00058
    00059    // Randomize some more
    00060    for (int i = 0; i < 37; i++) BRandom();
    00061 }
    00062
    00063
    00064 void CRandomMersenne::RandomInitByArray(int const seeds[], int NumSeeds) {
    00065    // Seed by more than 32 bits
    00066    int i, j, k;
    00067
    00068    // Initialize
    00069    Init0(19650218);
    00070
    00071    if (NumSeeds <= 0) return;
    00072
    00073    // Randomize mt[] using whole seeds[] array
    00074    i = 1;  j = 0;
    00075    k = (MERS_N > NumSeeds ? MERS_N : NumSeeds);
    00076    for (; k; k--) {
    00077       mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL)) + (uint32_t)seeds[j] + j;
    00078       i++; j++;
    00079       if (i >= MERS_N) {mt[0] = mt[MERS_N-1]; i=1;}
    00080       if (j >= NumSeeds) j=0;}
    00081    for (k = MERS_N-1; k; k--) {
    00082       mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL)) - i;
    00083       if (++i >= MERS_N) {mt[0] = mt[MERS_N-1]; i=1;}}
    00084    mt[0] = 0x80000000UL;  // MSB is 1; assuring non-zero initial array
    00085
    00086    // Randomize some more
    00087    mti = 0;
    00088    for (int i = 0; i <= MERS_N; i++) BRandom();
    00089 }
    00090
    00091
    00092 uint32_t CRandomMersenne::BRandom() {
    00093    // Generate 32 random bits
    00094    uint32_t y;
    00095
    00096    if (mti >= MERS_N) {
    00097       // Generate MERS_N words at one time
    00098       const uint32_t LOWER_MASK = (1LU << MERS_R) - 1;       // Lower MERS_R bits
    00099       const uint32_t UPPER_MASK = 0xFFFFFFFF << MERS_R;      // Upper (32 - MERS_R) bits
    00100       static const uint32_t mag01[2] = {0, MERS_A};
    00101
    00102       int kk;
    00103       for (kk=0; kk < MERS_N-MERS_M; kk++) {
    00104          y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
    00105          mt[kk] = mt[kk+MERS_M] ^ (y >> 1) ^ mag01[y & 1];}
    00106
    00107       for (; kk < MERS_N-1; kk++) {
    00108          y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
    00109          mt[kk] = mt[kk+(MERS_M-MERS_N)] ^ (y >> 1) ^ mag01[y & 1];}
    00110
    00111       y = (mt[MERS_N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
    00112       mt[MERS_N-1] = mt[MERS_M-1] ^ (y >> 1) ^ mag01[y & 1];
    00113       mti = 0;
    00114    }
    00115    y = mt[mti++];
    00116
    00117    // Tempering (May be omitted):
    00118    y ^=  y >> MERS_U;
    00119    y ^= (y << MERS_S) & MERS_B;
    00120    y ^= (y << MERS_T) & MERS_C;
    00121    y ^=  y >> MERS_L;
    00122
    00123    return y;
    00124 }
    00125
    00126
    00127 double CRandomMersenne::Random() {
    00128    // Output random float number in the interval 0 <= x < 1
    00129    // Multiply by 2^(-32)
    00130    return (double)BRandom() * (1./(65536.*65536.));
    00131 }
    00132
    00133
    00134 int CRandomMersenne::IRandom(int min, int max) {
    00135    // Output random integer in the interval min <= x <= max
    00136    // Relative error on frequencies < 2^-32
    00137    if (max <= min) {
    00138       if (max == min) return min; else return 0x80000000;
    00139    }
    00140    // Multiply interval with random and truncate
    00141    int r = int((double)(uint32_t)(max - min + 1) * Random() + min);
    00142    if (r > max) r = max;
    00143    return r;
    00144 }
    00145
    00146
    00147 int CRandomMersenne::IRandomX(int min, int max) {
    00148    // Output random integer in the interval min <= x <= max
    00149    // Each output value has exactly the same probability.
    00150    // This is obtained by rejecting certain bit values so that the number
    00151    // of possible bit values is divisible by the interval length
    00152    if (max <= min) {
    00153       if (max == min) return min; else return 0x80000000;
    00154    }
    00155 #ifdef  INT64_SUPPORTED
    00156    // 64 bit integers available. Use multiply and shift method
    00157    uint32_t interval;                    // Length of interval
    00158    uint64_t longran;                     // Random bits * interval
    00159    uint32_t iran;                        // Longran / 2^32
    00160    uint32_t remainder;                   // Longran % 2^32
    00161
    00162    interval = uint32_t(max - min + 1);
    00163    if (interval != LastInterval) {
    00164       // Interval length has changed. Must calculate rejection limit
    00165       // Reject when remainder >= 2^32 / interval * interval
    00166       // RLimit will be 0 if interval is a power of 2. No rejection then
    00167       RLimit = uint32_t(((uint64_t)1 << 32) / interval) * interval - 1;
    00168       LastInterval = interval;
    00169    }
    00170    do { // Rejection loop
    00171       longran  = (uint64_t)BRandom() * interval;
    00172       iran = (uint32_t)(longran >> 32);
    00173       remainder = (uint32_t)longran;
    00174    } while (remainder > RLimit);
    00175    // Convert back to signed and return result
    00176    return (int32_t)iran + min;
    00177
    00178 #else
    00179    // 64 bit integers not available. Use modulo method
    00180    uint32_t interval;                    // Length of interval
    00181    uint32_t bran;                        // Random bits
    00182    uint32_t iran;                        // bran / interval
    00183    uint32_t remainder;                   // bran % interval
    00184
    00185    interval = uint32_t(max - min + 1);
    00186    if (interval != LastInterval) {
    00187       // Interval length has changed. Must calculate rejection limit
    00188       // Reject when iran = 2^32 / interval
    00189       // We can't make 2^32 so we use 2^32-1 and correct afterwards
    00190       RLimit = (uint32_t)0xFFFFFFFF / interval;
    00191       if ((uint32_t)0xFFFFFFFF % interval == interval - 1) RLimit++;
    00192    }
    00193    do { // Rejection loop
    00194       bran = BRandom();
    00195       iran = bran / interval;
    00196       remainder = bran % interval;
    00197    } while (iran >= RLimit);
    00198    // Convert back to signed and return result
    00199    return (int32_t)remainder + min;
    00200
    00201 #endif
    00202 }