1 // (C) Copyright Nick Thompson, 2019
2 // Use, modification and distribution are subject to the
3 // Boost Software License, Version 1.0. (See accompanying file
4 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
5
6 #define BOOST_TEST_MODULE condition_number_test
7
8 #include <cmath>
9 #include <limits>
10 #include <iostream>
11 #include <boost/math/constants/constants.hpp>
12 #include <boost/math/special_functions/lambert_w.hpp>
13 #include <boost/test/included/unit_test.hpp>
14 #include <boost/multiprecision/cpp_bin_float.hpp>
15 #include <boost/math/tools/condition_numbers.hpp>
16
17 using std::abs;
18 using boost::math::constants::half;
19 using boost::math::constants::ln_two;
20 using boost::multiprecision::cpp_bin_float_50;
21 using boost::math::tools::summation_condition_number;
22 using boost::math::tools::evaluation_condition_number;
23
24 template<class Real>
test_summation_condition_number()25 void test_summation_condition_number()
26 {
27 Real tol = 1000*std::numeric_limits<float>::epsilon();
28 auto cond = summation_condition_number<Real>();
29 // I've checked that the condition number increases with max_n,
30 // and that the computed sum gets more accurate with increasing max_n.
31 // But the CI system would die with more terms.
32 Real max_n = 10000;
33 for (Real n = 1; n < max_n; n += 2)
34 {
35 cond += 1/n;
36 cond -= 1/(n+1);
37 }
38
39 BOOST_CHECK_CLOSE_FRACTION(cond.sum(), ln_two<Real>(), tol);
40 BOOST_TEST(cond() > 14);
41 }
42
43 template<class Real>
test_exponential_sum()44 void test_exponential_sum()
45 {
46 using std::exp;
47 using std::abs;
48 Real eps = std::numeric_limits<float>::epsilon();
49 for (Real x = -20; x <= -1; x += 0.5)
50 {
51 auto cond = summation_condition_number<Real>(1);
52 size_t n = 1;
53 Real term = x;
54 while(n++ < 1000)
55 {
56 cond += term;
57 term *= (x/n);
58 }
59 BOOST_CHECK_CLOSE_FRACTION(exp(x), cond.sum(), eps*cond());
60 BOOST_CHECK_CLOSE_FRACTION(exp(2*abs(x)), cond(), eps*cond());
61 }
62 }
63
64
65
66 template<class Real>
test_evaluation_condition_number()67 void test_evaluation_condition_number()
68 {
69 using std::abs;
70 using std::log;
71 using std::sqrt;
72 using std::exp;
73 using std::sin;
74 using std::tan;
75 Real tol = sqrt(std::numeric_limits<Real>::epsilon());
76
77 auto f1 = [](auto x) { return log(x); };
78 for (Real x = 1.125; x < 8; x += 0.125)
79 {
80 Real cond = evaluation_condition_number(f1, x);
81 BOOST_CHECK_CLOSE_FRACTION(cond, 1/log(x), tol);
82 }
83
84 auto f2 = [](auto x) { return exp(x); };
85 for (Real x = 1.125; x < 8; x += 0.125)
86 {
87 Real cond = evaluation_condition_number(f2, x);
88 BOOST_CHECK_CLOSE_FRACTION(cond, x, tol);
89 }
90
91 auto f3 = [](auto x) { return sin(x); };
92 for (Real x = 1.125; x < 8; x += 0.125)
93 {
94 Real cond = evaluation_condition_number(f3, x);
95 BOOST_CHECK_CLOSE_FRACTION(cond, abs(x/tan(x)), tol);
96 }
97
98 // Test a function which right differentiable:
99 using boost::math::constants::e;
100 auto f4 = [](Real x) { return boost::math::lambert_w0(x); };
101 Real cond = evaluation_condition_number(f4, -1/e<Real>());
102 if (std::is_same_v<Real, float>)
103 {
104 BOOST_CHECK_GE(cond, 30);
105 }
106 else
107 {
108 BOOST_CHECK_GE(cond, 4900);
109 }
110 }
111
112
BOOST_AUTO_TEST_CASE(numerical_differentiation_test)113 BOOST_AUTO_TEST_CASE(numerical_differentiation_test)
114 {
115 test_summation_condition_number<float>();
116 test_summation_condition_number<cpp_bin_float_50>();
117 test_evaluation_condition_number<float>();
118 test_evaluation_condition_number<double>();
119 test_evaluation_condition_number<long double>();
120 test_evaluation_condition_number<cpp_bin_float_50>();
121 test_exponential_sum<double>();
122 }
123