[ceph.git] / ceph / src / boost / libs / math / test / test_polynomial.cpp

//  (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 <boost/config.hpp>
#define BOOST_TEST_MAIN
#include <boost/array.hpp>
#include <boost/math/tools/polynomial.hpp>
#include <boost/integer/common_factor_rt.hpp>
#include <boost/mpl/list.hpp>
#include <boost/mpl/joint_view.hpp>
#include <boost/test/test_case_template.hpp>
#include <boost/test/unit_test.hpp>
#include <boost/multiprecision/cpp_int.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>
#include <utility>

#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 <typename T>
struct answer
{
    answer(std::pair< polynomial<T>, polynomial<T> > const &x) :
    quotient(x.first), remainder(x.second) {}

    polynomial<T> quotient;
    polynomial<T> remainder;
};

boost::array<double, 4> const d3a = {{10, -6, -4, 3}};
boost::array<double, 4> const d3b = {{-7, 5, 6, 1}};

boost::array<double, 2> const d1a = {{-2, 1}};
boost::array<double, 1> const d0a = {{6}};
boost::array<double, 2> const d0a1 = {{0, 6}};
boost::array<double, 6> const d0a5 = {{0, 0, 0, 0, 0, 6}};


boost::array<int, 9> const d8 = {{-5, 2, 8, -3, -3, 0, 1, 0, 1}};
boost::array<int, 9> 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<double, 4> const d3c = {{10.0/3.0, -2.0, -4.0/3.0, 1.0}};
boost::array<double, 3> const d2a = {{-2, 2, 3}};
boost::array<double, 3> const d2b = {{-7, 5, 6}};
boost::array<double, 3> const d2c = {{31, -21, -22}};
boost::array<double, 1> const d0b = {{3}};
boost::array<int, 7> const d6 = {{21, -9, -4, 0, 5, 0, 3}};
boost::array<int, 3> const d2 = {{-6, 0, 9}};
boost::array<int, 6> const d5 = {{-9, 0, 3, 0, -15}};


BOOST_AUTO_TEST_CASE( test_construction )
{
    polynomial<double> const a(d3a.begin(), d3a.end());
    polynomial<double> const b(d3a.begin(), 3);
    BOOST_CHECK_EQUAL(a, b);
}

#ifdef BOOST_MATH_HAS_IS_CONST_ITERABLE

#include <list>
#include <array>

