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