QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
faurersg.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) 2004 Ferdinando Ametrano
5 Copyright (C) 2004 Gianni Piolanti
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
23
24namespace QuantLib {
25
26 FaureRsg::FaureRsg(Size dimensionality)
27 : dimensionality_(dimensionality),
28// sequenceCounter_(0),
29 sequence_(std::vector<Real> (dimensionality), 1.0),
30 integerSequence_(dimensionality, 0) {
31
32 QL_REQUIRE(dimensionality>0,
33 "dimensionality must be greater than 0");
34
35 // base is the lowest prime number >= dimensionality_
36 Size i, j, k=1;
37 base_=2;
38 while (base_<dimensionality_) {
40 k++;
41 }
42
43
44 mbit_=(Size)(std::log((double)std::numeric_limits<long int>::max())/
45 std::log((double)base_));
46 gray_ = std::vector<std::vector<long int> >(dimensionality_,
47 std::vector<long int>(mbit_+1, 0));
48 bary_ = std::vector<long int>(mbit_+1, 0);
49
50
51 //setMatrixValues();
52 powBase_ = std::vector<std::vector<long int> >(mbit_,
53 std::vector<long int>(2*base_-1, 0));
54 powBase_[mbit_-1][base_] = 1;
55 for (int i2=mbit_-2; i2>=0; --i2)
56 powBase_[i2][base_] = powBase_[i2+1][base_] * base_;
57 for (int ii=0; ii<(int)mbit_; ii++) {
58 for (int j1=base_+1; j1<2*(int)base_-1; j1++ )
59 powBase_[ii][j1] = powBase_[ii][j1-1] + powBase_[ii][base_];
60 for (int j2=base_-1; j2>=0; --j2)
61 powBase_[ii][j2] = powBase_[ii][j2+1] - powBase_[ii][base_];
62 }
63
64 addOne_.resize(base_);
65 for (j=0; j<base_ ; j++)
66 addOne_[j] = (j+1) % base_;
67
68
69 //setPascalMatrix();
70 for (k=0; k<mbit_; k++) {
71 std::vector<std::vector<long int> > mm(dimensionality_+1,
72 std::vector<long int>(k+1, 0));
73 pascal3D.push_back(mm);
74 pascal3D[k][0][k] = 1;
75 pascal3D[k][1][0] = 1;
76 pascal3D[k][1][k] = 1;
77 }
78
79 long int p1, p2;
80 for (k=2; k<mbit_ ; k++) {
81 for (i=1; i<k ; i++) {
82 p1 = pascal3D[k-1][1][i-1];
83 p2 = pascal3D[k-1][1][i];
84 pascal3D[k][1][i] = (p1+p2) % base_;
85 }
86 }
87
88 long int fact = 1, diag;
89 for (j=2; j<dimensionality_; j++) {
90 for (long int kk=mbit_-1; kk>=0 ; --kk) {
91 diag = mbit_ - kk - 1;
92 if (diag==0)
93 fact = 1;
94 else
95 fact = (fact*j) % base_;
96 for (long int ii=0; ii<=kk; ii++)
97 pascal3D[diag+ii][j][ii] = (fact*
98 pascal3D[diag+ii][1][ii]) % base_;
99 }
100 }
101
102
103 normalizationFactor_ = (double)base_ * (double)powBase_[0][base_];
104 // std::cout << IntegerFormatter::toString(dimensionality_) << ", " ;
105 // std::cout << IntegerFormatter::toString(normalizationFactor_);
106 // std::cout << std::endl;
107 }
108
110 // sequenceCounter_++;
111
112 Size bit = 0;
113 bary_[bit] = addOne_[bary_[bit]];
114 while (bary_[bit] == 0) {
115 bit++;
116 bary_[bit] = addOne_[bary_[bit]];
117 };
118 QL_REQUIRE(bit != mbit_,
119 "Error processing Faure sequence." );
120
121 long int tmp, g1, g2;
122 for (Size i=0; i<dimensionality_ ; i++) {
123 for (Size j=0; j<=bit ; j++) {
124 tmp = gray_[i][j];
125 gray_[i][j] = (pascal3D[bit][i][j] + tmp) % base_;
126 g1 = gray_[i][j];
127 g2 = base_ - 1 + g1 - tmp;
128 integerSequence_[i] += powBase_[j][g2];
129 }
130 }
131 }
132
133}
134
135
double normalizationFactor_
Definition: faurersg.hpp:77
std::vector< std::vector< long int > > gray_
Definition: faurersg.hpp:72
FaureRsg(Size dimensionality)
Definition: faurersg.cpp:26
std::vector< long int > addOne_
Definition: faurersg.hpp:75
std::vector< long int > integerSequence_
Definition: faurersg.hpp:70
std::vector< std::vector< long int > > powBase_
Definition: faurersg.hpp:74
std::vector< std::vector< std::vector< long int > > > pascal3D
Definition: faurersg.hpp:76
void generateNextIntSequence() const
Definition: faurersg.cpp:109
std::vector< long int > bary_
Definition: faurersg.hpp:71
static BigNatural get(Size absoluteIndex)
Get and store one after another.
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
Faure low-discrepancy sequence generator.
QL_REAL Real
real number
Definition: types.hpp:50
std::size_t Size
size of a container
Definition: types.hpp:58
Definition: any.hpp:35
STL namespace.
Prime numbers calculator.