BOOST_AUTO_TEST_CASE(test_range_construction)
{
   std::list<double> l{ 1, 2, 3, 4 };
   std::array<double, 4> a{ 3, 4, 5, 6 };
   polynomial<double> p1{ 1, 2, 3, 4 };
   polynomial<double> p2{ 3, 4, 5, 6 };

   polynomial<double> p3(l);
   polynomial<double> 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<double> a(begin(d3a), end(d3a));
    polynomial<double> b = {10, -6, -4, 3};
    polynomial<double> c{10, -6, -4, 3};
    polynomial<double> 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<double> a(begin(d3a), end(d3a));
    polynomial<double> 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<double> const zero;
    polynomial<double> 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<double> const a(d3a.begin(), d3a.end());
    polynomial<double> const b(d1a.begin(), d1a.end());
    polynomial<double> const q(d2a.begin(), d2a.end());
    polynomial<double> const r(d0a.begin(), d0a.end());
    polynomial<double> const c(d3b.begin(), d3b.end());
    polynomial<double> const d(d2b.begin(), d2b.end());
    polynomial<double> const e(d2c.begin(), d2c.end());
    polynomial<double> const f(d0b.begin(), d0b.end());
    polynomial<double> const g(d3c.begin(), d3c.end());
    polynomial<double> const zero;
    polynomial<double> const one(1.0);

    answer<double> 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<int> const zero;
    polynomial<int> const one(1);
    polynomial<int> const aa(d8.begin(), d8.end());
    polynomial<int> const bb(d6.begin(), d6.end());
    polynomial<int> const q(d2.begin(), d2.end());
    polynomial<int> const r(d5.begin(), d5.end());

    answer<int> 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 <typename T>
struct FM2GP_Ex_8_3__1
{
    polynomial<T> x;
    polynomial<T> y;
    polynomial<T> z;

    FM2GP_Ex_8_3__1()
    {
        boost::array<T, 5> const x_data = {{105, 278, -88, -56, 16}};
        boost::array<T, 5> const y_data = {{70, 232, -44, -64, 16}};
        boost::array<T, 3> const z_data = {{35, -24, 4}};
        x = polynomial<T>(x_data.begin(), x_data.end());
        y = polynomial<T>(y_data.begin(), y_data.end());
        z = polynomial<T>(z_data.begin(), z_data.end());
    }
};

template <typename T>
struct FM2GP_Ex_8_3__2
{
    polynomial<T> x;
    polynomial<T> y;
    polynomial<T> z;

    FM2GP_Ex_8_3__2()
    {
        boost::array<T, 5> const x_data = {{1, -6, -8, 6, 7}};
        boost::array<T, 5> const y_data = {{1, -5, -2, 15, 11}};
        boost::array<T, 3> const z_data = {{1, 2, 1}};
        x = polynomial<T>(x_data.begin(), x_data.end());
        y = polynomial<T>(y_data.begin(), y_data.end());
        z = polynomial<T>(z_data.begin(), z_data.end());
    }
};


template <typename T>
struct FM2GP_mixed
{
    polynomial<T> x;
    polynomial<T> y;
    polynomial<T> z;

    FM2GP_mixed()
    {
        boost::array<T, 4> const x_data = {{-2.2, -3.3, 0, 1}};
        boost::array<T, 3> const y_data = {{-4.4, 0, 1}};
        boost::array<T, 2> const z_data= {{-2, 1}};
        x = polynomial<T>(x_data.begin(), x_data.end());
        y = polynomial<T>(y_data.begin(), y_data.end());
        z = polynomial<T>(z_data.begin(), z_data.end());
    }
};


template <typename T>
struct FM2GP_trivial
{
    polynomial<T> x;
    polynomial<T> y;
    polynomial<T> z;

    FM2GP_trivial()
    {
        boost::array<T, 4> const x_data = {{-2, -3, 0, 1}};
        boost::array<T, 3> const y_data = {{-4, 0, 1}};
        boost::array<T, 2> const z_data= {{-2, 1}};
        x = polynomial<T>(x_data.begin(), x_data.end());
        y = polynomial<T>(y_data.begin(), y_data.end());
        z = polynomial<T>(z_data.begin(), z_data.end());
    }
};

// Sanity checks to make sure I didn't break it.
#ifdef TEST1
typedef boost::mpl::list<char, short, int, long> integral_test_types;
typedef boost::mpl::list<int, long> 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<double, long double> 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<integral_test_types, non_integral_test_types> all_test_types;


template <typename T>
void normalize(polynomial<T> &p)
{
    if (leading_coefficient(p) < T(0))
        std::transform(p.data().begin(), p.data().end(), p.data().begin(), std::negate<T>());
}

/**
 * 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<polynomial<T>>(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<T> )
{
    typedef FM2GP_Ex_8_3__1<T> fixture_type;
    polynomial<T> 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<polynomial<T>>(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<T> )
{
    typedef FM2GP_Ex_8_3__1<T> fixture_type;
    polynomial<T> 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<T> )
{
    typedef FM2GP_Ex_8_3__1<T> fixture_type;
    polynomial<T> 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<T> )
{
    typedef FM2GP_Ex_8_3__2<T> fixture_type;
    polynomial<T> 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<T> )
{
    typedef FM2GP_trivial<T> fixture_type;
    polynomial<T> 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<T> const a(d3a.begin(), d3a.end());
    polynomial<T> const b(d1a.begin(), d1a.end());
    polynomial<T> const zero;

    polynomial<T> result = a + b; // different degree
    boost::array<T, 4> tmp = {{8, -5, -4, 3}};
    polynomial<T> 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<T> const a(d3a.begin(), d3a.end());
    polynomial<T> 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<T> const a(d3a.begin(), d3a.end());
    polynomial<T> const b(d1a.begin(), d1a.end());
    polynomial<T> const zero;
    boost::array<T, 7> const d3a_sq = {{100, -120, -44, 108, -20, -24, 9}};
    polynomial<T> 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<T> 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<T> const a(d8b.begin(), d8b.end());
    polynomial<T> 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<T> const a(d8b.begin(), d8b.end());
    polynomial<T> 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<polynomial<T>, 4> const q={{
        polynomial<T>(d8.begin(), d8.end()),
        polynomial<T>(d8b.begin(), d8b.end()),
        polynomial<T>(d3a.begin(), d3a.end()),
        polynomial<T>(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<T> 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<T> a(d3a.begin(), d3a.end());
    polynomial<T> const b(a);
    boost::array<double, 7> const d3a_sq = {{100, -120, -44, 108, -20, -24, 9}};
    polynomial<T> 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<T> a(d8b.begin(), d8b.end());
    polynomial<T> const aa(a);
    polynomial<T> const b(d8b.begin() + 1, d8b.end());
    polynomial<T> 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<T> a(d0a.begin(), d0a.end());
    polynomial<T> const aa(a);
    polynomial<T> const b(d0a1.begin(), d0a1.end());
    polynomial<T> 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<T> zero;
    // Multiplying zero by x should still be zero.
    zero <<= 1u;
    BOOST_CHECK_EQUAL(zero, zero_element(multiplies< polynomial<T> >()));
}


BOOST_AUTO_TEST_CASE_TEMPLATE(test_odd_even, T, all_test_types)
{
    polynomial<T> const zero;
    BOOST_CHECK_EQUAL(odd(zero), false);
    BOOST_CHECK_EQUAL(even(zero), true);
    polynomial<T> const a(d0a.begin(), d0a.end());
    BOOST_CHECK_EQUAL(odd(a), true);
    BOOST_CHECK_EQUAL(even(a), false);
    polynomial<T> 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<T>::digits < 32)
      return;   // Invokes undefined behaviour
    polynomial<T> a(d3a.begin(), d3a.end());
    polynomial<T> const one(T(1));
    boost::array<double, 7> const d3a_sqr = {{100, -120, -44, 108, -20, -24, 9}};
    boost::array<double, 10> const d3a_cub =
        {{1000, -1800, -120, 2124, -1032, -684, 638, -18, -108, 27}};
    polynomial<T> const asqr(d3a_sqr.begin(), d3a_sqr.end());
    polynomial<T> 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<T> const zero;
    polynomial<T> 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<T> const zero;
    polynomial<T> 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<T> const zero;
    BOOST_CHECK_EQUAL(leading_coefficient(zero), T(0));
    polynomial<T> 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<T> d{1,1,1,1,1};
    polynomial<T> p(std::move(d));
    polynomial<T> 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<T> 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<T>(i));
    }

    polynomial<T> empty;
    q = empty.prime();
    BOOST_CHECK_EQUAL(q.size(), 0);

}
#endif
Commit	Line	Data
7c673cae FG	1	// (C) Copyright Jeremy Murphy 2015.
	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	#include <boost/config.hpp>
	7	#define BOOST_TEST_MAIN
	8	#include <boost/array.hpp>
	9	#include <boost/math/tools/polynomial.hpp>
11fdf7f2	10	#include <boost/integer/common_factor_rt.hpp>
7c673cae FG	11	#include <boost/mpl/list.hpp>
	12	#include <boost/mpl/joint_view.hpp>
	13	#include <boost/test/test_case_template.hpp>
	14	#include <boost/test/unit_test.hpp>
	15	#include <boost/multiprecision/cpp_int.hpp>
	16	#include <boost/multiprecision/cpp_bin_float.hpp>
	17	#include <boost/multiprecision/cpp_dec_float.hpp>
	18	#include <utility>
	19
b32b8144 FG	20	#if !defined(TEST1) && !defined(TEST2) && !defined(TEST3)
	21	# define TEST1
	22	# define TEST2
	23	# define TEST3
	24	#endif
	25
	26	using namespace boost::math;
11fdf7f2	27	using boost::integer::gcd;
7c673cae FG	28	using namespace boost::math::tools;
7c673cae FG	29	using namespace std;
b32b8144 FG	30	using boost::integer::gcd_detail::Euclid_gcd;
b32b8144 FG	31	using boost::math::tools::subresultant_gcd;
7c673cae FG	32
	33	template <typename T>
	34	struct answer
	35	{
	36	answer(std::pair< polynomial<T>, polynomial<T> > const &x) :
	37	quotient(x.first), remainder(x.second) {}
11fdf7f2	38
7c673cae FG	39	polynomial<T> quotient;
	40	polynomial<T> remainder;
	41	};
	42
	43	boost::array<double, 4> const d3a = {{10, -6, -4, 3}};
	44	boost::array<double, 4> const d3b = {{-7, 5, 6, 1}};
92f5a8d4	45
7c673cae	46	boost::array<double, 2> const d1a = {{-2, 1}};
7c673cae FG	47	boost::array<double, 1> const d0a = {{6}};
	48	boost::array<double, 2> const d0a1 = {{0, 6}};
	49	boost::array<double, 6> const d0a5 = {{0, 0, 0, 0, 0, 6}};
92f5a8d4	50
7c673cae FG	51
	52	boost::array<int, 9> const d8 = {{-5, 2, 8, -3, -3, 0, 1, 0, 1}};
	53	boost::array<int, 9> const d8b = {{0, 2, 8, -3, -3, 0, 1, 0, 1}};
92f5a8d4 TL	54
92f5a8d4 TL	55
7c673cae	56
b32b8144 FG	57	BOOST_AUTO_TEST_CASE(trivial)
	58	{
	59	/* 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 */
	60	}
	61
	62
	63	#ifdef TEST1
	64
