/ray/src/lib/numerics/rng/mt19937.cc

Go to the documentation of this file.

/*
 * lib/numerics/rng/mt19937.cc
 * 
 * Implementation of fast Mersenne Twister PRNG with period 2^19937-1. 
 * 
 * Converted to class 2004 by Wolfgang Wieser ] wwieser (a) gmx <*> de [ based on 
 * mt19937ar-cok.c by Takuji Nishimura and Makoto Matsumoto; see below. 
 * 
 * This file may be distributed and/or modified under the terms of the 
 * GNU General Public License version 2 as published by the Free Software 
 * Foundation. (See COPYING.GPL for details.)
 * 
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 * 
 * Original copyright and legal info of mt19937ar-cok.c (the C version). 
--------------------------------------------------------------------------------
   A C-program for MT19937, with initialization improved 2002/2/10.
   Coded by Takuji Nishimura and Makoto Matsumoto.
   This is a faster version by taking Shawn Cokus's optimization,
   Matthe Bellew's simplification, Isaku Wada's real version.

   Before using, initialize the state by using init_genrand(seed) 
   or init_by_array(init_key, key_length).

   Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
   All rights reserved.                          

   Redistribution and use in source and binary forms, with or without
   modification, are permitted provided that the following conditions
   are met:

     1. Redistributions of source code must retain the above copyright
        notice, this list of conditions and the following disclaimer.

     2. Redistributions in binary form must reproduce the above copyright
        notice, this list of conditions and the following disclaimer in the
        documentation and/or other materials provided with the distribution.

     3. The names of its contributors may not be used to endorse or promote 
        products derived from this software without specific prior written 
        permission.

   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.


   Any feedback is very welcome.
   http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
   email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
*/

#include "mt19937.h"

namespace NUM   // numerics
{

void RandomNumberGenerator_MersenneTwister19937::seed(uint32 s)
{
    state[0] = s & 0xffffffffU;
    
    for(int j=1; j<N; j++)
    {
        /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
        /* In the previous versions, MSBs of the seed affect   */
        /* only MSBs of the array state[].                     */
        /* 2002/01/09 modified by Makoto Matsumoto             */
        state[j] = (1812433253U * (state[j-1] ^ (state[j-1] >> 30)) + j); 
        // The following should not be needed since we're using uint32 (WW). 
        // (But never mind since it gets optimized away (verified gcc-3.4.2).)
        state[j] &= 0xffffffffU;  /* for >32 bit machines */
    }
    
    left=1;
}


void RandomNumberGenerator_MersenneTwister19937::seed(const uint32 *init_key,
    int key_length)
{
    seed(19650218U);
    
    int i=1;
    int j=0;
    int k = (N>key_length ? N : key_length);
    for(; k; k--)
    {
        state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1664525U))
                   + init_key[j] + j;   /* non linear */
        // The following should not be needed since we're using uint32 (WW). 
        // (But never mind since it gets optimized away (verified gcc-3.4.2).)
        state[i] &= 0xffffffffU; /* for WORDSIZE > 32 machines */
        
        if(++i>=N)
        {  state[0] = state[N-1];  i=1;  }
        if(++j>=key_length)  j=0;
    }
    
    for(k=N-1; k; k--)
    {
        state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) 
                    * 1566083941U)) - i; /* non linear */
        // The following should not be needed since we're using uint32 (WW). 
        // (But never mind since it gets optimized away (verified gcc-3.4.2).)
        state[i] &= 0xffffffffU; /* for WORDSIZE > 32 machines */
        
        if(++i>=N)  { state[0] = state[N-1];  i=1;  }
    }
    
    state[0] = 0x80000000U; /* MSB is 1; assuring non-zero initial array */ 
    left=1;
}


