[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/math/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 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, 4> const d3c = {{10.0/3.0, -2.0, -4.0/3.0, 1.0}};
boost::array<double, 2> const d1a = {{-2, 1}};
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 d0a = {{6}};
boost::array<double, 2> const d0a1 = {{0, 6}};
boost::array<double, 6> const d0a5 = {{0, 0, 0, 0, 0, 6}};
boost::array<double, 1> const d0b = {{3}};

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::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(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_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);
}


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