92f5a8d4 TL	65	boost::array<double, 4> const d3c = {{10.0/3.0, -2.0, -4.0/3.0, 1.0}};
	66	boost::array<double, 3> const d2a = {{-2, 2, 3}};
	67	boost::array<double, 3> const d2b = {{-7, 5, 6}};
	68	boost::array<double, 3> const d2c = {{31, -21, -22}};
	69	boost::array<double, 1> const d0b = {{3}};
	70	boost::array<int, 7> const d6 = {{21, -9, -4, 0, 5, 0, 3}};
	71	boost::array<int, 3> const d2 = {{-6, 0, 9}};
	72	boost::array<int, 6> const d5 = {{-9, 0, 3, 0, -15}};
	73
7c673cae FG	74
	75	BOOST_AUTO_TEST_CASE( test_construction )
	76	{
	77	polynomial<double> const a(d3a.begin(), d3a.end());
	78	polynomial<double> const b(d3a.begin(), 3);
	79	BOOST_CHECK_EQUAL(a, b);
	80	}
	81
92f5a8d4 TL	82	#ifdef BOOST_MATH_HAS_IS_CONST_ITERABLE
	83
	84	#include <list>
	85	#include <array>
	86
	87	BOOST_AUTO_TEST_CASE(test_range_construction)
	88	{
	89	std::list<double> l{ 1, 2, 3, 4 };
	90	std::array<double, 4> a{ 3, 4, 5, 6 };
	91	polynomial<double> p1{ 1, 2, 3, 4 };
	92	polynomial<double> p2{ 3, 4, 5, 6 };
	93
	94	polynomial<double> p3(l);
	95	polynomial<double> p4(a);
	96
	97	BOOST_CHECK_EQUAL(p1, p3);
	98	BOOST_CHECK_EQUAL(p2, p4);
	99	}
	100	#endif
