1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2003 Ferdinando Ametrano
5 Copyright (C) 2014 Klaus Spanderen
6 Copyright (C) 2015 Johannes Göttker-Schnetmann
7
8 This file is part of QuantLib, a free-software/open-source library
9 for financial quantitative analysts and developers - http://quantlib.org/
10
11 QuantLib is free software: you can redistribute it and/or modify it
12 under the terms of the QuantLib license. You should have received a
13 copy of the license along with this program; if not, please email
14 <quantlib-dev@lists.sf.net>. The license is also available online at
15 <http://quantlib.org/license.shtml>.
16
17 This program is distributed in the hope that it will be useful, but WITHOUT
18 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 FOR A PARTICULAR PURPOSE. See the license for more details.
20*/
21
22#include "functions.hpp"
23#include "utilities.hpp"
24#include <ql/math/comparison.hpp>
25#include <ql/math/factorial.hpp>
26#include <ql/math/distributions/gammadistribution.hpp>
27#include <ql/math/modifiedbessel.hpp>
28
29using namespace QuantLib;
30using namespace boost::unit_test_framework;
31
32using std::exp;
33
34void FunctionsTest::testFactorial() {
35
36 BOOST_TEST_MESSAGE("Testing factorial numbers...");
37
38 Real expected = 1.0;
39 Real calculated = Factorial::get(n: 0);
40 if (calculated!=expected)
41 BOOST_FAIL("Factorial(0) = " << calculated);
42
43 for (Natural i=1; i<171; ++i) {
44 expected *= i;
45 calculated = Factorial::get(n: i);
46 if (std::fabs(x: calculated-expected)/expected > 1.0e-9)
47 BOOST_FAIL("Factorial(" << i << ")" <<
48 std::setprecision(16) << std::scientific <<
49 "\n calculated: " << calculated <<
50 "\n expected: " << expected <<
51 "\n rel. error: " <<
52 std::fabs(calculated-expected)/expected);
53 }
54}
55
56void FunctionsTest::testGammaFunction() {
57
58 BOOST_TEST_MESSAGE("Testing Gamma function...");
59
60 Real expected = 0.0;
61 Real calculated = GammaFunction().logValue(x: 1);
62 if (std::fabs(x: calculated) > 1.0e-15)
63 BOOST_ERROR("GammaFunction(1)\n"
64 << std::setprecision(16) << std::scientific
65 << " calculated: " << calculated << "\n"
66 << " expected: " << expected);
67
68 for (Size i=2; i<9000; i++) {
69 expected += std::log(x: Real(i));
70 calculated = GammaFunction().logValue(x: static_cast<Real>(i+1));
71 if (std::fabs(x: calculated-expected)/expected > 1.0e-9)
72 BOOST_ERROR("GammaFunction(" << i << ")\n"
73 << std::setprecision(16) << std::scientific
74 << " calculated: " << calculated << "\n"
75 << " expected: " << expected << "\n"
76 << " rel. error: "
77 << std::fabs(calculated-expected)/expected);
78 }
79}
80
81void FunctionsTest::testGammaValues() {
82
83 BOOST_TEST_MESSAGE("Testing Gamma values...");
84
85 // reference results are calculated with R
86 Real tasks[][3] = {
87 { 0.0001, 9999.422883231624, 1e3},
88 { 1.2, 0.9181687423997607, 1e3},
89 { 7.3, 1271.4236336639089586, 1e3},
90 {-1.1, 9.7148063829028946, 1e3},
91 {-4.001,-41.6040228304425312, 1e3},
92 {-4.999, -8.347576090315059, 1e3},
93 {-19.000001, 8.220610833201313e-12, 1e8},
94 {-19.5, 5.811045977502255e-18, 1e3},
95 {-21.000001, 1.957288098276488e-14, 1e8},
96 {-21.5, 1.318444918321553e-20, 1e6}
97 };
98
99 for (auto& task : tasks) {
100 const Real x = task[0];
101 const Real expected = task[1];
102 const Real calculated = GammaFunction().value(x);
103 const Real tol = task[2] * QL_EPSILON * std::fabs(x: expected);
104
105 if (std::fabs(x: calculated - expected) > tol) {
106 BOOST_ERROR("GammaFunction(" << x << ")\n"
107 << std::setprecision(16) << std::scientific
108 << " calculated: " << calculated << "\n"
109 << " expected: " << expected << "\n"
110 << " rel. error: "
111 << std::fabs(calculated-expected)/expected);
112 }
113 }
114}
115
116void FunctionsTest::testModifiedBesselFunctions() {
117 BOOST_TEST_MESSAGE("Testing modified Bessel function of first and second kind...");
118
119 /* reference values are computed with R and the additional package Bessel
120 * http://cran.r-project.org/web/packages/Bessel
121 */
122
123 Real r[][4] = {
124 {-1.3, 2.0, 1.2079888436539505, 0.1608243636110430},
125 { 1.3, 2.0, 1.2908192151358788, 0.1608243636110430},
126 { 0.001, 2.0, 2.2794705965773794, 0.1138938963603362},
127 { 1.2, 0.5, 0.1768918783499572, 2.1086579232338192},
128 { 2.3, 0.1, 0.00037954958988425198, 572.096866928290183},
129 {-2.3, 1.1, 1.07222017902746969, 1.88152553684107371},
130 {-10.0001, 1.1, 13857.7715614282552, 69288858.9474423379}
131 };
132
133 for (auto& i : r) {
134 const Real nu = i[0];
135 const Real x = i[1];
136 const Real expected_i = i[2];
137 const Real expected_k = i[3];
138 const Real tol_i = 5e4 * QL_EPSILON*std::fabs(x: expected_i);
139 const Real tol_k = 5e4 * QL_EPSILON*std::fabs(x: expected_k);
140
141 const Real calculated_i = modifiedBesselFunction_i(nu, x);
142 const Real calculated_k = modifiedBesselFunction_k(nu, x);
143
144 if (std::fabs(x: expected_i - calculated_i) > tol_i) {
145 BOOST_ERROR("failed to reproduce modified Bessel "
146 << "function of first kind"
147 << "\n order : " << nu
148 << "\n argument : " << x
149 << "\n calculated: " << calculated_i
150 << "\n expected : " << expected_i);
151 }
152 if (std::fabs(x: expected_k - calculated_k) > tol_k) {
153 BOOST_ERROR("failed to reproduce modified Bessel "
154 << "function of second kind"
155 << "\n order : " << nu
156 << "\n argument : " << x
157 << "\n calculated: " << calculated_k
158 << "\n expected : " << expected_k);
159 }
160 }
161
162 Real c[][7] = {
163 {-1.3, 2.0, 0.0, 1.2079888436539505, 0.0,
164 0.1608243636110430, 0.0},
165 { 1.2, 1.5, 0.3, 0.7891550871263575, 0.2721408731632123,
166 0.275126507673411, -0.1316314405663727},
167 { 1.2, -1.5,0.0,-0.6650597524355781, -0.4831941938091643,
168 -0.251112360556051, -2.400130904230102},
169 {-11.2, 1.5, 0.3,12780719.20252659, 16401053.26770633,
170 -34155172.65672453, -43830147.36759921},
171 { 1.2, -1.5,2.0,-0.3869803778520574, 0.9756701796853728,
172 -3.111629716783005, 0.6307859871879062},
173 { 1.2, 0.0, 9.9999,-0.03507838078252647, 0.1079601550451466,
174 -0.05979939995451453, 0.3929814473878203},
175 { 1.2, 0.0, 10.1, -0.02782046891519293, 0.08562259917678558,
176 -0.02035685034691133, 0.3949834389686676},
177 { 1.2, 0.0, 12.1, 0.07092110620741207, -0.2182727210128104,
178 0.3368505862966958, -0.1299038064313366},
179 { 1.2, 0.0, 14.1,-0.03014378676768797, 0.09277303628303372,
180 -0.237531022649052, -0.2351923034581644},
181 { 1.2, 0.0, 16.1,-0.03823210284792657, 0.1176663135266562,
182 -0.1091239402448228, 0.2930535651966139},
183 { 1.2, 0.0, 18.1,0.05626742394733754, -0.173173324361983,
184 0.2941636588154642, -0.02023355577954348},
185 { 1.2, 0.0, 180.1,-0.001230682086826484, 0.003787649998122361,
186 0.02284509628723454, 0.09055419580980778},
187 { 1.2, 0.0, 21.0,-0.04746415965014021, 0.1460796627610969,
188 -0.2693825171336859, -0.04830804448126782},
189 { 1.2, 10.0, 0.0, 2609.784936867044, 0, 1.904394919838336e-05, 0},
190 { 1.2, 14.0, 0.0, 122690.4873454286, 0, 2.902060692576643e-07, 0},
191 { 1.2, 20.0, 10.0, -37452017.91168936, -13917587.22151363,
192 -3.821534367487143e-10, 4.083211255351664e-10},
193 { 1.2, 9.0, 9.0, -621.7335051293694, 618.1455736670332,
194 -4.480795479964915e-05, -3.489034389148745e-08}
195 };
196
197 for (auto& i : c) {
198 const Real nu = i[0];
199 const std::complex<Real> z = std::complex<Real>(i[1], i[2]);
200 const std::complex<Real> expected_i = std::complex<Real>(i[3], i[4]);
201 const std::complex<Real> expected_k = std::complex<Real>(i[5], i[6]);
202
203 const Real tol_i = 5e4*QL_EPSILON*std::abs(z: expected_i);
204 const Real tol_k = 1e6*QL_EPSILON*std::abs(z: expected_k);
205
206 const std::complex<Real> calculated_i=modifiedBesselFunction_i(nu, z);
207 const std::complex<Real> calculated_k=modifiedBesselFunction_k(nu, z);
208
209 if (std::abs(z: expected_i - calculated_i) > tol_i) {
210 BOOST_ERROR("failed to reproduce modified Bessel "
211 << "function of first kind"
212 << "\n order : " << nu
213 << "\n argument : " << z
214 << "\n calculated: " << calculated_i
215 << "\n expected : " << expected_i);
216 }
217 if ( std::abs(z: expected_k) > 1e-4 // do not check small values
218 && std::abs(z: expected_k - calculated_k) > tol_k) {
219 BOOST_ERROR("failed to reproduce modified Bessel "
220 << "function of second kind"
221 << "\n order : " << nu
222 << "\n argument : " << z
223 << "\n diff : " << calculated_k-expected_k
224 << "\n calculated: " << calculated_k
225 << "\n expected : " << expected_k);
226 }
227 }
228}
229
230void FunctionsTest::testWeightedModifiedBesselFunctions() {
231 BOOST_TEST_MESSAGE("Testing weighted modified Bessel functions...");
232 for (Real nu = -5.0; nu <= 5.0; nu += 0.5) {
233 for (Real x = 0.1; x <= 15.0; x += 0.5) {
234 Real calculated_i = modifiedBesselFunction_i_exponentiallyWeighted(nu, x);
235 Real expected_i = modifiedBesselFunction_i(nu, x) * exp(x: -x);
236 Real calculated_k = modifiedBesselFunction_k_exponentiallyWeighted(nu, x);
237 Real expected_k =
238 M_PI_2 * (modifiedBesselFunction_i(nu: -nu, x) - modifiedBesselFunction_i(nu, x)) *
239 exp(x: -x) / std::sin(M_PI * nu);
240 Real tol_i = 1e3 * QL_EPSILON * std::fabs(x: expected_i) * std::max(a: exp(x: x), b: 1.0);
241 Real tol_k = std::max(QL_EPSILON, b: 1e3 * QL_EPSILON * std::fabs(x: expected_k) *
242 std::max(a: exp(x: x), b: 1.0));
243 if (std::abs(x: expected_i - calculated_i) > tol_i) {
244 BOOST_ERROR("failed to verify exponentially weighted"
245 << "modified Bessel function of first kind"
246 << "\n order : " << nu << "\n argument : " << x
247 << "\n calculated : " << calculated_i << "\n expected : "
248 << expected_i << "\n difference : " << (expected_i - calculated_i));
249 }
250 if (std::abs(x: expected_k - calculated_k) > tol_k) {
251 BOOST_ERROR("failed to verify exponentially weighted"
252 << "modified Bessel function of second kind"
253 << "\n order : " << nu << "\n argument : " << x
254 << "\n calculated : " << calculated_k << "\n expected : "
255 << expected_k << "\n difference : " << (expected_k - calculated_k));
256 }
257 }
258 }
259 for (Real nu = -5.0; nu <= 5.0; nu += 0.5) {
260 for (Real x = -5.0; x <= 5.0; x += 0.5) {
261 for (Real y = -5.0; y <= 5.0; y += 0.5) {
262 std::complex<Real> z(x, y);
263 std::complex<Real> calculated_i =
264 modifiedBesselFunction_i_exponentiallyWeighted(nu, z);
265 std::complex<Real> expected_i = modifiedBesselFunction_i(nu, z) * exp(z: -z);
266 std::complex<Real> calculated_k =
267 modifiedBesselFunction_k_exponentiallyWeighted(nu, z);
268 std::complex<Real> expected_k = M_PI_2 *
269 (modifiedBesselFunction_i(nu: -nu, z) * exp(z: -z) -
270 modifiedBesselFunction_i(nu, z) * exp(z: -z)) /
271 std::sin(M_PI * nu);
272 Real tol_i = 1e3 * QL_EPSILON * std::abs(z: calculated_i);
273 Real tol_k = 1e3 * QL_EPSILON * std::abs(z: calculated_k);
274 if (std::abs(z: calculated_i - expected_i) > tol_i) {
275 BOOST_ERROR("failed to verify exponentially weighted"
276 << "modified Bessel function of first kind"
277 << "\n order : " << nu << "\n argument : " << x
278 << "\n calculated : " << calculated_i << "\n expected : "
279 << expected_i << "\n difference : " << (expected_i - calculated_i));
280 }
281 if (std::abs(z: expected_k - calculated_k) > tol_k) {
282 BOOST_ERROR("failed to verify exponentially weighted"
283 << "modified Bessel function of second kind"
284 << "\n order : " << nu << "\n argument : " << x
285 << "\n calculated : " << calculated_k << "\n expected : "
286 << expected_k << "\n difference : " << (expected_k - calculated_k));
287 }
288 }
289 }
290 }
291}
292
293test_suite* FunctionsTest::suite() {
294 auto* suite = BOOST_TEST_SUITE("Factorial tests");
295 suite->add(QUANTLIB_TEST_CASE(&FunctionsTest::testFactorial));
296 suite->add(QUANTLIB_TEST_CASE(&FunctionsTest::testGammaFunction));
297 suite->add(QUANTLIB_TEST_CASE(&FunctionsTest::testGammaValues));
298 suite->add(QUANTLIB_TEST_CASE(
299 &FunctionsTest::testModifiedBesselFunctions));
300 suite->add(QUANTLIB_TEST_CASE(
301 &FunctionsTest::testWeightedModifiedBesselFunctions));
302 return suite;
303}
304

source code of quantlib/test-suite/functions.cpp