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/integrals/exponentialintegrals.hpp>
27
28#include <cmath>
29
30#ifndef M_EULER_MASCHERONI
31 #define M_EULER_MASCHERONI 0.5772156649015328606065121
32#endif
33
34namespace QuantLib {
35 namespace exponential_integrals_helper {
36
37 // Reference:
38 // Rowe et al: GALSIM: The modular galaxy image simulation toolkit
39 // https://arxiv.org/abs/1407.7676
40
41 Real f(Real x) {
42 const Real x2 = 1.0/(x*x);
43
44 return (
45 1 + x2*(7.44437068161936700618e2 + x2*(1.96396372895146869801e5
46 + x2*(2.37750310125431834034e7 + x2*(1.43073403821274636888e9
47 + x2*(4.33736238870432522765e10 + x2*(6.40533830574022022911e11
48 + x2*(4.20968180571076940208e12 + x2*(1.00795182980368574617e13
49 + x2*(4.94816688199951963482e12 - x2*4.94701168645415959931e11)))))))))
50 )/(x *(
51 1 + x2*(7.46437068161927678031e2 + x2*(1.97865247031583951450e5
52 + x2*(2.41535670165126845144e7 + x2*(1.47478952192985464958e9
53 + x2*(4.58595115847765779830e10 + x2*(7.08501308149515401563e11
54 + x2*(5.06084464593475076774e12 + x2*(1.43468549171581016479e13
55 + x2*1.11535493509914254097e13))))))))
56 ) );
57 }
58
59 Real g(Real x) {
60 const Real x2 = 1.0/(x*x);
61
62 return x2*(
63 1 + x2*(8.1359520115168615e2 + x2*(2.35239181626478200e5
64 + x2*(3.12557570795778731e7 + x2*(2.06297595146763354e9
65 + x2*(6.83052205423625007e10 + x2*(1.09049528450362786e12
66 + x2*(7.57664583257834349e12 + x2*(1.81004487464664575e13
67 + x2*(6.43291613143049485e12 - x2*1.36517137670871689e12)))))))))
68 )/(
69 1 + x2*(8.19595201151451564e2 + x2*(2.40036752835578777e5
70 + x2*(3.26026661647090822e7 + x2*(2.23355543278099360e9
71 + x2*(7.87465017341829930e10 + x2*(1.39866710696414565e12
72 + x2*(1.17164723371736605e13 + x2*(4.01839087307656620e13
73 + x2*3.99653257887490811e13))))))))
74 );
75 }
76 }
77
78 namespace ExponentialIntegral {
79 Real Si(Real x) {
80 if (x < 0)
81 return -Si(x: Real(-x));
82 else if (x <= 4.0) {
83 const Real x2 = x*x;
84
85 return x*
86 ( 1 + x2*(-4.54393409816329991e-2 + x2*(1.15457225751016682e-3
87 + x2*(-1.41018536821330254e-5 + x2*(9.43280809438713025e-8
88 + x2*(-3.53201978997168357e-10 + x2*(7.08240282274875911e-13
89 - x2*6.05338212010422477e-16))))))
90 ) / (
91 1 + x2*(1.01162145739225565e-2 + x2*(4.99175116169755106e-5
92 + x2*(1.55654986308745614e-7 + x2*(3.28067571055789734e-10
93 + x2*(4.5049097575386581e-13 + x2*3.21107051193712168e-16)))))
94 );
95 }
96 else {
97 using namespace exponential_integrals_helper;
98 return M_PI_2 - f(x)*std::cos(x: x) - g(x)*std::sin(x: x);
99 }
100 }
101
102 Real Ci(Real x) {
103 QL_REQUIRE(x >= 0, "x < 0 => Ci(x) = Ci(-x) + i*pi");
104
105 if (x <= 4.0) {
106 const Real x2 = x*x;
107
108 return M_EULER_MASCHERONI + std::log(x: x) +
109 x2* ( -0.25 + x2*(7.51851524438898291e-3 +x2*(-1.27528342240267686e-4
110 + x2*(1.05297363846239184e-6 +x2*(-4.68889508144848019e-9
111 + x2*(1.06480802891189243e-11 - x2*9.93728488857585407e-15)))))
112 ) / (
113 1 + x2*(1.1592605689110735e-2 + x2*(6.72126800814254432e-5
114 + x2*(2.55533277086129636e-7 + x2*(6.97071295760958946e-10
115 + x2*(1.38536352772778619e-12 + x2*(1.89106054713059759e-15
116 + x2*1.39759616731376855e-18))))))
117 );
118 }
119 else {
120 using namespace exponential_integrals_helper;
121 return f(x)*std::sin(x: x) - g(x)*std::cos(x: x);
122 }
123 }
124
125
126 std::complex<Real> E1(std::complex<Real> z) {
127 QL_REQUIRE(std::abs(z) <= 25.0, "Insufficient precision for |z| > 25.0");
128
129 std::complex<Real> s(0.0), sn(-z);
130
131 Size n;
132 for (n=2; n < 1000 && s + sn/Real(n-1) != s; ++n) {
133 s+=sn/Real(n-1);
134 sn *= -z/Real(n);
135 }
136
137 QL_REQUIRE(n < 1000, "series conversion issue");
138
139 return -M_EULER_MASCHERONI - std::log(z: z) -s;
140 }
141
142 std::complex<Real> Ei(std::complex<Real> z) {
143 QL_REQUIRE(std::abs(z) <= 25.0, "Insufficient precision for |z| > 25.0");
144
145 std::complex<Real> s(0.0), sn=z;
146
147 Real nn=1.0;
148
149 Size n;
150 for (n=2; n < 1000 && s+sn*nn != s; ++n) {
151 s+=sn*nn;
152
153 if ((n & 1) != 0U)
154 nn += 1/(2.0*(n/2) + 1); // NOLINT(bugprone-integer-division)
155
156 sn *= -z / Real(2*n);
157 }
158
159 QL_REQUIRE(n < 1000, "series conversion issue");
160
161 return M_EULER_MASCHERONI + std::log(z: z) + std::exp(z: 0.5*z)*s;
162 }
163
164 // Reference:
165 // https://functions.wolfram.com/GammaBetaErf/ExpIntegralEi/introductions/ExpIntegrals/ShowAll.html
166 std::complex<Real> Si(std::complex<Real> z) {
167 const std::complex<Real> i(0.0, 1.0);
168
169 return 0.25*i*(2.0*(Ei(z: -i*z) - Ei(z: i*z))
170 + std::log(z: i/z) - std::log(z: -i/z) - std::log(z: -i*z)
171 + std::log(z: i*z));
172 }
173
174 std::complex<Real> Ci(std::complex<Real> z) {
175 const std::complex<Real> i(0.0, 1.0);
176
177 return 0.25*(2.0*(Ei(z: -i*z) + Ei(z: i*z))
178 + std::log(z: i/z) + std::log(z: -i/z) - std::log(z: -i*z)
179 - std::log(z: i*z)) + std::log(z: z);
180 }
181 }
182}
183

source code of quantlib/ql/math/integrals/exponentialintegrals.cpp

Morty Proxy This is a proxified and sanitized view of the page, visit original site.