7c673cae FG	101
	102	#if !defined(BOOST_NO_CXX11_HDR_INITIALIZER_LIST) && !BOOST_WORKAROUND(BOOST_GCC_VERSION, < 40500)
	103	BOOST_AUTO_TEST_CASE( test_initializer_list_construction )
	104	{
	105	polynomial<double> a(begin(d3a), end(d3a));
	106	polynomial<double> b = {10, -6, -4, 3};
92f5a8d4 TL	107	polynomial<double> c{10, -6, -4, 3};
92f5a8d4 TL	108	polynomial<double> d{10, -6, -4, 3, 0, 0};
7c673cae FG	109	BOOST_CHECK_EQUAL(a, b);
	110	BOOST_CHECK_EQUAL(b, c);
	111	BOOST_CHECK_EQUAL(d.degree(), 3u);
	112	}
	113
	114	BOOST_AUTO_TEST_CASE( test_initializer_list_assignment )
	115	{
	116	polynomial<double> a(begin(d3a), end(d3a));
	117	polynomial<double> b;
	118	b = {10, -6, -4, 3, 0, 0};
	119	BOOST_CHECK_EQUAL(b.degree(), 3u);
	120	BOOST_CHECK_EQUAL(a, b);
	121	}
	122	#endif
	123
	124
	125	BOOST_AUTO_TEST_CASE( test_degree )
	126	{
	127	polynomial<double> const zero;
	128	polynomial<double> const a(d3a.begin(), d3a.end());
	129	BOOST_CHECK_THROW(zero.degree(), std::logic_error);
	130	BOOST_CHECK_EQUAL(a.degree(), 3u);
	131	}
	132
	133
	134	BOOST_AUTO_TEST_CASE( test_division_over_field )
	135	{
	136	polynomial<double> const a(d3a.begin(), d3a.end());
	137	polynomial<double> const b(d1a.begin(), d1a.end());
	138	polynomial<double> const q(d2a.begin(), d2a.end());
	139	polynomial<double> const r(d0a.begin(), d0a.end());
	140	polynomial<double> const c(d3b.begin(), d3b.end());
	141	polynomial<double> const d(d2b.begin(), d2b.end());
	142	polynomial<double> const e(d2c.begin(), d2c.end());
	143	polynomial<double> const f(d0b.begin(), d0b.end());
	144	polynomial<double> const g(d3c.begin(), d3c.end());
	145	polynomial<double> const zero;
	146	polynomial<double> const one(1.0);
	147
	148	answer<double> result = quotient_remainder(a, b);
	149	BOOST_CHECK_EQUAL(result.quotient, q);
	150	BOOST_CHECK_EQUAL(result.remainder, r);
	151	BOOST_CHECK_EQUAL(a, q * b + r); // Sanity check.
11fdf7f2	152
7c673cae FG	153	result = quotient_remainder(a, c);
	154	BOOST_CHECK_EQUAL(result.quotient, f);
	155	BOOST_CHECK_EQUAL(result.remainder, e);
	156	BOOST_CHECK_EQUAL(a, f * c + e); // Sanity check.
11fdf7f2	157
7c673cae FG	158	result = quotient_remainder(a, f);
	159	BOOST_CHECK_EQUAL(result.quotient, g);
	160	BOOST_CHECK_EQUAL(result.remainder, zero);
	161	BOOST_CHECK_EQUAL(a, g * f + zero); // Sanity check.
	162	// Check that division by a regular number gives the same result.
	163	BOOST_CHECK_EQUAL(a / 3.0, g);
	164	BOOST_CHECK_EQUAL(a % 3.0, zero);
	165
	166	// Sanity checks.
	167	BOOST_CHECK_EQUAL(a / a, one);
	168	BOOST_CHECK_EQUAL(a % a, zero);
	169	// BOOST_CHECK_EQUAL(zero / zero, zero); // TODO
	170	}
	171
	172	BOOST_AUTO_TEST_CASE( test_division_over_ufd )
	173	{
	174	polynomial<int> const zero;
	175	polynomial<int> const one(1);
	176	polynomial<int> const aa(d8.begin(), d8.end());
	177	polynomial<int> const bb(d6.begin(), d6.end());
	178	polynomial<int> const q(d2.begin(), d2.end());
	179	polynomial<int> const r(d5.begin(), d5.end());
11fdf7f2	180
7c673cae FG	181	answer<int> result = quotient_remainder(aa, bb);
	182	BOOST_CHECK_EQUAL(result.quotient, q);
	183	BOOST_CHECK_EQUAL(result.remainder, r);
	184
	185	// Sanity checks.
	186	BOOST_CHECK_EQUAL(aa / aa, one);
	187	BOOST_CHECK_EQUAL(aa % aa, zero);
	188	}
	189
