ceph/src/boost/libs/math/test/test_polynomial.cpp

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