QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
gausslaguerrecosinepolynomial.hpp
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) 2020 Klaus Spanderen
5
6 This file is part of QuantLib, a free-software/open-source library
7 for financial quantitative analysts and developers - http://quantlib.org/
8
9 QuantLib is free software: you can redistribute it and/or modify it
10 under the terms of the QuantLib license. You should have received a
11 copy of the license along with this program; if not, please email
12 <quantlib-dev@lists.sf.net>. The license is also available online at
13 <http://quantlib.org/license.shtml>.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the license for more details.
18*/
19
20/*! \file gausslaguerrecosinepolynomial.hpp
21 \brief Laguerre-Cosine Gaussian quadrature
22*/
23
24#ifndef quantlib_gauss_laguerre_cosine_polynomial_hpp
25#define quantlib_gauss_laguerre_cosine_polynomial_hpp
26
29#include <cmath>
30
31namespace QuantLib {
32
33 template <class mp_real>
35 : public MomentBasedGaussianPolynomial<mp_real> {
36 public:
38 : u_(u) { }
39
40 protected:
41 virtual mp_real m0() const = 0;
42 virtual mp_real m1() const = 0;
43
44 mp_real moment_(Size n) const { //NOLINT(bugprone-virtual-near-miss)
45 if (m_.size() <= n)
46 m_.resize(n+1, std::numeric_limits<mp_real>::quiet_NaN());
47
48 if (std::isnan(m_[n])) {
49 if (n == 0)
50 m_[0] = m0();
51 else if (n == 1)
52 m_[1] = m1();
53 else
54 m_[n] = (2*n*moment_(n-1)
55 - n*(n-1)*moment_(n-2))/(1+u_*u_);
56 }
57
58 return m_[n];
59 }
60 mp_real fact(Size n) const {
61 if (f_.size() <= n)
62 f_.resize(n+1, std::numeric_limits<mp_real>::quiet_NaN());
63
64 if (std::isnan(f_[n])) {
65 if (n == 0)
66 f_[0] = 1.0;
67 else
68 f_[n] = n*fact(n-1);
69 }
70 return f_[n];
71
72 }
73 const Real u_;
74
75 private:
76 mutable std::vector<mp_real> m_, f_;
77 };
78
79 //! Gauss-Laguerre Cosine integration
80
81 /*! This class performs a 1-dimensional Gauss-Laguerre-Cosine integration.
82 \f[
83 \int_{0}^{\inf} f(x) \mathrm{d}x
84 \f]
85 The weighting function is
86 \f[
87 w(x;u)=e^{-x}*\cos{u*x}
88 \f]
89 */
90 template <class mp_real>
92 : public GaussLaguerreTrigonometricBase<mp_real> {
93 public:
96 m0_(1.0+1.0/(1.0+u*u)) { }
97
98 mp_real moment(Size n) const override { return (this->moment_(n) + this->fact(n)) / m0_; }
99 Real w(Real x) const override { return std::exp(-x) * (1 + std::cos(this->u_ * x)) / m0_; }
100
101 protected:
102 mp_real m0() const override { return 1 / (1 + this->u_ * this->u_); }
103 mp_real m1() const override {
104 return (1 - this->u_*this->u_) / squared(1 + this->u_*this->u_);
105 }
106
107 private:
108 const Real m0_;
109 };
110
111 //! Gauss-Laguerre Sine integration
112
113 /*! This class performs a 1-dimensional Gauss-Laguerre-Cosine integration.
114 \f[
115 \int_{0}^{\inf} f(x) \mathrm{d}x
116 \f]
117 The weighting function is
118 \f[
119 w(x;u)=e^{-x}*\sin{u*x}
120 \f]
121 */
122 template <class mp_real>
124 : public GaussLaguerreTrigonometricBase<mp_real> {
125 public:
127 : GaussLaguerreTrigonometricBase<mp_real>(u),
128 m0_(1.0+u/(1.0+u*u)) { }
129
130 mp_real moment(Size n) const override { return (this->moment_(n) + this->fact(n)) / m0_; }
131 Real w(Real x) const override { return std::exp(-x) * (1 + std::sin(this->u_ * x)) / m0_; }
132
133 protected:
134 mp_real m0() const override { return this->u_ / (1 + this->u_ * this->u_); }
135 mp_real m1() const override {
136 return 2*this->u_ / squared(1 + this->u_*this->u_);
137 }
138
139 private:
140 const Real m0_;
141 };
142}
143
144#endif
QL_REAL Real
real number
Definition: types.hpp:50
std::size_t Size
size of a container
Definition: types.hpp:58
functionals and combinators not included in the STL
Gaussian quadrature defined by the moments of the distribution.
Definition: any.hpp:35
T squared(T x)
Definition: functional.hpp:37