b32b8144 FG	190	#endif
	191
	192	template <typename T>
	193	struct FM2GP_Ex_8_3__1
	194	{
	195	polynomial<T> x;
	196	polynomial<T> y;
	197	polynomial<T> z;
11fdf7f2	198
b32b8144 FG	199	FM2GP_Ex_8_3__1()
	200	{
	201	boost::array<T, 5> const x_data = {{105, 278, -88, -56, 16}};
	202	boost::array<T, 5> const y_data = {{70, 232, -44, -64, 16}};
	203	boost::array<T, 3> const z_data = {{35, -24, 4}};
	204	x = polynomial<T>(x_data.begin(), x_data.end());
	205	y = polynomial<T>(y_data.begin(), y_data.end());
	206	z = polynomial<T>(z_data.begin(), z_data.end());
	207	}
	208	};
	209
	210	template <typename T>
	211	struct FM2GP_Ex_8_3__2
	212	{
	213	polynomial<T> x;
	214	polynomial<T> y;
	215	polynomial<T> z;
11fdf7f2	216
b32b8144 FG	217	FM2GP_Ex_8_3__2()
	218	{
	219	boost::array<T, 5> const x_data = {{1, -6, -8, 6, 7}};
	220	boost::array<T, 5> const y_data = {{1, -5, -2, 15, 11}};
	221	boost::array<T, 3> const z_data = {{1, 2, 1}};
	222	x = polynomial<T>(x_data.begin(), x_data.end());
	223	y = polynomial<T>(y_data.begin(), y_data.end());
	224	z = polynomial<T>(z_data.begin(), z_data.end());
	225	}
	226	};
7c673cae	227
b32b8144 FG	228
	229	template <typename T>
	230	struct FM2GP_mixed
	231	{
	232	polynomial<T> x;
	233	polynomial<T> y;
	234	polynomial<T> z;
11fdf7f2	235
b32b8144 FG	236	FM2GP_mixed()
	237	{
	238	boost::array<T, 4> const x_data = {{-2.2, -3.3, 0, 1}};
	239	boost::array<T, 3> const y_data = {{-4.4, 0, 1}};
	240	boost::array<T, 2> const z_data= {{-2, 1}};
	241	x = polynomial<T>(x_data.begin(), x_data.end());
	242	y = polynomial<T>(y_data.begin(), y_data.end());
	243	z = polynomial<T>(z_data.begin(), z_data.end());
	244	}
	245	};
	246
	247
	248	template <typename T>
	249	struct FM2GP_trivial
	250	{
	251	polynomial<T> x;
	252	polynomial<T> y;
	253	polynomial<T> z;
11fdf7f2	254
b32b8144 FG	255	FM2GP_trivial()
	256	{
	257	boost::array<T, 4> const x_data = {{-2, -3, 0, 1}};
	258	boost::array<T, 3> const y_data = {{-4, 0, 1}};
	259	boost::array<T, 2> const z_data= {{-2, 1}};
	260	x = polynomial<T>(x_data.begin(), x_data.end());
	261	y = polynomial<T>(y_data.begin(), y_data.end());
	262	z = polynomial<T>(z_data.begin(), z_data.end());
	263	}
	264	};
7c673cae FG	265
7c673cae FG	266	// Sanity checks to make sure I didn't break it.
b32b8144 FG	267	#ifdef TEST1
	268	typedef boost::mpl::list<char, short, int, long> integral_test_types;
	269	typedef boost::mpl::list<int, long> large_integral_test_types;
	270	typedef boost::mpl::list<> mp_integral_test_types;
	271	#elif defined(TEST2)
	272	typedef boost::mpl::list<
7c673cae	273	#if !BOOST_WORKAROUND(BOOST_MSVC, <= 1500)
b32b8144	274	boost::multiprecision::cpp_int
7c673cae FG	275	#endif
7c673cae FG	276	> integral_test_types;
b32b8144 FG	277	typedef integral_test_types large_integral_test_types;
	278	typedef large_integral_test_types mp_integral_test_types;
	279	#elif defined(TEST3)
	280	typedef boost::mpl::list<> large_integral_test_types;
	281	typedef boost::mpl::list<> integral_test_types;
	282	typedef large_integral_test_types mp_integral_test_types;
	283	#endif
	284
	285	#ifdef TEST1
	286	typedef boost::mpl::list<double, long double> non_integral_test_types;
	287	#elif defined(TEST2)
	288	typedef boost::mpl::list<
7c673cae	289	#if !BOOST_WORKAROUND(BOOST_MSVC, <= 1500)
b32b8144	290	boost::multiprecision::cpp_rational
7c673cae FG	291	#endif
7c673cae FG	292	> non_integral_test_types;
b32b8144 FG	293	#elif defined(TEST3)
	294	typedef boost::mpl::list<
	295	#if !BOOST_WORKAROUND(BOOST_MSVC, <= 1500)
	296	boost::multiprecision::cpp_bin_float_single, boost::multiprecision::cpp_dec_float_50
	297	#endif
	298	> non_integral_test_types;
	299	#endif
	300
