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
|
// (C) Copyright Nick Thompson, John Maddock 2023.
// 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)
#include "math_unit_test.hpp"
#include <boost/math/tools/condition_numbers.hpp>
#include <boost/math/tools/estrin.hpp>
#include <numeric>
#include <random>
#if BOOST_MATH_TEST_FLOAT128
#include <boost/multiprecision/float128.hpp>
#endif
using boost::math::tools::evaluate_polynomial_estrin;
using boost::math::tools::summation_condition_number;
template <typename Real, typename Complex> void test_symmetric_coefficients(size_t seed) {
std::mt19937_64 gen(seed);
std::uniform_real_distribution<Real> dis(0, 1);
std::vector<Real> coeffs(0);
for (int i = 0; i < 10; ++i) {
auto p = evaluate_polynomial_estrin(coeffs, Complex(dis(gen)));
if (!CHECK_LE(std::norm(p), Real(0))) {
std::cerr << " If there are zero coefficients, the polynomial should evaluate to zero; seed = " << seed << "\n";
}
}
// Now a single coefficient:
coeffs.resize(1);
coeffs[0] = dis(gen);
for (int i = 0; i < 10; ++i) {
auto p = evaluate_polynomial_estrin(coeffs, Complex(dis(gen)));
if (!CHECK_ULP_CLOSE(p.real(), coeffs[0], 0)) {
std::cerr << " For a polynomial with 1 coefficient, the polynomial should evaluate to c0 for all z; seed = " << seed << "\n";
}
}
}
template <typename Real, typename Complex> void test_polynomial_properties(size_t seed) {
std::mt19937_64 gen(seed);
std::uniform_real_distribution<Real> dis(0, 1);
// We're already testing the n=0 case above:
size_t n = (seed % 128) + 1;
std::vector<Real> coeffs(n);
for (auto &c : coeffs) {
c = dis(gen);
}
// Test evaluation at zero:
auto p0 = evaluate_polynomial_estrin(coeffs, Complex(0));
if (!CHECK_ULP_CLOSE(p0.real(), coeffs[0], 0)) {
std::cerr << " p(0) != c0 with seed " << seed << "\n";
}
// Just given the structure of the calculation, I believe the complex part should be identically zero:
if (!CHECK_ULP_CLOSE(p0.imag(), Real(0), 0)) {
std::cerr << " p(0) != c0 with seed " << seed << "\n";
}
auto p1_computed = evaluate_polynomial_estrin(coeffs, Complex(1));
// The recursive nature of Estrin's method makes it more accurate than a naive sum.
// Use Kahan summation to make sure the expected value is accurate:
auto s = summation_condition_number<Real>(0);
for (auto c : coeffs) {
s += c;
}
if (!CHECK_ULP_CLOSE(s.sum(), p1_computed.real(), coeffs.size())) {
std::cerr << " p(1) != sum(coeffs) with seed " << seed << "\n";
}
if (!CHECK_ULP_CLOSE(Real(0), p1_computed.imag(), coeffs.size())) {
std::cerr << " p(1) != sum(coeffs) with seed " << seed << "\n";
}
}
template <typename Real> void test_std_array_overload(size_t seed) {
std::mt19937_64 gen(seed);
std::uniform_real_distribution<Real> dis(0, 1);
std::array<Real, 12> coeffs;
for (auto &c : coeffs) {
c = dis(gen);
}
// Test evaluation at zero:
auto p0 = evaluate_polynomial_estrin(coeffs, Real(0));
if (!CHECK_ULP_CLOSE(p0, coeffs[0], 0)) {
std::cerr << " p(0) != c0 with seed " << seed << "\n";
}
auto p0_real = evaluate_polynomial_estrin(coeffs, Real(0));
if (!CHECK_ULP_CLOSE(p0_real, coeffs[0], 0)) {
std::cerr << " p(0) != c0 with seed " << seed << "\n";
}
auto p1_computed = evaluate_polynomial_estrin(coeffs, Real(1));
auto s = summation_condition_number<Real>(0);
for (auto c : coeffs) {
s += c;
}
if (!CHECK_ULP_CLOSE(s.sum(), p1_computed, coeffs.size())) {
std::cerr << " p(1) != sum(coeffs) with seed " << seed << "\n";
}
}
int main() {
std::random_device rd;
auto seed = rd();
test_symmetric_coefficients<float, std::complex<float>>(seed);
test_symmetric_coefficients<double, std::complex<double>>(seed);
test_polynomial_properties<float, std::complex<float>>(seed);
test_polynomial_properties<double, std::complex<double>>(seed);
test_std_array_overload<float>(seed);
test_std_array_overload<double>(seed);
#if BOOST_MATH_TEST_FLOAT128
test_std_array_overload<boost::multiprecision::float128>(seed);
#endif
return boost::math::test::report_errors();
}
|
