From 483eb2f56657e8e7f419ab1a4fab8dce9ade8609 Mon Sep 17 00:00:00 2001 From: Daniel Baumann Date: Sat, 27 Apr 2024 20:24:20 +0200 Subject: Adding upstream version 14.2.21. Signed-off-by: Daniel Baumann --- src/boost/libs/math/test/test_polynomial.cpp | 614 +++++++++++++++++++++++++++ 1 file changed, 614 insertions(+) create mode 100644 src/boost/libs/math/test/test_polynomial.cpp (limited to 'src/boost/libs/math/test/test_polynomial.cpp') diff --git a/src/boost/libs/math/test/test_polynomial.cpp b/src/boost/libs/math/test/test_polynomial.cpp new file mode 100644 index 00000000..503d0e2a --- /dev/null +++ b/src/boost/libs/math/test/test_polynomial.cpp @@ -0,0 +1,614 @@ +// (C) Copyright Jeremy Murphy 2015. +// 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 +#define BOOST_TEST_MAIN +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#if !defined(TEST1) && !defined(TEST2) && !defined(TEST3) +# define TEST1 +# define TEST2 +# define TEST3 +#endif + +using namespace boost::math; +using boost::integer::gcd; +using namespace boost::math::tools; +using namespace std; +using boost::integer::gcd_detail::Euclid_gcd; +using boost::math::tools::subresultant_gcd; + +template +struct answer +{ + answer(std::pair< polynomial, polynomial > const &x) : + quotient(x.first), remainder(x.second) {} + + polynomial quotient; + polynomial remainder; +}; + +boost::array const d3a = {{10, -6, -4, 3}}; +boost::array const d3b = {{-7, 5, 6, 1}}; + +boost::array const d1a = {{-2, 1}}; +boost::array const d0a = {{6}}; +boost::array const d0a1 = {{0, 6}}; +boost::array const d0a5 = {{0, 0, 0, 0, 0, 6}}; + + +boost::array const d8 = {{-5, 2, 8, -3, -3, 0, 1, 0, 1}}; +boost::array const d8b = {{0, 2, 8, -3, -3, 0, 1, 0, 1}}; + + + +BOOST_AUTO_TEST_CASE(trivial) +{ + /* We have one empty test case here, so that there is always something for Boost.Test to do even if the tests below are #if'ed out */ +} + + +#ifdef TEST1 + +boost::array const d3c = {{10.0/3.0, -2.0, -4.0/3.0, 1.0}}; +boost::array const d2a = {{-2, 2, 3}}; +boost::array const d2b = {{-7, 5, 6}}; +boost::array const d2c = {{31, -21, -22}}; +boost::array const d0b = {{3}}; +boost::array const d6 = {{21, -9, -4, 0, 5, 0, 3}}; +boost::array const d2 = {{-6, 0, 9}}; +boost::array const d5 = {{-9, 0, 3, 0, -15}}; + + +BOOST_AUTO_TEST_CASE( test_construction ) +{ + polynomial const a(d3a.begin(), d3a.end()); + polynomial const b(d3a.begin(), 3); + BOOST_CHECK_EQUAL(a, b); +} + +#ifdef BOOST_MATH_HAS_IS_CONST_ITERABLE + +#include +#include + +BOOST_AUTO_TEST_CASE(test_range_construction) +{ + std::list l{ 1, 2, 3, 4 }; + std::array a{ 3, 4, 5, 6 }; + polynomial p1{ 1, 2, 3, 4 }; + polynomial p2{ 3, 4, 5, 6 }; + + polynomial p3(l); + polynomial p4(a); + + BOOST_CHECK_EQUAL(p1, p3); + BOOST_CHECK_EQUAL(p2, p4); +} +#endif + +#if !defined(BOOST_NO_CXX11_HDR_INITIALIZER_LIST) && !BOOST_WORKAROUND(BOOST_GCC_VERSION, < 40500) +BOOST_AUTO_TEST_CASE( test_initializer_list_construction ) +{ + polynomial a(begin(d3a), end(d3a)); + polynomial b = {10, -6, -4, 3}; + polynomial c{10, -6, -4, 3}; + polynomial d{10, -6, -4, 3, 0, 0}; + BOOST_CHECK_EQUAL(a, b); + BOOST_CHECK_EQUAL(b, c); + BOOST_CHECK_EQUAL(d.degree(), 3u); +} + +BOOST_AUTO_TEST_CASE( test_initializer_list_assignment ) +{ + polynomial a(begin(d3a), end(d3a)); + polynomial b; + b = {10, -6, -4, 3, 0, 0}; + BOOST_CHECK_EQUAL(b.degree(), 3u); + BOOST_CHECK_EQUAL(a, b); +} +#endif + + +BOOST_AUTO_TEST_CASE( test_degree ) +{ + polynomial const zero; + polynomial const a(d3a.begin(), d3a.end()); + BOOST_CHECK_THROW(zero.degree(), std::logic_error); + BOOST_CHECK_EQUAL(a.degree(), 3u); +} + + +BOOST_AUTO_TEST_CASE( test_division_over_field ) +{ + polynomial const a(d3a.begin(), d3a.end()); + polynomial const b(d1a.begin(), d1a.end()); + polynomial const q(d2a.begin(), d2a.end()); + polynomial const r(d0a.begin(), d0a.end()); + polynomial const c(d3b.begin(), d3b.end()); + polynomial const d(d2b.begin(), d2b.end()); + polynomial const e(d2c.begin(), d2c.end()); + polynomial const f(d0b.begin(), d0b.end()); + polynomial const g(d3c.begin(), d3c.end()); + polynomial const zero; + polynomial const one(1.0); + + answer result = quotient_remainder(a, b); + BOOST_CHECK_EQUAL(result.quotient, q); + BOOST_CHECK_EQUAL(result.remainder, r); + BOOST_CHECK_EQUAL(a, q * b + r); // Sanity check. + + result = quotient_remainder(a, c); + BOOST_CHECK_EQUAL(result.quotient, f); + BOOST_CHECK_EQUAL(result.remainder, e); + BOOST_CHECK_EQUAL(a, f * c + e); // Sanity check. + + result = quotient_remainder(a, f); + BOOST_CHECK_EQUAL(result.quotient, g); + BOOST_CHECK_EQUAL(result.remainder, zero); + BOOST_CHECK_EQUAL(a, g * f + zero); // Sanity check. + // Check that division by a regular number gives the same result. + BOOST_CHECK_EQUAL(a / 3.0, g); + BOOST_CHECK_EQUAL(a % 3.0, zero); + + // Sanity checks. + BOOST_CHECK_EQUAL(a / a, one); + BOOST_CHECK_EQUAL(a % a, zero); + // BOOST_CHECK_EQUAL(zero / zero, zero); // TODO +} + +BOOST_AUTO_TEST_CASE( test_division_over_ufd ) +{ + polynomial const zero; + polynomial const one(1); + polynomial const aa(d8.begin(), d8.end()); + polynomial const bb(d6.begin(), d6.end()); + polynomial const q(d2.begin(), d2.end()); + polynomial const r(d5.begin(), d5.end()); + + answer result = quotient_remainder(aa, bb); + BOOST_CHECK_EQUAL(result.quotient, q); + BOOST_CHECK_EQUAL(result.remainder, r); + + // Sanity checks. + BOOST_CHECK_EQUAL(aa / aa, one); + BOOST_CHECK_EQUAL(aa % aa, zero); +} + +#endif + +template +struct FM2GP_Ex_8_3__1 +{ + polynomial x; + polynomial y; + polynomial z; + + FM2GP_Ex_8_3__1() + { + boost::array const x_data = {{105, 278, -88, -56, 16}}; + boost::array const y_data = {{70, 232, -44, -64, 16}}; + boost::array const z_data = {{35, -24, 4}}; + x = polynomial(x_data.begin(), x_data.end()); + y = polynomial(y_data.begin(), y_data.end()); + z = polynomial(z_data.begin(), z_data.end()); + } +}; + +template +struct FM2GP_Ex_8_3__2 +{ + polynomial x; + polynomial y; + polynomial z; + + FM2GP_Ex_8_3__2() + { + boost::array const x_data = {{1, -6, -8, 6, 7}}; + boost::array const y_data = {{1, -5, -2, 15, 11}}; + boost::array const z_data = {{1, 2, 1}}; + x = polynomial(x_data.begin(), x_data.end()); + y = polynomial(y_data.begin(), y_data.end()); + z = polynomial(z_data.begin(), z_data.end()); + } +}; + + +template +struct FM2GP_mixed +{ + polynomial x; + polynomial y; + polynomial z; + + FM2GP_mixed() + { + boost::array const x_data = {{-2.2, -3.3, 0, 1}}; + boost::array const y_data = {{-4.4, 0, 1}}; + boost::array const z_data= {{-2, 1}}; + x = polynomial(x_data.begin(), x_data.end()); + y = polynomial(y_data.begin(), y_data.end()); + z = polynomial(z_data.begin(), z_data.end()); + } +}; + + +template +struct FM2GP_trivial +{ + polynomial x; + polynomial y; + polynomial z; + + FM2GP_trivial() + { + boost::array const x_data = {{-2, -3, 0, 1}}; + boost::array const y_data = {{-4, 0, 1}}; + boost::array const z_data= {{-2, 1}}; + x = polynomial(x_data.begin(), x_data.end()); + y = polynomial(y_data.begin(), y_data.end()); + z = polynomial(z_data.begin(), z_data.end()); + } +}; + +// Sanity checks to make sure I didn't break it. +#ifdef TEST1 +typedef boost::mpl::list integral_test_types; +typedef boost::mpl::list large_integral_test_types; +typedef boost::mpl::list<> mp_integral_test_types; +#elif defined(TEST2) +typedef boost::mpl::list< +#if !BOOST_WORKAROUND(BOOST_MSVC, <= 1500) + boost::multiprecision::cpp_int +#endif +> integral_test_types; +typedef integral_test_types large_integral_test_types; +typedef large_integral_test_types mp_integral_test_types; +#elif defined(TEST3) +typedef boost::mpl::list<> large_integral_test_types; +typedef boost::mpl::list<> integral_test_types; +typedef large_integral_test_types mp_integral_test_types; +#endif + +#ifdef TEST1 +typedef boost::mpl::list non_integral_test_types; +#elif defined(TEST2) +typedef boost::mpl::list< +#if !BOOST_WORKAROUND(BOOST_MSVC, <= 1500) + boost::multiprecision::cpp_rational +#endif +> non_integral_test_types; +#elif defined(TEST3) +typedef boost::mpl::list< +#if !BOOST_WORKAROUND(BOOST_MSVC, <= 1500) + boost::multiprecision::cpp_bin_float_single, boost::multiprecision::cpp_dec_float_50 +#endif +> non_integral_test_types; +#endif + +typedef boost::mpl::joint_view all_test_types; + + +template +void normalize(polynomial &p) +{ + if (leading_coefficient(p) < T(0)) + std::transform(p.data().begin(), p.data().end(), p.data().begin(), std::negate()); +} + +/** + * Note that we do not expect 'pure' gcd algorithms to normalize the result. + * However, the usual public interface function gcd() will do that. + */ + +BOOST_AUTO_TEST_SUITE(test_subresultant_gcd) + +// This test is just to show that gcd>(u, v) is defined (and works) when T is integral and multiprecision. +BOOST_FIXTURE_TEST_CASE_TEMPLATE( gcd_interface, T, mp_integral_test_types, FM2GP_Ex_8_3__1 ) +{ + typedef FM2GP_Ex_8_3__1 fixture_type; + polynomial w; + w = gcd(fixture_type::x, fixture_type::y); + normalize(w); + BOOST_CHECK_EQUAL(w, fixture_type::z); + w = gcd(fixture_type::y, fixture_type::x); + normalize(w); + BOOST_CHECK_EQUAL(w, fixture_type::z); +} + +// This test is just to show that gcd>(u, v) is defined (and works) when T is floating point. +BOOST_FIXTURE_TEST_CASE_TEMPLATE( gcd_float_interface, T, non_integral_test_types, FM2GP_Ex_8_3__1 ) +{ + typedef FM2GP_Ex_8_3__1 fixture_type; + polynomial w; + w = gcd(fixture_type::x, fixture_type::y); + normalize(w); + BOOST_CHECK_EQUAL(w, fixture_type::z); + w = gcd(fixture_type::y, fixture_type::x); + normalize(w); + BOOST_CHECK_EQUAL(w, fixture_type::z); +} + +// The following tests call subresultant_gcd explicitly to remove any ambiguity +// and to permit testing on single-precision integral types. +BOOST_FIXTURE_TEST_CASE_TEMPLATE( Ex_8_3__1, T, large_integral_test_types, FM2GP_Ex_8_3__1 ) +{ + typedef FM2GP_Ex_8_3__1 fixture_type; + polynomial w; + w = subresultant_gcd(fixture_type::x, fixture_type::y); + normalize(w); + BOOST_CHECK_EQUAL(w, fixture_type::z); + w = subresultant_gcd(fixture_type::y, fixture_type::x); + normalize(w); + BOOST_CHECK_EQUAL(w, fixture_type::z); +} + +BOOST_FIXTURE_TEST_CASE_TEMPLATE( Ex_8_3__2, T, large_integral_test_types, FM2GP_Ex_8_3__2 ) +{ + typedef FM2GP_Ex_8_3__2 fixture_type; + polynomial w; + w = subresultant_gcd(fixture_type::x, fixture_type::y); + normalize(w); + BOOST_CHECK_EQUAL(w, fixture_type::z); + w = subresultant_gcd(fixture_type::y, fixture_type::x); + normalize(w); + BOOST_CHECK_EQUAL(w, fixture_type::z); +} + +BOOST_FIXTURE_TEST_CASE_TEMPLATE( trivial_int, T, large_integral_test_types, FM2GP_trivial ) +{ + typedef FM2GP_trivial fixture_type; + polynomial w; + w = subresultant_gcd(fixture_type::x, fixture_type::y); + normalize(w); + BOOST_CHECK_EQUAL(w, fixture_type::z); + w = subresultant_gcd(fixture_type::y, fixture_type::x); + normalize(w); + BOOST_CHECK_EQUAL(w, fixture_type::z); +} + +BOOST_AUTO_TEST_SUITE_END() + + +BOOST_AUTO_TEST_CASE_TEMPLATE( test_addition, T, all_test_types ) +{ + polynomial const a(d3a.begin(), d3a.end()); + polynomial const b(d1a.begin(), d1a.end()); + polynomial const zero; + + polynomial result = a + b; // different degree + boost::array tmp = {{8, -5, -4, 3}}; + polynomial expected(tmp.begin(), tmp.end()); + BOOST_CHECK_EQUAL(result, expected); + BOOST_CHECK_EQUAL(a + zero, a); + BOOST_CHECK_EQUAL(a + b, b + a); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE( test_subtraction, T, all_test_types ) +{ + polynomial const a(d3a.begin(), d3a.end()); + polynomial const zero; + + BOOST_CHECK_EQUAL(a - T(0), a); + BOOST_CHECK_EQUAL(T(0) - a, -a); + BOOST_CHECK_EQUAL(a - zero, a); + BOOST_CHECK_EQUAL(zero - a, -a); + BOOST_CHECK_EQUAL(a - a, zero); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE( test_multiplication, T, all_test_types ) +{ + polynomial const a(d3a.begin(), d3a.end()); + polynomial const b(d1a.begin(), d1a.end()); + polynomial const zero; + boost::array const d3a_sq = {{100, -120, -44, 108, -20, -24, 9}}; + polynomial const a_sq(d3a_sq.begin(), d3a_sq.end()); + + BOOST_CHECK_EQUAL(a * T(0), zero); + BOOST_CHECK_EQUAL(a * zero, zero); + BOOST_CHECK_EQUAL(zero * T(0), zero); + BOOST_CHECK_EQUAL(zero * zero, zero); + BOOST_CHECK_EQUAL(a * b, b * a); + polynomial aa(a); + aa *= aa; + BOOST_CHECK_EQUAL(aa, a_sq); + BOOST_CHECK_EQUAL(aa, a * a); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE( test_arithmetic_relations, T, all_test_types ) +{ + polynomial const a(d8b.begin(), d8b.end()); + polynomial const b(d1a.begin(), d1a.end()); + + BOOST_CHECK_EQUAL(a * T(2), a + a); + BOOST_CHECK_EQUAL(a - b, -b + a); + BOOST_CHECK_EQUAL(a, (a * a) / a); + BOOST_CHECK_EQUAL(a, (a / a) * a); +} + + +BOOST_AUTO_TEST_CASE_TEMPLATE(test_non_integral_arithmetic_relations, T, non_integral_test_types ) +{ + polynomial const a(d8b.begin(), d8b.end()); + polynomial const b(d1a.begin(), d1a.end()); + + BOOST_CHECK_EQUAL(a * T(0.5), a / T(2)); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(test_cont_and_pp, T, integral_test_types) +{ + boost::array, 4> const q={{ + polynomial(d8.begin(), d8.end()), + polynomial(d8b.begin(), d8b.end()), + polynomial(d3a.begin(), d3a.end()), + polynomial(d3b.begin(), d3b.end()) + }}; + for (std::size_t i = 0; i < q.size(); i++) + { + BOOST_CHECK_EQUAL(q[i], content(q[i]) * primitive_part(q[i])); + BOOST_CHECK_EQUAL(primitive_part(q[i]), primitive_part(q[i], content(q[i]))); + } + + polynomial const zero; + BOOST_CHECK_EQUAL(primitive_part(zero), zero); + BOOST_CHECK_EQUAL(content(zero), T(0)); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE( test_self_multiply_assign, T, all_test_types ) +{ + polynomial a(d3a.begin(), d3a.end()); + polynomial const b(a); + boost::array const d3a_sq = {{100, -120, -44, 108, -20, -24, 9}}; + polynomial const asq(d3a_sq.begin(), d3a_sq.end()); + + a *= a; + + BOOST_CHECK_EQUAL(a, b*b); + BOOST_CHECK_EQUAL(a, asq); + + a *= a; + + BOOST_CHECK_EQUAL(a, b*b*b*b); +} + + +BOOST_AUTO_TEST_CASE_TEMPLATE(test_right_shift, T, all_test_types ) +{ + polynomial a(d8b.begin(), d8b.end()); + polynomial const aa(a); + polynomial const b(d8b.begin() + 1, d8b.end()); + polynomial const c(d8b.begin() + 5, d8b.end()); + a >>= 0u; + BOOST_CHECK_EQUAL(a, aa); + a >>= 1u; + BOOST_CHECK_EQUAL(a, b); + a = a >> 4u; + BOOST_CHECK_EQUAL(a, c); +} + + +BOOST_AUTO_TEST_CASE_TEMPLATE(test_left_shift, T, all_test_types ) +{ + polynomial a(d0a.begin(), d0a.end()); + polynomial const aa(a); + polynomial const b(d0a1.begin(), d0a1.end()); + polynomial const c(d0a5.begin(), d0a5.end()); + a <<= 0u; + BOOST_CHECK_EQUAL(a, aa); + a <<= 1u; + BOOST_CHECK_EQUAL(a, b); + a = a << 4u; + BOOST_CHECK_EQUAL(a, c); + polynomial zero; + // Multiplying zero by x should still be zero. + zero <<= 1u; + BOOST_CHECK_EQUAL(zero, zero_element(multiplies< polynomial >())); +} + + +BOOST_AUTO_TEST_CASE_TEMPLATE(test_odd_even, T, all_test_types) +{ + polynomial const zero; + BOOST_CHECK_EQUAL(odd(zero), false); + BOOST_CHECK_EQUAL(even(zero), true); + polynomial const a(d0a.begin(), d0a.end()); + BOOST_CHECK_EQUAL(odd(a), true); + BOOST_CHECK_EQUAL(even(a), false); + polynomial const b(d0a1.begin(), d0a1.end()); + BOOST_CHECK_EQUAL(odd(b), false); + BOOST_CHECK_EQUAL(even(b), true); +} + +// NOTE: Slightly unexpected: this unit test passes even when T = char. +BOOST_AUTO_TEST_CASE_TEMPLATE( test_pow, T, all_test_types ) +{ + if (std::numeric_limits::digits < 32) + return; // Invokes undefined behaviour + polynomial a(d3a.begin(), d3a.end()); + polynomial const one(T(1)); + boost::array const d3a_sqr = {{100, -120, -44, 108, -20, -24, 9}}; + boost::array const d3a_cub = + {{1000, -1800, -120, 2124, -1032, -684, 638, -18, -108, 27}}; + polynomial const asqr(d3a_sqr.begin(), d3a_sqr.end()); + polynomial const acub(d3a_cub.begin(), d3a_cub.end()); + + BOOST_CHECK_EQUAL(pow(a, 0), one); + BOOST_CHECK_EQUAL(pow(a, 1), a); + BOOST_CHECK_EQUAL(pow(a, 2), asqr); + BOOST_CHECK_EQUAL(pow(a, 3), acub); + BOOST_CHECK_EQUAL(pow(a, 4), pow(asqr, 2)); + BOOST_CHECK_EQUAL(pow(a, 5), asqr * acub); + BOOST_CHECK_EQUAL(pow(a, 6), pow(acub, 2)); + BOOST_CHECK_EQUAL(pow(a, 7), acub * acub * a); + + BOOST_CHECK_THROW(pow(a, -1), std::domain_error); + BOOST_CHECK_EQUAL(pow(one, 137), one); +} + + +BOOST_AUTO_TEST_CASE_TEMPLATE(test_bool, T, all_test_types) +{ + polynomial const zero; + polynomial const a(d0a.begin(), d0a.end()); + BOOST_CHECK_EQUAL(bool(zero), false); + BOOST_CHECK_EQUAL(bool(a), true); +} + + +BOOST_AUTO_TEST_CASE_TEMPLATE(test_set_zero, T, all_test_types) +{ + polynomial const zero; + polynomial a(d0a.begin(), d0a.end()); + a.set_zero(); + BOOST_CHECK_EQUAL(a, zero); + a.set_zero(); // Ensure that setting zero to zero is a no-op. + BOOST_CHECK_EQUAL(a, zero); +} + + +BOOST_AUTO_TEST_CASE_TEMPLATE(test_leading_coefficient, T, all_test_types) +{ + polynomial const zero; + BOOST_CHECK_EQUAL(leading_coefficient(zero), T(0)); + polynomial a(d0a.begin(), d0a.end()); + BOOST_CHECK_EQUAL(leading_coefficient(a), T(d0a.back())); +} + +#if !defined(BOOST_NO_CXX11_RVALUE_REFERENCES) && !defined(BOOST_NO_CXX11_UNIFIED_INITIALIZATION_SYNTAX) +BOOST_AUTO_TEST_CASE_TEMPLATE(test_prime, T, all_test_types) +{ + std::vector d{1,1,1,1,1}; + polynomial p(std::move(d)); + polynomial q = p.prime(); + BOOST_CHECK_EQUAL(q(0), T(1)); + + for (size_t i = 0; i < q.size(); ++i) + { + BOOST_CHECK_EQUAL(q[i], i+1); + } + + polynomial P = p.integrate(); + BOOST_CHECK_EQUAL(P(0), T(0)); + for (size_t i = 1; i < P.size(); ++i) + { + BOOST_CHECK_EQUAL(P[i], 1/static_cast(i)); + } + + polynomial empty; + q = empty.prime(); + BOOST_CHECK_EQUAL(q.size(), 0); + +} +#endif -- cgit v1.2.3