7c673cae FG	301	typedef boost::mpl::joint_view<integral_test_types, non_integral_test_types> all_test_types;
7c673cae FG	302
b32b8144 FG	303
	304	template <typename T>
	305	void normalize(polynomial<T> &p)
	306	{
	307	if (leading_coefficient(p) < T(0))
	308	std::transform(p.data().begin(), p.data().end(), p.data().begin(), std::negate<T>());
	309	}
	310
	311	/**
	312	* Note that we do not expect 'pure' gcd algorithms to normalize the result.
	313	* However, the usual public interface function gcd() will do that.
	314	*/
	315
	316	BOOST_AUTO_TEST_SUITE(test_subresultant_gcd)
	317
	318	// This test is just to show that gcd<polynomial<T>>(u, v) is defined (and works) when T is integral and multiprecision.
	319	BOOST_FIXTURE_TEST_CASE_TEMPLATE( gcd_interface, T, mp_integral_test_types, FM2GP_Ex_8_3__1<T> )
	320	{
	321	typedef FM2GP_Ex_8_3__1<T> fixture_type;
	322	polynomial<T> w;
	323	w = gcd(fixture_type::x, fixture_type::y);
	324	normalize(w);
	325	BOOST_CHECK_EQUAL(w, fixture_type::z);
	326	w = gcd(fixture_type::y, fixture_type::x);
	327	normalize(w);
	328	BOOST_CHECK_EQUAL(w, fixture_type::z);
	329	}
	330
	331	// This test is just to show that gcd<polynomial<T>>(u, v) is defined (and works) when T is floating point.
	332	BOOST_FIXTURE_TEST_CASE_TEMPLATE( gcd_float_interface, T, non_integral_test_types, FM2GP_Ex_8_3__1<T> )
	333	{
	334	typedef FM2GP_Ex_8_3__1<T> fixture_type;
	335	polynomial<T> w;
	336	w = gcd(fixture_type::x, fixture_type::y);
	337	normalize(w);
	338	BOOST_CHECK_EQUAL(w, fixture_type::z);
	339	w = gcd(fixture_type::y, fixture_type::x);
	340	normalize(w);
	341	BOOST_CHECK_EQUAL(w, fixture_type::z);
	342	}
	343
	344	// The following tests call subresultant_gcd explicitly to remove any ambiguity
	345	// and to permit testing on single-precision integral types.
	346	BOOST_FIXTURE_TEST_CASE_TEMPLATE( Ex_8_3__1, T, large_integral_test_types, FM2GP_Ex_8_3__1<T> )
	347	{
	348	typedef FM2GP_Ex_8_3__1<T> fixture_type;
	349	polynomial<T> w;
	350	w = subresultant_gcd(fixture_type::x, fixture_type::y);
	351	normalize(w);
	352	BOOST_CHECK_EQUAL(w, fixture_type::z);
	353	w = subresultant_gcd(fixture_type::y, fixture_type::x);
	354	normalize(w);
	355	BOOST_CHECK_EQUAL(w, fixture_type::z);
	356	}
	357
	358	BOOST_FIXTURE_TEST_CASE_TEMPLATE( Ex_8_3__2, T, large_integral_test_types, FM2GP_Ex_8_3__2<T> )
	359	{
	360	typedef FM2GP_Ex_8_3__2<T> fixture_type;
	361	polynomial<T> w;
	362	w = subresultant_gcd(fixture_type::x, fixture_type::y);
	363	normalize(w);
	364	BOOST_CHECK_EQUAL(w, fixture_type::z);
	365	w = subresultant_gcd(fixture_type::y, fixture_type::x);
	366	normalize(w);
367	BOOST_CHECK_EQUAL(w, fixture_type::z);
368	}
369
370	BOOST_FIXTURE_TEST_CASE_TEMPLATE( trivial_int, T, large_integral_test_types, FM2GP_trivial<T> )
371	{
372	typedef FM2GP_trivial<T> fixture_type;
373	polynomial<T> w;
374	w = subresultant_gcd(fixture_type::x, fixture_type::y);
375	normalize(w);
376	BOOST_CHECK_EQUAL(w, fixture_type::z);
377	w = subresultant_gcd(fixture_type::y, fixture_type::x);
378	normalize(w);
379	BOOST_CHECK_EQUAL(w, fixture_type::z);
380	}
381
382	BOOST_AUTO_TEST_SUITE_END()
383
384
7c673cae FG	385	BOOST_AUTO_TEST_CASE_TEMPLATE( test_addition, T, all_test_types )
	386	{
	387	polynomial<T> const a(d3a.begin(), d3a.end());
	388	polynomial<T> const b(d1a.begin(), d1a.end());
	389	polynomial<T> const zero;
11fdf7f2	390
7c673cae FG	391	polynomial<T> result = a + b; // different degree
	392	boost::array<T, 4> tmp = {{8, -5, -4, 3}};
	393	polynomial<T> expected(tmp.begin(), tmp.end());
	394	BOOST_CHECK_EQUAL(result, expected);
	395	BOOST_CHECK_EQUAL(a + zero, a);
	396	BOOST_CHECK_EQUAL(a + b, b + a);
	397	}
	398
	399	BOOST_AUTO_TEST_CASE_TEMPLATE( test_subtraction, T, all_test_types )
	400	{
	401	polynomial<T> const a(d3a.begin(), d3a.end());
	402	polynomial<T> const zero;
	403
	404	BOOST_CHECK_EQUAL(a - T(0), a);
	405	BOOST_CHECK_EQUAL(T(0) - a, -a);
	406	BOOST_CHECK_EQUAL(a - zero, a);
	407	BOOST_CHECK_EQUAL(zero - a, -a);
	408	BOOST_CHECK_EQUAL(a - a, zero);
	409	}
	410
	411	BOOST_AUTO_TEST_CASE_TEMPLATE( test_multiplication, T, all_test_types )
	412	{
	413	polynomial<T> const a(d3a.begin(), d3a.end());
	414	polynomial<T> const b(d1a.begin(), d1a.end());
	415	polynomial<T> const zero;
	416	boost::array<T, 7> const d3a_sq = {{100, -120, -44, 108, -20, -24, 9}};
	417	polynomial<T> const a_sq(d3a_sq.begin(), d3a_sq.end());
11fdf7f2	418
7c673cae FG	419	BOOST_CHECK_EQUAL(a * T(0), zero);
	420	BOOST_CHECK_EQUAL(a * zero, zero);
	421	BOOST_CHECK_EQUAL(zero * T(0), zero);
	422	BOOST_CHECK_EQUAL(zero * zero, zero);
	423	BOOST_CHECK_EQUAL(a * b, b * a);
	424	polynomial<T> aa(a);
	425	aa *= aa;
	426	BOOST_CHECK_EQUAL(aa, a_sq);
	427	BOOST_CHECK_EQUAL(aa, a * a);
	428	}
	429
	430	BOOST_AUTO_TEST_CASE_TEMPLATE( test_arithmetic_relations, T, all_test_types )
	431	{
	432	polynomial<T> const a(d8b.begin(), d8b.end());
	433	polynomial<T> const b(d1a.begin(), d1a.end());
	434
	435	BOOST_CHECK_EQUAL(a * T(2), a + a);
	436	BOOST_CHECK_EQUAL(a - b, -b + a);
	437	BOOST_CHECK_EQUAL(a, (a * a) / a);
	438	BOOST_CHECK_EQUAL(a, (a / a) * a);
	439	}
	440
	441
	442	BOOST_AUTO_TEST_CASE_TEMPLATE(test_non_integral_arithmetic_relations, T, non_integral_test_types )
	443	{
	444	polynomial<T> const a(d8b.begin(), d8b.end());
	445	polynomial<T> const b(d1a.begin(), d1a.end());
11fdf7f2	446
7c673cae FG	447	BOOST_CHECK_EQUAL(a * T(0.5), a / T(2));
	448	}
	449
