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);}
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); }
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); }
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) {
53 /GammaFunction().value(1.0+
nu);
58 while (std::abs(B_k*=Y/(k*(k+
nu)))>std::abs(sum)*
QL_EPSILON) {
60 QL_REQUIRE(++k < 1000,
"max iterations exceeded");
62 return sum * W<T>().weightSmallX(x);
65 Real na_k=1.0, sign=1.0;
68 T s1=
T(1.0), s2=
T(1.0);
69 for (
Size k=1; k < 30; ++k) {
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;
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);
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)) /
97 QL_REQUIRE(x >= 0.0,
"negative argument requires complex version of "
98 "modifiedBesselFunction");
99 return modifiedBesselFunction_i_impl<Real, Unweighted>(
nu, x);
103 const std::complex<Real> &z) {
104 if (z.imag() == 0.0 && z.real() >= 0.0)
107 return modifiedBesselFunction_i_impl<
108 std::complex<Real>, Unweighted>(
nu, z);
112 return modifiedBesselFunction_k_impl<Real, Unweighted>(
nu, x);
116 const std::complex<Real> &z) {
117 if (z.imag() == 0.0 && z.real() >= 0.0)
120 return modifiedBesselFunction_k_impl<
121 std::complex<Real>, Unweighted>(
nu, z);
125 QL_REQUIRE(x >= 0.0,
"negative argument requires complex version of "
126 "modifiedBesselFunction");
127 return modifiedBesselFunction_i_impl<Real, ExponentiallyWeighted>(
132 Real nu,
const std::complex<Real> &z) {
134 if (z.imag() == 0.0 && z.real() >= 0.0)
135 return std::complex<Real>(
138 return modifiedBesselFunction_i_impl<
139 std::complex<Real>, ExponentiallyWeighted>(
nu, z);
143 return modifiedBesselFunction_k_impl<Real, ExponentiallyWeighted>(
148 Real nu,
const std::complex<Real> &z) {
150 if (z.imag() == 0.0 && z.real() >= 0.0)
151 return std::complex<Real>(
154 return modifiedBesselFunction_k_impl<
155 std::complex<Real>, ExponentiallyWeighted>(
nu, z);
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
std::size_t Size
size of a container
modified Bessel functions of first and second kind
Real modifiedBesselFunction_k_exponentiallyWeighted(Real nu, Real x)
Real modifiedBesselFunction_k(Real nu, Real x)
Real modifiedBesselFunction_i(Real nu, Real x)
Real modifiedBesselFunction_i_exponentiallyWeighted(Real nu, Real x)