1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
|
/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
/*
Copyright (C) 2000, 2001, 2002, 2003 RiskMap srl
This file is part of QuantLib, a free-software/open-source library
for financial quantitative analysts and developers - http://quantlib.org/
QuantLib is free software: you can redistribute it and/or modify it
under the terms of the QuantLib license. You should have received a
copy of the license along with this program; if not, please email
<quantlib-dev@lists.sf.net>. The license is also available online at
<http://quantlib.org/license.shtml>.
This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the license for more details.
*/
#include <ql/math/randomnumbers/lecuyeruniformrng.hpp>
#include <ql/math/randomnumbers/seedgenerator.hpp>
namespace QuantLib {
const long LecuyerUniformRng::m1 = 2147483563L;
const long LecuyerUniformRng::a1 = 40014L;
const long LecuyerUniformRng::q1 = 53668L;
const long LecuyerUniformRng::r1 = 12211L;
const long LecuyerUniformRng::m2 = 2147483399L;
const long LecuyerUniformRng::a2 = 40692L;
const long LecuyerUniformRng::q2 = 52774L;
const long LecuyerUniformRng::r2 = 3791L;
const int LecuyerUniformRng::bufferSize = 32;
// int(1+m1/bufferSize) = int(1+(m1-1)/bufferSize)
const long LecuyerUniformRng::bufferNormalizer = 67108862L;
const long double LecuyerUniformRng::maxRandom = 1.0-QL_EPSILON;
LecuyerUniformRng::LecuyerUniformRng(long seed)
: buffer(LecuyerUniformRng::bufferSize) {
// Need to prevent seed=0, so use seed=0 to have a "random" seed
temp2 = temp1 = (seed != 0 ? seed : SeedGenerator::instance().get());
// Load the shuffle table (after 8 warm-ups)
for (int j=bufferSize+7; j>=0; j--) {
long k = temp1/q1;
temp1 = a1*(temp1-k*q1)-k*r1;
if (temp1 < 0)
temp1 += m1;
if (j < bufferSize)
buffer[j] = temp1;
}
y = buffer[0];
}
LecuyerUniformRng::sample_type LecuyerUniformRng::next() const {
long k = temp1/q1;
// Compute temp1=(a1*temp1) % m1
// without overflows (Schrage's method)
temp1 = a1*(temp1-k*q1)-k*r1;
if (temp1 < 0)
temp1 += m1;
k = temp2/q2;
// Compute temp2=(a2*temp2) % m2
// without overflows (Schrage's method)
temp2 = a2*(temp2-k*q2)-k*r2;
if (temp2 < 0)
temp2 += m2;
// Will be in the range 0..bufferSize-1
int j = y/bufferNormalizer;
// Here temp1 is shuffled, temp1 and temp2 are
// combined to generate output
y = buffer[j]-temp2;
buffer[j] = temp1;
if (y < 1)
y += m1-1;
double result = y/double(m1);
// users don't expect endpoint values
if (result > maxRandom)
result = (double) maxRandom;
return sample_type(result,1.0);
}
}
|