b32b8144 FG	450	BOOST_AUTO_TEST_CASE_TEMPLATE(test_cont_and_pp, T, integral_test_types)
	451	{
	452	boost::array<polynomial<T>, 4> const q={{
11fdf7f2	453	polynomial<T>(d8.begin(), d8.end()),
b32b8144 FG	454	polynomial<T>(d8b.begin(), d8b.end()),
	455	polynomial<T>(d3a.begin(), d3a.end()),
	456	polynomial<T>(d3b.begin(), d3b.end())
	457	}};
	458	for (std::size_t i = 0; i < q.size(); i++)
	459	{
	460	BOOST_CHECK_EQUAL(q[i], content(q[i]) * primitive_part(q[i]));
	461	BOOST_CHECK_EQUAL(primitive_part(q[i]), primitive_part(q[i], content(q[i])));
	462	}
	463
	464	polynomial<T> const zero;
	465	BOOST_CHECK_EQUAL(primitive_part(zero), zero);
	466	BOOST_CHECK_EQUAL(content(zero), T(0));
	467	}
	468
7c673cae FG	469	BOOST_AUTO_TEST_CASE_TEMPLATE( test_self_multiply_assign, T, all_test_types )
	470	{
	471	polynomial<T> a(d3a.begin(), d3a.end());
	472	polynomial<T> const b(a);
	473	boost::array<double, 7> const d3a_sq = {{100, -120, -44, 108, -20, -24, 9}};
	474	polynomial<T> const asq(d3a_sq.begin(), d3a_sq.end());
	475
	476	a *= a;
	477
	478	BOOST_CHECK_EQUAL(a, b*b);
	479	BOOST_CHECK_EQUAL(a, asq);
	480
	481	a *= a;
	482
	483	BOOST_CHECK_EQUAL(a, bbb*b);
	484	}
	485
b32b8144	486
7c673cae FG	487	BOOST_AUTO_TEST_CASE_TEMPLATE(test_right_shift, T, all_test_types )
	488	{
	489	polynomial<T> a(d8b.begin(), d8b.end());
	490	polynomial<T> const aa(a);
	491	polynomial<T> const b(d8b.begin() + 1, d8b.end());
	492	polynomial<T> const c(d8b.begin() + 5, d8b.end());
	493	a >>= 0u;
	494	BOOST_CHECK_EQUAL(a, aa);
	495	a >>= 1u;
	496	BOOST_CHECK_EQUAL(a, b);
	497	a = a >> 4u;
	498	BOOST_CHECK_EQUAL(a, c);
	499	}
	500
	501
	502	BOOST_AUTO_TEST_CASE_TEMPLATE(test_left_shift, T, all_test_types )
	503	{
	504	polynomial<T> a(d0a.begin(), d0a.end());
	505	polynomial<T> const aa(a);
	506	polynomial<T> const b(d0a1.begin(), d0a1.end());
	507	polynomial<T> const c(d0a5.begin(), d0a5.end());
	508	a <<= 0u;
11fdf7f2	509	BOOST_CHECK_EQUAL(a, aa);
7c673cae FG	510	a <<= 1u;
	511	BOOST_CHECK_EQUAL(a, b);
	512	a = a << 4u;
	513	BOOST_CHECK_EQUAL(a, c);
	514	polynomial<T> zero;
	515	// Multiplying zero by x should still be zero.
	516	zero <<= 1u;
	517	BOOST_CHECK_EQUAL(zero, zero_element(multiplies< polynomial<T> >()));
	518	}
	519
	520
	521	BOOST_AUTO_TEST_CASE_TEMPLATE(test_odd_even, T, all_test_types)
	522	{
	523	polynomial<T> const zero;
	524	BOOST_CHECK_EQUAL(odd(zero), false);
	525	BOOST_CHECK_EQUAL(even(zero), true);
	526	polynomial<T> const a(d0a.begin(), d0a.end());
	527	BOOST_CHECK_EQUAL(odd(a), true);
	528	BOOST_CHECK_EQUAL(even(a), false);
	529	polynomial<T> const b(d0a1.begin(), d0a1.end());
	530	BOOST_CHECK_EQUAL(odd(b), false);
	531	BOOST_CHECK_EQUAL(even(b), true);
	532	}
	533