void RandomNumberGenerator_MersenneTwister19937::next_state()
{
    left = N;
    next = state;
    
    uint32 *p=state;
    for(int j=N-M+1; --j; p++) 
        *p = p[M] ^ TWIST(p[0], p[1]);
    
    for(int j=M; --j; p++) 
        *p = p[M-N] ^ TWIST(p[0], p[1]);
    
    *p = p[M-N] ^ TWIST(p[0], state[0]);
}

}  // end of namespace NUM


#if 0
#include <stdio.h>  /* test program */

/*
 * This test program should produce the following: 
 * ===============================================
 * 
 * (Verified 11/2004 by Wolfgang Wieser when converting the C code into 
 * the class RandomNumberGenerator_MersenneTwister19937.) 
 * 
 
1000 outputs of genrand_int32()
1067595299  955945823  477289528 4107218783 4228976476 
3344332714 3355579695  227628506  810200273 2591290167 
2560260675 3242736208  646746669 1479517882 4245472273 
1143372638 3863670494 3221021970 1773610557 1138697238 
1421897700 1269916527 2859934041 1764463362 3874892047 
3965319921   72549643 2383988930 2600218693 3237492380 
2792901476  725331109  605841842  271258942  715137098 
3297999536 1322965544 4229579109 1395091102 3735697720 
2101727825 3730287744 2950434330 1661921839 2895579582 
2370511479 1004092106 2247096681 2111242379 3237345263 
4082424759  219785033 2454039889 3709582971  835606218 
2411949883 2735205030  756421180 2175209704 1873865952 
2762534237 4161807854 3351099340  181129879 3269891896 
 776029799 2218161979 3001745796 1866825872 2133627728 
  34862734 1191934573 3102311354 2916517763 1012402762 
2184831317 4257399449 2899497138 3818095062 3030756734 
[...]

1000 outputs of genrand_real2()
0.76275443 0.99000644 0.98670464 0.10143112 0.27933125 
0.69867227 0.94218740 0.03427201 0.78842173 0.28180608 
0.92179002 0.20785655 0.54534773 0.69644020 0.38107718 
0.23978165 0.65286910 0.07514568 0.22765211 0.94872929 
0.74557914 0.62664415 0.54708246 0.90959343 0.42043116 
0.86334511 0.19189126 0.14718544 0.70259889 0.63426346 
0.77408121 0.04531601 0.04605807 0.88595519 0.69398270 
0.05377184 0.61711170 0.05565708 0.10133577 0.41500776 
0.91810699 0.22320679 0.23353705 0.92871862 0.98897234 
0.19786706 0.80558809 0.06961067 0.55840445 0.90479405 
0.63288060 0.95009721 0.54948447 0.20645042 0.45000959 
0.87050869 0.70806991 0.19406895 0.79286390 0.49332866 
0.78483914 0.75145146 0.12341941 0.42030252 0.16728160 
0.59906494 0.37575460 0.97815160 0.39815952 0.43595080 
0.04952478 0.33917805 0.76509902 0.61034321 0.90654701 
0.92915732 0.85365931 0.18812377 0.65913428 0.28814566 
[...]
*/ 

int main()
{
    CritAssert(sizeof(uint32)==4);  // test app
    uint32 _x=0xffffffff;
    CritAssert(_x>=0);   // test app
    
    uint32 init[4]={0x123, 0x234, 0x345, 0x456}, length=4;
    /* This is an example of initializing by an array. */
    /* You may use seed(seed) with any 32bit integer   */
    /* as a seed for a simpler initialization          */
    NUM::RandomNumberGenerator_MersenneTwister19937 rng(init, length);
    
    printf("1000 outputs of genrand_int32()\n");
    for(int i=0; i<1000; i++)
    {
        printf("%10u ", rng.rand_uint32());
        if(i%5==4) printf("\n");
    }
    
    printf("\n1000 outputs of genrand_real2()\n");
    for(int i=0; i<1000; i++)
    {
        printf("%10.8f ", rng.rand_dbl32o());
        if(i%5==4) printf("\n");
    }
    
    return(0);
}
#endif

---