/src/quantlib/ql/math/integrals/exponentialintegrals.cpp
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) 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 exponentialintegrals.cpp |
21 | | */ |
22 | | |
23 | | |
24 | | #include <ql/errors.hpp> |
25 | | #include <ql/mathconstants.hpp> |
26 | | #include <ql/math/comparison.hpp> |
27 | | #include <ql/math/integrals/exponentialintegrals.hpp> |
28 | | |
29 | | #include <boost/math/special_functions/sign.hpp> |
30 | | #include <cmath> |
31 | | |
32 | | namespace QuantLib { |
33 | | namespace exponential_integrals_helper { |
34 | | |
35 | | // Reference: |
36 | | // Rowe et al: GALSIM: The modular galaxy image simulation toolkit |
37 | | // https://arxiv.org/abs/1407.7676 |
38 | | |
39 | 0 | Real f(Real x) { |
40 | 0 | const Real x2 = 1.0/(x*x); |
41 | |
|
42 | 0 | return ( |
43 | 0 | 1 + x2*(7.44437068161936700618e2 + x2*(1.96396372895146869801e5 |
44 | 0 | + x2*(2.37750310125431834034e7 + x2*(1.43073403821274636888e9 |
45 | 0 | + x2*(4.33736238870432522765e10 + x2*(6.40533830574022022911e11 |
46 | 0 | + x2*(4.20968180571076940208e12 + x2*(1.00795182980368574617e13 |
47 | 0 | + x2*(4.94816688199951963482e12 - x2*4.94701168645415959931e11))))))))) |
48 | 0 | )/(x *( |
49 | 0 | 1 + x2*(7.46437068161927678031e2 + x2*(1.97865247031583951450e5 |
50 | 0 | + x2*(2.41535670165126845144e7 + x2*(1.47478952192985464958e9 |
51 | 0 | + x2*(4.58595115847765779830e10 + x2*(7.08501308149515401563e11 |
52 | 0 | + x2*(5.06084464593475076774e12 + x2*(1.43468549171581016479e13 |
53 | 0 | + x2*1.11535493509914254097e13)))))))) |
54 | 0 | ) ); |
55 | 0 | } |
56 | | |
57 | 0 | Real g(Real x) { |
58 | 0 | const Real x2 = 1.0/(x*x); |
59 | |
|
60 | 0 | return x2*( |
61 | 0 | 1 + x2*(8.1359520115168615e2 + x2*(2.35239181626478200e5 |
62 | 0 | + x2*(3.12557570795778731e7 + x2*(2.06297595146763354e9 |
63 | 0 | + x2*(6.83052205423625007e10 + x2*(1.09049528450362786e12 |
64 | 0 | + x2*(7.57664583257834349e12 + x2*(1.81004487464664575e13 |
65 | 0 | + x2*(6.43291613143049485e12 - x2*1.36517137670871689e12))))))))) |
66 | 0 | )/( |
67 | 0 | 1 + x2*(8.19595201151451564e2 + x2*(2.40036752835578777e5 |
68 | 0 | + x2*(3.26026661647090822e7 + x2*(2.23355543278099360e9 |
69 | 0 | + x2*(7.87465017341829930e10 + x2*(1.39866710696414565e12 |
70 | 0 | + x2*(1.17164723371736605e13 + x2*(4.01839087307656620e13 |
71 | 0 | + x2*3.99653257887490811e13)))))))) |
72 | 0 | ); |
73 | 0 | } |
74 | | } |
75 | | |
76 | | namespace ExponentialIntegral { |
77 | 0 | Real Si(Real x) { |
78 | 0 | if (x < 0) |
79 | 0 | return -Si(Real(-x)); |
80 | 0 | else if (x <= 4.0) { |
81 | 0 | const Real x2 = x*x; |
82 | |
|
83 | 0 | return x* |
84 | 0 | ( 1 + x2*(-4.54393409816329991e-2 + x2*(1.15457225751016682e-3 |
85 | 0 | + x2*(-1.41018536821330254e-5 + x2*(9.43280809438713025e-8 |
86 | 0 | + x2*(-3.53201978997168357e-10 + x2*(7.08240282274875911e-13 |
87 | 0 | - x2*6.05338212010422477e-16)))))) |
88 | 0 | ) / ( |
89 | 0 | 1 + x2*(1.01162145739225565e-2 + x2*(4.99175116169755106e-5 |
90 | 0 | + x2*(1.55654986308745614e-7 + x2*(3.28067571055789734e-10 |
91 | 0 | + x2*(4.5049097575386581e-13 + x2*3.21107051193712168e-16))))) |
92 | 0 | ); |
93 | 0 | } |
94 | 0 | else { |
95 | 0 | using namespace exponential_integrals_helper; |
96 | 0 | return M_PI_2 - f(x)*std::cos(x) - g(x)*std::sin(x); |
97 | 0 | } |
98 | 0 | } |
99 | | |
100 | 0 | Real Ci(Real x) { |
101 | 0 | QL_REQUIRE(x >= 0, "x < 0 => Ci(x) = Ci(-x) + i*pi"); |
102 | | |
103 | 0 | if (x <= 4.0) { |
104 | 0 | const Real x2 = x*x; |
105 | |
|
106 | 0 | return M_EULER_MASCHERONI + std::log(x) + |
107 | 0 | x2* ( -0.25 + x2*(7.51851524438898291e-3 +x2*(-1.27528342240267686e-4 |
108 | 0 | + x2*(1.05297363846239184e-6 +x2*(-4.68889508144848019e-9 |
109 | 0 | + x2*(1.06480802891189243e-11 - x2*9.93728488857585407e-15))))) |
110 | 0 | ) / ( |
111 | 0 | 1 + x2*(1.1592605689110735e-2 + x2*(6.72126800814254432e-5 |
112 | 0 | + x2*(2.55533277086129636e-7 + x2*(6.97071295760958946e-10 |
113 | 0 | + x2*(1.38536352772778619e-12 + x2*(1.89106054713059759e-15 |
114 | 0 | + x2*1.39759616731376855e-18)))))) |
115 | 0 | ); |
116 | 0 | } |
117 | 0 | else { |
118 | 0 | using namespace exponential_integrals_helper; |
119 | 0 | return f(x)*std::sin(x) - g(x)*std::cos(x); |
120 | 0 | } |
121 | 0 | } |
122 | | |
123 | | std::complex<Real> Ei( |
124 | 0 | const std::complex<Real>& z, const std::complex<Real>& acc) { |
125 | |
|
126 | 0 | if (z.real() == 0.0 && z.imag() == 0.0 |
127 | 0 | && std::numeric_limits<Real>::has_infinity) { |
128 | 0 | return std::complex<Real>(-std::numeric_limits<Real>::infinity()); |
129 | 0 | } |
130 | | |
131 | | |
132 | 0 | constexpr double DIST = 4.5; |
133 | 0 | constexpr double MAX_ERROR = 5.0 * QL_EPSILON; |
134 | |
|
135 | 0 | const Real z_inf = std::log(0.01*QL_MAX_REAL) + std::log(100.0); |
136 | 0 | QL_REQUIRE(z.real() < z_inf, "argument error " << z); |
137 | | |
138 | 0 | const Real z_asym = 2.0 - 1.035*std::log(MAX_ERROR); |
139 | |
|
140 | 0 | using boost::math::sign; |
141 | 0 | const Real abs_z = std::abs(z); |
142 | |
|
143 | 0 | const auto match = [=]( |
144 | 0 | const std::complex<Real>& z1, const std::complex<Real>& z2) |
145 | 0 | -> bool { |
146 | 0 | const std::complex<Real> d = z1 - z2; |
147 | 0 | return std::abs(d.real()) <= MAX_ERROR*std::abs(z1.real()) |
148 | 0 | && std::abs(d.imag()) <= MAX_ERROR*std::abs(z1.imag()); |
149 | 0 | }; |
150 | |
|
151 | 0 | if (z.real() > z_inf) |
152 | 0 | return std::complex<Real>(std::exp(z)/z) + acc; |
153 | | |
154 | 0 | if (abs_z > 1.1*z_asym) { |
155 | 0 | std::complex<Real> ei = acc + std::complex<Real>(Real(0.0), sign(z.imag())*M_PI); |
156 | 0 | std::complex<Real> s(std::exp(z)/z); |
157 | 0 | for (Size i=1; i <= std::floor(abs_z)+1; ++i) { |
158 | 0 | if (match(ei+s, ei)) |
159 | 0 | return ei+s; |
160 | 0 | ei += s; |
161 | 0 | s *= Real(i)/z; |
162 | 0 | } |
163 | 0 | QL_FAIL("series conversion issue for Ei(" << z << ")"); |
164 | 0 | } |
165 | | |
166 | 0 | if (abs_z > DIST && (z.real() < 0 || std::abs(z.imag()) > DIST)) { |
167 | 0 | std::complex<Real> ei(0.0); |
168 | 0 | for (Size k = 47; k >=1; --k) { |
169 | 0 | ei = - Real(k*k)/(2.0*k + 1.0 - z + ei); |
170 | 0 | } |
171 | 0 | return (acc + std::complex<Real>(0.0, sign(z.imag())*M_PI)) |
172 | 0 | - std::exp(z)/ (1.0 - z + ei); |
173 | 0 | } |
174 | | |
175 | 0 | std::complex<Real> s(0.0), sn=z; |
176 | 0 | Real nn=1.0; |
177 | |
|
178 | 0 | Size n; |
179 | 0 | for (n=2; n < 1000 && s+sn*nn != s; ++n) { |
180 | 0 | s+=sn*nn; |
181 | |
|
182 | 0 | if ((n & 1) != 0U) |
183 | 0 | nn += 1/(2.0*(n/2) + 1); // NOLINT(bugprone-integer-division) |
184 | |
|
185 | 0 | sn *= -z / Real(2*n); |
186 | 0 | } |
187 | |
|
188 | 0 | QL_REQUIRE(n < 1000, "series conversion issue for Ei(" << z << ")"); |
189 | | |
190 | 0 | const std::complex<Real> r |
191 | 0 | = (M_EULER_MASCHERONI + acc) + std::log(z) + std::exp(0.5*z)*s; |
192 | |
|
193 | 0 | if (z.imag() != Real(0.0)) |
194 | 0 | return r; |
195 | 0 | else |
196 | 0 | return std::complex<Real>(r.real(), acc.imag()); |
197 | 0 | } |
198 | | |
199 | 0 | std::complex<Real> Ei(const std::complex<Real>&z) { |
200 | 0 | return Ei(z, std::complex<Real>(0.0)); |
201 | 0 | } |
202 | | |
203 | 0 | std::complex<Real> E1(const std::complex<Real>& z) { |
204 | 0 | if (z.imag() < 0.0) { |
205 | 0 | return -Ei(-z, std::complex<Real>(0.0, -M_PI)); |
206 | 0 | } |
207 | 0 | else if (z.imag() > 0.0 || z.real() < 0.0) { |
208 | 0 | return -Ei(-z, std::complex<Real>(0.0, M_PI)); |
209 | 0 | } |
210 | 0 | else { |
211 | 0 | return -Ei(-z); |
212 | 0 | } |
213 | 0 | } |
214 | | |
215 | | // Reference: |
216 | | // https://functions.wolfram.com/GammaBetaErf/ExpIntegralEi/introductions/ExpIntegrals/ShowAll.html |
217 | 0 | std::complex<Real> Si(const std::complex<Real>& z) { |
218 | 0 | if (std::abs(z) <= 0.2) { |
219 | 0 | std::complex<Real> s(0), nn(z); |
220 | 0 | Size k; |
221 | 0 | for (k=2; k < 100 && s != s+nn; ++k) { |
222 | 0 | s += nn; |
223 | 0 | nn *= -z*z/((2.0*k-2)*(2*k-1)*(2*k-1))*(2.0*k-3); |
224 | 0 | } |
225 | 0 | QL_REQUIRE(k < 100, "series conversion issue for Si(" << z << ")"); |
226 | | |
227 | 0 | return s; |
228 | 0 | } |
229 | 0 | else { |
230 | 0 | const std::complex<Real> i(0.0, 1.0); |
231 | 0 | return 0.5*i*(E1(-i*z) - E1(i*z) |
232 | 0 | - std::complex<Real>(0.0, ((z.real() >= 0 && z.imag() >= 0) |
233 | 0 | || (z.real() > 0 && z.imag() < 0) )? M_PI : -M_PI)); |
234 | 0 | } |
235 | 0 | } |
236 | | |
237 | 0 | std::complex<Real> Ci(const std::complex<Real>& z) { |
238 | 0 | const std::complex<Real> i(0.0, 1.0); |
239 | |
|
240 | 0 | std::complex<Real> acc(0.0); |
241 | 0 | if (z.real() < 0.0 && z.imag() >= 0.0) |
242 | 0 | acc.imag(M_PI); |
243 | 0 | else if (z.real() <= 0.0 && z.imag() <= 0.0) |
244 | 0 | acc.imag(-M_PI); |
245 | |
|
246 | 0 | return -0.5*(E1(-i*z) + E1(i*z)) + acc; |
247 | 0 | } |
248 | | } |
249 | | } |