QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
mt19937uniformrng.cpp
Go to the documentation of this file.
1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2003 Ferdinando Ametrano
5 Copyright (C) 2010 Kakhkhor Abdijalilov
6
7 This file is part of QuantLib, a free-software/open-source library
8 for financial quantitative analysts and developers - http://quantlib.org/
9
10 QuantLib is free software: you can redistribute it and/or modify it
11 under the terms of the QuantLib license. You should have received a
12 copy of the license along with this program; if not, please email
13 <quantlib-dev@lists.sf.net>. The license is also available online at
14 <http://quantlib.org/license.shtml>.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the license for more details.
19*/
20
21
22// NOTE: The following copyright notice applies to
23// the original C implementation that has been used for this class
24
25/*
26 A C-program for MT19937, with initialization improved 2002/1/26.
27 Coded by Takuji Nishimura and Makoto Matsumoto.
28
29 Before using, initialize the state by using init_genrand(seed)
30 or init_by_array(init_key, key_length).
31
32 Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
33 All rights reserved.
34
35 Redistribution and use in source and binary forms, with or without
36 modification, are permitted provided that the following conditions
37 are met:
38
39 1. Redistributions of source code must retain the above copyright
40 notice, this list of conditions and the following disclaimer.
41
42 2. Redistributions in binary form must reproduce the above copyright
43 notice, this list of conditions and the following disclaimer in the
44 documentation and/or other materials provided with the distribution.
45
46 3. The names of its contributors may not be used to endorse or promote
47 products derived from this software without specific prior written
48 permission.
49
50 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
51 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
52 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
53 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
54 OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
55 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
56 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
57 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
58 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
59 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
60 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
61
62
63 Any feedback is very welcome.
64 http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
65 email: matumoto@math.keio.ac.jp
66*/
67
68
71
72namespace QuantLib {
73
74 // constant vector a
75 const unsigned long MersenneTwisterUniformRng::MATRIX_A = 0x9908b0dfUL;
76 // most significant w-r bits
77 const unsigned long MersenneTwisterUniformRng::UPPER_MASK=0x80000000UL;
78 // least significant r bits
79 const unsigned long MersenneTwisterUniformRng::LOWER_MASK=0x7fffffffUL;
80
81
84 }
85
87 /* initializes mt with a seed */
88 unsigned long s = (seed != 0 ? seed : SeedGenerator::instance().get());
89 mt[0]= s & 0xffffffffUL;
90 for (mti=1; mti<N; mti++) {
91 mt[mti] =
92 (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
93 /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
94 /* In the previous versions, MSBs of the seed affect */
95 /* only MSBs of the array mt[]. */
96 /* 2002/01/09 modified by Makoto Matsumoto */
97 mt[mti] &= 0xffffffffUL;
98 /* for >32 bit machines */
99 }
100 }
101
103 const std::vector<unsigned long>& seeds) {
104 seedInitialization(19650218UL);
105 Size i=1, j=0, k = (N>seeds.size() ? N : seeds.size());
106 for (; k != 0U; k--) {
107 mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL))
108 + seeds[j] + j; /* non linear */
109 mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
110 i++; j++;
111 if (i>=N) { mt[0] = mt[N-1]; i=1; }
112 if (j>=seeds.size()) j=0;
113 }
114 for (k = N - 1; k != 0U; k--) {
115 mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL))
116 - i; /* non linear */
117 mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
118 i++;
119 if (i>=N) { mt[0] = mt[N-1]; i=1; }
120 }
121
122 mt[0] = UPPER_MASK; /*MSB is 1; assuring non-zero initial array*/
123 }
124
126 static const unsigned long mag01[2]={0x0UL, MATRIX_A};
127 /* mag01[x] = x * MATRIX_A for x=0,1 */
128 Size kk;
129 unsigned long y;
130
131 for (kk=0;kk<N-M;kk++) {
132 y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
133 mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL];
134 }
135 for (;kk<N-1;kk++) {
136 y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
137 mt[kk] = mt[(kk+M)-N] ^ (y >> 1) ^ mag01[y & 0x1UL];
138 }
139 y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
140 mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];
141
142 mti = 0;
143 }
144
145}
void seedInitialization(unsigned long seed)
static const unsigned long UPPER_MASK
static const unsigned long LOWER_MASK
static const unsigned long MATRIX_A
MersenneTwisterUniformRng(unsigned long seed=0)
static SeedGenerator & instance()
access to the unique instance
Definition: singleton.hpp:104
std::size_t Size
size of a container
Definition: types.hpp:58
Mersenne Twister uniform random number generator.
Definition: any.hpp:35
Random seed generator.