b32b8144	534	// NOTE: Slightly unexpected: this unit test passes even when T = char.
7c673cae FG	535	BOOST_AUTO_TEST_CASE_TEMPLATE( test_pow, T, all_test_types )
7c673cae FG	536	{
b32b8144 FG	537	if (std::numeric_limits<T>::digits < 32)
b32b8144 FG	538	return; // Invokes undefined behaviour
7c673cae FG	539	polynomial<T> a(d3a.begin(), d3a.end());
	540	polynomial<T> const one(T(1));
	541	boost::array<double, 7> const d3a_sqr = {{100, -120, -44, 108, -20, -24, 9}};
	542	boost::array<double, 10> const d3a_cub =
	543	{{1000, -1800, -120, 2124, -1032, -684, 638, -18, -108, 27}};
	544	polynomial<T> const asqr(d3a_sqr.begin(), d3a_sqr.end());
	545	polynomial<T> const acub(d3a_cub.begin(), d3a_cub.end());
	546
	547	BOOST_CHECK_EQUAL(pow(a, 0), one);
	548	BOOST_CHECK_EQUAL(pow(a, 1), a);
	549	BOOST_CHECK_EQUAL(pow(a, 2), asqr);
	550	BOOST_CHECK_EQUAL(pow(a, 3), acub);
	551	BOOST_CHECK_EQUAL(pow(a, 4), pow(asqr, 2));
	552	BOOST_CHECK_EQUAL(pow(a, 5), asqr * acub);
	553	BOOST_CHECK_EQUAL(pow(a, 6), pow(acub, 2));
	554	BOOST_CHECK_EQUAL(pow(a, 7), acub * acub * a);
	555
	556	BOOST_CHECK_THROW(pow(a, -1), std::domain_error);
	557	BOOST_CHECK_EQUAL(pow(one, 137), one);
	558	}
	559
	560
	561	BOOST_AUTO_TEST_CASE_TEMPLATE(test_bool, T, all_test_types)
	562	{
	563	polynomial<T> const zero;
	564	polynomial<T> const a(d0a.begin(), d0a.end());
	565	BOOST_CHECK_EQUAL(bool(zero), false);
	566	BOOST_CHECK_EQUAL(bool(a), true);
	567	}
	568
	569
	570	BOOST_AUTO_TEST_CASE_TEMPLATE(test_set_zero, T, all_test_types)
	571	{
	572	polynomial<T> const zero;
	573	polynomial<T> a(d0a.begin(), d0a.end());
	574	a.set_zero();
	575	BOOST_CHECK_EQUAL(a, zero);
	576	a.set_zero(); // Ensure that setting zero to zero is a no-op.
	577	BOOST_CHECK_EQUAL(a, zero);
	578	}
b32b8144 FG	579
	580
	581	BOOST_AUTO_TEST_CASE_TEMPLATE(test_leading_coefficient, T, all_test_types)
	582	{
	583	polynomial<T> const zero;
	584	BOOST_CHECK_EQUAL(leading_coefficient(zero), T(0));
	585	polynomial<T> a(d0a.begin(), d0a.end());
	586	BOOST_CHECK_EQUAL(leading_coefficient(a), T(d0a.back()));
	587	}
92f5a8d4 TL	588
	589	#if !defined(BOOST_NO_CXX11_RVALUE_REFERENCES) && !defined(BOOST_NO_CXX11_UNIFIED_INITIALIZATION_SYNTAX)
	590	BOOST_AUTO_TEST_CASE_TEMPLATE(test_prime, T, all_test_types)
	591	{
	592	std::vector<T> d{1,1,1,1,1};
	593	polynomial<T> p(std::move(d));
	594	polynomial<T> q = p.prime();
	595	BOOST_CHECK_EQUAL(q(0), T(1));
	596
	597	for (size_t i = 0; i < q.size(); ++i)
	598	{
	599	BOOST_CHECK_EQUAL(q[i], i+1);
	600	}
	601
	602	polynomial<T> P = p.integrate();
	603	BOOST_CHECK_EQUAL(P(0), T(0));
	604	for (size_t i = 1; i < P.size(); ++i)
	605	{
	606	BOOST_CHECK_EQUAL(P[i], 1/static_cast<T>(i));
	607	}
	608
	609	polynomial<T> empty;
	610	q = empty.prime();
	611	BOOST_CHECK_EQUAL(q.size(), 0);
	612
	613	}
	614	#endif