1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2014 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 modifiedbessel.cpp
21 \brief modified Bessel functions of first and second kind
22*/
23
24#include <ql/math/modifiedbessel.hpp>
25#include <ql/math/distributions/gammadistribution.hpp>
26
27#include <cmath>
28
29namespace QuantLib {
30
31 namespace {
32
33 template <class T> struct I {};
34 template <> struct I<Real> { Real value() { return 0.0;} };
35 template <> struct I<std::complex<Real> > {
36 std::complex<Real> value() { return std::complex<Real>(0.0,1.0);}
37 };
38 template <class T> struct Unweighted {
39 T weightSmallX(const T& x) { return 1.0; }
40 T weight1LargeX(const T& x) { return std::exp(x); }
41 T weight2LargeX(const T& x) { return std::exp(-x); }
42 };
43 template <class T> struct ExponentiallyWeighted {
44 T weightSmallX(const T& x) { return std::exp(-x); }
45 T weight1LargeX(const T& x) { return 1.0; }
46 T weight2LargeX(const T& x) { return std::exp(-2.0*x); }
47 };
48
49 template <class T, template <class> class W>
50 T modifiedBesselFunction_i_impl(Real nu, const T& x) {
51 if (std::abs(x) < 13.0) {
52 const T alpha = std::pow(0.5*x, nu)
53 /GammaFunction().value(x: 1.0+nu);
54 const T Y = 0.25*x*x;
55 Size k=1;
56 T sum=alpha, B_k=alpha;
57
58 while (std::abs(B_k*=Y/(k*(k+nu)))>std::abs(sum)*QL_EPSILON) {
59 sum += B_k;
60 QL_REQUIRE(++k < 1000, "max iterations exceeded");
61 }
62 return sum * W<T>().weightSmallX(x);
63 }
64 else {
65 Real na_k=1.0, sign=1.0;
66 T da_k=T(1.0);
67
68 T s1=T(1.0), s2=T(1.0);
69 for (Size k=1; k < 30; ++k) {
70 sign*=-1;
71 na_k *= (4.0 * nu * nu -
72 (2.0 * static_cast<Real>(k) - 1.0) *
73 (2.0 * static_cast<Real>(k) - 1.0));
74 da_k *= (8.0 * k) * x;
75 const T a_k = na_k/da_k;
76
77 s2+=a_k;
78 s1+=sign*a_k;
79 }
80
81 const T i = I<T>().value();
82 return 1.0 / std::sqrt(2 * M_PI * x) *
83 (W<T>().weight1LargeX(x) * s1 +
84 i * std::exp(i * nu * M_PI) * W<T>().weight2LargeX(x) * s2);
85 }
86 }
87
88 template <class T, template <class> class W>
89 T modifiedBesselFunction_k_impl(Real nu, const T& x) {
90 return M_PI_2 * (modifiedBesselFunction_i_impl<T,W>(-nu, x) -
91 modifiedBesselFunction_i_impl<T,W>(nu, x)) /
92 std::sin(M_PI * nu);
93 }
94 }
95
96 Real modifiedBesselFunction_i(Real nu, Real x) {
97 QL_REQUIRE(x >= 0.0, "negative argument requires complex version of "
98 "modifiedBesselFunction");
99 return modifiedBesselFunction_i_impl<Real, Unweighted>(nu, x);
100 }
101
102 std::complex<Real> modifiedBesselFunction_i(Real nu,
103 const std::complex<Real> &z) {
104 if (z.imag() == 0.0 && z.real() >= 0.0)
105 return std::complex<Real>(modifiedBesselFunction_i(nu, x: z.real()));
106
107 return modifiedBesselFunction_i_impl<
108 std::complex<Real>, Unweighted>(nu, x: z);
109 }
110
111 Real modifiedBesselFunction_k(Real nu, Real x) {
112 return modifiedBesselFunction_k_impl<Real, Unweighted>(nu, x);
113 }
114
115 std::complex<Real> modifiedBesselFunction_k(Real nu,
116 const std::complex<Real> &z) {
117 if (z.imag() == 0.0 && z.real() >= 0.0)
118 return std::complex<Real>(modifiedBesselFunction_k(nu, x: z.real()));
119
120 return modifiedBesselFunction_k_impl<
121 std::complex<Real>, Unweighted>(nu, x: z);
122 }
123
124 Real modifiedBesselFunction_i_exponentiallyWeighted(Real nu, Real x) {
125 QL_REQUIRE(x >= 0.0, "negative argument requires complex version of "
126 "modifiedBesselFunction");
127 return modifiedBesselFunction_i_impl<Real, ExponentiallyWeighted>(
128 nu, x);
129 }
130
131 std::complex<Real> modifiedBesselFunction_i_exponentiallyWeighted(
132 Real nu, const std::complex<Real> &z) {
133
134 if (z.imag() == 0.0 && z.real() >= 0.0)
135 return std::complex<Real>(
136 modifiedBesselFunction_i_exponentiallyWeighted(nu, x: z.real()));
137
138 return modifiedBesselFunction_i_impl<
139 std::complex<Real>, ExponentiallyWeighted>(nu, x: z);
140 }
141
142 Real modifiedBesselFunction_k_exponentiallyWeighted(Real nu, Real x) {
143 return modifiedBesselFunction_k_impl<Real, ExponentiallyWeighted>(
144 nu, x);
145 }
146
147 std::complex<Real> modifiedBesselFunction_k_exponentiallyWeighted(
148 Real nu, const std::complex<Real> &z) {
149
150 if (z.imag() == 0.0 && z.real() >= 0.0)
151 return std::complex<Real>(
152 modifiedBesselFunction_k_exponentiallyWeighted(nu, x: z.real()));
153
154 return modifiedBesselFunction_k_impl<
155 std::complex<Real>, ExponentiallyWeighted>(nu, x: z);
156 }
157
158}
159

source code of quantlib/ql/math/modifiedbessel.cpp