/src/quantlib/ql/math/modifiedbessel.cpp

Source (jump to first uncovered line)
/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*
 Copyright (C) 2014 Klaus Spanderen

 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.
*/

/*! \file modifiedbessel.cpp
    \brief modified Bessel functions of first and second kind
*/

#include <ql/math/modifiedbessel.hpp>
#include <ql/math/distributions/gammadistribution.hpp>

#include <cmath>

namespace QuantLib {

    namespace {

        template <class T>  struct I {};
        template <> struct I<Real> { Real value() { return 0.0;} };
        template <> struct I<std::complex<Real> > {
            std::complex<Real> value() { return std::complex<Real>(0.0,1.0);}
        };
        template <class T> struct Unweighted {
            T weightSmallX(const T& x) { return 1.0; }
            T weight1LargeX(const T& x) { return std::exp(x); }
            T weight2LargeX(const T& x) { return std::exp(-x); }
        };
        template <class T> struct ExponentiallyWeighted {
            T weightSmallX(const T& x) { return std::exp(-x); }
            T weight1LargeX(const T& x) { return 1.0; }
            T weight2LargeX(const T& x) { return std::exp(-2.0*x); }
        };

        template <class T, template <class> class W>
        T modifiedBesselFunction_i_impl(Real nu, const T& x) {
            if (std::abs(x) < 13.0) {
                const T alpha = std::pow(0.5*x, nu)
                               /GammaFunction().value(1.0+nu);
                const T Y = 0.25*x*x;
                Size k=1;
                T sum=alpha, B_k=alpha;

                while (std::abs(B_k*=Y/(k*(k+nu)))>std::abs(sum)*QL_EPSILON) {
                    sum += B_k;
                    QL_REQUIRE(++k < 1000, "max iterations exceeded");
                }
                return sum * W<T>().weightSmallX(x);
            }
            else {
                Real na_k=1.0, sign=1.0;
                T da_k=T(1.0);

                T s1=T(1.0), s2=T(1.0);
                for (Size k=1; k < 30; ++k) {
                    sign*=-1;
                    na_k *= (4.0 * nu * nu -
                             (2.0 * static_cast<Real>(k) - 1.0) *
                                 (2.0 * static_cast<Real>(k) - 1.0));
                    da_k *= (8.0 * k) * x;
                    const T a_k = na_k/da_k;

                    s2+=a_k;
                    s1+=sign*a_k;
                }

                const T i = I<T>().value();
                return 1.0 / std::sqrt(2 * M_PI * x) *
                    (W<T>().weight1LargeX(x) * s1 +
                     i * std::exp(i * nu * M_PI) * W<T>().weight2LargeX(x) * s2);
            }
        }

        template <class T, template <class> class W>
        T modifiedBesselFunction_k_impl(Real nu, const T& x) {
            return M_PI_2 * (modifiedBesselFunction_i_impl<T,W>(-nu, x) -
                             modifiedBesselFunction_i_impl<T,W>(nu, x)) /
                             std::sin(M_PI * nu);
        }
    }

    Real modifiedBesselFunction_i(Real nu, Real x) {
        QL_REQUIRE(x >= 0.0, "negative argument requires complex version of "
                             "modifiedBesselFunction");
        return modifiedBesselFunction_i_impl<Real, Unweighted>(nu, x);
    }

    std::complex<Real> modifiedBesselFunction_i(Real nu,
                                                const std::complex<Real> &z) {
        if (z.imag() == 0.0 && z.real() >= 0.0)
            return std::complex<Real>(modifiedBesselFunction_i(nu, z.real()));

        return modifiedBesselFunction_i_impl<
            std::complex<Real>, Unweighted>(nu, z);
    }

    Real modifiedBesselFunction_k(Real nu, Real x) {
        return modifiedBesselFunction_k_impl<Real, Unweighted>(nu, x);
    }

    std::complex<Real> modifiedBesselFunction_k(Real nu,
                                                const std::complex<Real> &z) {
        if (z.imag() == 0.0 && z.real() >= 0.0)
            return std::complex<Real>(modifiedBesselFunction_k(nu, z.real()));

        return modifiedBesselFunction_k_impl<
            std::complex<Real>, Unweighted>(nu, z);
    }

    Real modifiedBesselFunction_i_exponentiallyWeighted(Real nu, Real x) {
        QL_REQUIRE(x >= 0.0, "negative argument requires complex version of "
                             "modifiedBesselFunction");
        return modifiedBesselFunction_i_impl<Real, ExponentiallyWeighted>(
            nu, x);
    }

    std::complex<Real> modifiedBesselFunction_i_exponentiallyWeighted(
        Real nu, const std::complex<Real> &z) {

        if (z.imag() == 0.0 && z.real() >= 0.0)
            return std::complex<Real>(
                modifiedBesselFunction_i_exponentiallyWeighted(nu, z.real()));

        return modifiedBesselFunction_i_impl<
            std::complex<Real>, ExponentiallyWeighted>(nu, z);
    }

    Real modifiedBesselFunction_k_exponentiallyWeighted(Real nu, Real x) {
        return modifiedBesselFunction_k_impl<Real, ExponentiallyWeighted>(
            nu, x);
    }

    std::complex<Real> modifiedBesselFunction_k_exponentiallyWeighted(
        Real nu, const std::complex<Real> &z) {

        if (z.imag() == 0.0 && z.real() >= 0.0)
            return std::complex<Real>(
                modifiedBesselFunction_k_exponentiallyWeighted(nu, z.real()));

        return modifiedBesselFunction_k_impl<
            std::complex<Real>, ExponentiallyWeighted>(nu, z);
    }

}

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

Coverage Report

Created: 2025-08-05 06:45