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/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
  20 using namespace boost::math::tools;
  21 using namespace std;
  22
  23 template <typename T>
  24 struct answer
  25 {
  26     answer(std::pair< polynomial<T>, polynomial<T> > const &x) :
  27     quotient(x.first), remainder(x.second) {}
  28
  29     polynomial<T> quotient;
  30     polynomial<T> remainder;
  31 };
  32
  33 boost::array<double, 4> const d3a = {{10, -6, -4, 3}};
  34 boost::array<double, 4> const d3b = {{-7, 5, 6, 1}};
  35 boost::array<double, 4> const d3c = {{10.0/3.0, -2.0, -4.0/3.0, 1.0}};
  36 boost::array<double, 2> const d1a = {{-2, 1}};
  37 boost::array<double, 3> const d2a = {{-2, 2, 3}};
  38 boost::array<double, 3> const d2b = {{-7, 5, 6}};
  39 boost::array<double, 3> const d2c = {{31, -21, -22}};
  40 boost::array<double, 1> const d0a = {{6}};
  41 boost::array<double, 2> const d0a1 = {{0, 6}};
  42 boost::array<double, 6> const d0a5 = {{0, 0, 0, 0, 0, 6}};
  43 boost::array<double, 1> const d0b = {{3}};
  44
  45 boost::array<int, 9> const d8 = {{-5, 2, 8, -3, -3, 0, 1, 0, 1}};
  46 boost::array<int, 9> const d8b = {{0, 2, 8, -3, -3, 0, 1, 0, 1}};
  47 boost::array<int, 7> const d6 = {{21, -9, -4, 0, 5, 0, 3}};
  48 boost::array<int, 3> const d2 = {{-6, 0, 9}};
  49 boost::array<int, 6> const d5 = {{-9, 0, 3, 0, -15}};
  50
  51
  52 BOOST_AUTO_TEST_CASE( test_construction )
  53 {
  54     polynomial<double> const a(d3a.begin(), d3a.end());
  55     polynomial<double> const b(d3a.begin(), 3);
  56     BOOST_CHECK_EQUAL(a, b);
  57 }
  58
  59
  60 #if !defined(BOOST_NO_CXX11_HDR_INITIALIZER_LIST) && !BOOST_WORKAROUND(BOOST_GCC_VERSION, < 40500)
  61 BOOST_AUTO_TEST_CASE( test_initializer_list_construction )
  62 {
  63     polynomial<double> a(begin(d3a), end(d3a));
  64     polynomial<double> b = {10, -6, -4, 3};
  65     polynomial<double> c{{10, -6, -4, 3}};
  66     polynomial<double> d{{10, -6, -4, 3, 0, 0}};
  67     BOOST_CHECK_EQUAL(a, b);
  68     BOOST_CHECK_EQUAL(b, c);
  69     BOOST_CHECK_EQUAL(d.degree(), 3u);
  70 }
  71
  72 BOOST_AUTO_TEST_CASE( test_initializer_list_assignment )
  73 {
  74     polynomial<double> a(begin(d3a), end(d3a));
  75     polynomial<double> b;
  76     b = {10, -6, -4, 3, 0, 0};
  77     BOOST_CHECK_EQUAL(b.degree(), 3u);
  78     BOOST_CHECK_EQUAL(a, b);
  79 }
  80 #endif
  81
  82
  83 BOOST_AUTO_TEST_CASE( test_degree )
  84 {
  85     polynomial<double> const zero;
  86     polynomial<double> const a(d3a.begin(), d3a.end());
  87     BOOST_CHECK_THROW(zero.degree(), std::logic_error);
  88     BOOST_CHECK_EQUAL(a.degree(), 3u);
  89 }
  90
  91
  92 BOOST_AUTO_TEST_CASE( test_division_over_field )
  93 {
  94     polynomial<double> const a(d3a.begin(), d3a.end());
  95     polynomial<double> const b(d1a.begin(), d1a.end());
  96     polynomial<double> const q(d2a.begin(), d2a.end());
  97     polynomial<double> const r(d0a.begin(), d0a.end());
  98     polynomial<double> const c(d3b.begin(), d3b.end());
  99     polynomial<double> const d(d2b.begin(), d2b.end());
 100     polynomial<double> const e(d2c.begin(), d2c.end());
 101     polynomial<double> const f(d0b.begin(), d0b.end());
 102     polynomial<double> const g(d3c.begin(), d3c.end());
 103     polynomial<double> const zero;
 104     polynomial<double> const one(1.0);
 105
 106     answer<double> result = quotient_remainder(a, b);
 107     BOOST_CHECK_EQUAL(result.quotient, q);
 108     BOOST_CHECK_EQUAL(result.remainder, r);
 109     BOOST_CHECK_EQUAL(a, q * b + r); // Sanity check.
 110
 111     result = quotient_remainder(a, c);
 112     BOOST_CHECK_EQUAL(result.quotient, f);
 113     BOOST_CHECK_EQUAL(result.remainder, e);
 114     BOOST_CHECK_EQUAL(a, f * c + e); // Sanity check.
 115
 116     result = quotient_remainder(a, f);
 117     BOOST_CHECK_EQUAL(result.quotient, g);
 118     BOOST_CHECK_EQUAL(result.remainder, zero);
 119     BOOST_CHECK_EQUAL(a, g * f + zero); // Sanity check.
 120     // Check that division by a regular number gives the same result.
 121     BOOST_CHECK_EQUAL(a / 3.0, g);
 122     BOOST_CHECK_EQUAL(a % 3.0, zero);
 123
 124     // Sanity checks.
 125     BOOST_CHECK_EQUAL(a / a, one);
 126     BOOST_CHECK_EQUAL(a % a, zero);
 127     // BOOST_CHECK_EQUAL(zero / zero, zero); // TODO
 128 }
 129
 130 BOOST_AUTO_TEST_CASE( test_division_over_ufd )
 131 {
 132     polynomial<int> const zero;
 133     polynomial<int> const one(1);
 134     polynomial<int> const aa(d8.begin(), d8.end());
 135     polynomial<int> const bb(d6.begin(), d6.end());
 136     polynomial<int> const q(d2.begin(), d2.end());
 137     polynomial<int> const r(d5.begin(), d5.end());
 138
 139     answer<int> result = quotient_remainder(aa, bb);
 140     BOOST_CHECK_EQUAL(result.quotient, q);
 141     BOOST_CHECK_EQUAL(result.remainder, r);
 142
 143     // Sanity checks.
 144     BOOST_CHECK_EQUAL(aa / aa, one);
 145     BOOST_CHECK_EQUAL(aa % aa, zero);
 146 }
 147
 148
 149 BOOST_AUTO_TEST_CASE( test_gcd )
 150 {
 151     /* NOTE: Euclidean gcd is not yet customized to return THE greatest
 152      * common polynomial divisor. If d is THE greatest common divisior of u and
 153      * v, then gcd(u, v) will return d or -d according to the algorithm.
 154      * By convention, it should return d, as for example Maxima and Wolfram
 155      * Alpha do.
 156      * This test is an example of the fact that it returns -d.
 157      */
 158     boost::array<double, 9> const d8 = {{105, 278, -88, -56, 16}};
 159     boost::array<double, 7> const d6 = {{70, 232, -44, -64, 16}};
 160     boost::array<double, 7> const d2 = {{-35, 24, -4}};
 161     polynomial<double> const u(d8.begin(), d8.end());
 162     polynomial<double> const v(d6.begin(), d6.end());
 163     polynomial<double> const w(d2.begin(), d2.end());
 164     polynomial<double> const d = boost::math::gcd(u, v);
 165     BOOST_CHECK_EQUAL(w, d);
 166 }
 167
 168 // Sanity checks to make sure I didn't break it.
 169 typedef boost::mpl::list<int, long
 170 #if !BOOST_WORKAROUND(BOOST_MSVC, <= 1500)
 171    , boost::multiprecision::cpp_int
 172 #endif
 173 > integral_test_types;
 174 typedef boost::mpl::list<double
 175 #if !BOOST_WORKAROUND(BOOST_MSVC, <= 1500)
 176    , boost::multiprecision::cpp_rational, boost::multiprecision::cpp_bin_float_single, boost::multiprecision::cpp_dec_float_50
 177 #endif
 178 > non_integral_test_types;
 179 typedef boost::mpl::joint_view<integral_test_types, non_integral_test_types> all_test_types;
 180
 181 BOOST_AUTO_TEST_CASE_TEMPLATE( test_addition, T, all_test_types )
 182 {
 183     polynomial<T> const a(d3a.begin(), d3a.end());
 184     polynomial<T> const b(d1a.begin(), d1a.end());
 185     polynomial<T> const zero;
 186
 187     polynomial<T> result = a + b; // different degree
 188     boost::array<T, 4> tmp = {{8, -5, -4, 3}};
 189     polynomial<T> expected(tmp.begin(), tmp.end());
 190     BOOST_CHECK_EQUAL(result, expected);
 191     BOOST_CHECK_EQUAL(a + zero, a);
 192     BOOST_CHECK_EQUAL(a + b, b + a);
 193 }
 194
 195 BOOST_AUTO_TEST_CASE_TEMPLATE( test_subtraction, T, all_test_types )
 196 {
 197     polynomial<T> const a(d3a.begin(), d3a.end());
 198     polynomial<T> const zero;
 199
 200     BOOST_CHECK_EQUAL(a - T(0), a);
 201     BOOST_CHECK_EQUAL(T(0) - a, -a);
 202     BOOST_CHECK_EQUAL(a - zero, a);
 203     BOOST_CHECK_EQUAL(zero - a, -a);
 204     BOOST_CHECK_EQUAL(a - a, zero);
 205 }
 206
 207 BOOST_AUTO_TEST_CASE_TEMPLATE( test_multiplication, T, all_test_types )
 208 {
 209     polynomial<T> const a(d3a.begin(), d3a.end());
 210     polynomial<T> const b(d1a.begin(), d1a.end());
 211     polynomial<T> const zero;
 212     boost::array<T, 7> const d3a_sq = {{100, -120, -44, 108, -20, -24, 9}};
 213     polynomial<T> const a_sq(d3a_sq.begin(), d3a_sq.end());
 214
 215     BOOST_CHECK_EQUAL(a * T(0), zero);
 216     BOOST_CHECK_EQUAL(a * zero, zero);
 217     BOOST_CHECK_EQUAL(zero * T(0), zero);
 218     BOOST_CHECK_EQUAL(zero * zero, zero);
 219     BOOST_CHECK_EQUAL(a * b, b * a);
 220     polynomial<T> aa(a);
 221     aa *= aa;
 222     BOOST_CHECK_EQUAL(aa, a_sq);
 223     BOOST_CHECK_EQUAL(aa, a * a);
 224 }
 225
 226 BOOST_AUTO_TEST_CASE_TEMPLATE( test_arithmetic_relations, T, all_test_types )
 227 {
 228     polynomial<T> const a(d8b.begin(), d8b.end());
 229     polynomial<T> const b(d1a.begin(), d1a.end());
 230
 231     BOOST_CHECK_EQUAL(a * T(2), a + a);
 232     BOOST_CHECK_EQUAL(a - b, -b + a);
 233     BOOST_CHECK_EQUAL(a, (a * a) / a);
 234     BOOST_CHECK_EQUAL(a, (a / a) * a);
 235 }
 236
 237
 238 BOOST_AUTO_TEST_CASE_TEMPLATE(test_non_integral_arithmetic_relations, T, non_integral_test_types )
 239 {
 240     polynomial<T> const a(d8b.begin(), d8b.end());
 241     polynomial<T> const b(d1a.begin(), d1a.end());
 242
 243     BOOST_CHECK_EQUAL(a * T(0.5), a / T(2));
 244 }
 245
 246 BOOST_AUTO_TEST_CASE_TEMPLATE( test_self_multiply_assign, T, all_test_types )
 247 {
 248     polynomial<T> a(d3a.begin(), d3a.end());
 249     polynomial<T> const b(a);
 250     boost::array<double, 7> const d3a_sq = {{100, -120, -44, 108, -20, -24, 9}};
 251     polynomial<T> const asq(d3a_sq.begin(), d3a_sq.end());
 252
 253     a *= a;
 254
 255     BOOST_CHECK_EQUAL(a, b*b);
 256     BOOST_CHECK_EQUAL(a, asq);
 257
 258     a *= a;
 259
 260     BOOST_CHECK_EQUAL(a, b*b*b*b);
 261 }
 262
 263 BOOST_AUTO_TEST_CASE_TEMPLATE(test_right_shift, T, all_test_types )
 264 {
 265     polynomial<T> a(d8b.begin(), d8b.end());
 266     polynomial<T> const aa(a);
 267     polynomial<T> const b(d8b.begin() + 1, d8b.end());
 268     polynomial<T> const c(d8b.begin() + 5, d8b.end());
 269     a >>= 0u;
 270     BOOST_CHECK_EQUAL(a, aa);
 271     a >>= 1u;
 272     BOOST_CHECK_EQUAL(a, b);
 273     a = a >> 4u;
 274     BOOST_CHECK_EQUAL(a, c);
 275 }
 276
 277
 278 BOOST_AUTO_TEST_CASE_TEMPLATE(test_left_shift, T, all_test_types )
 279 {
 280     polynomial<T> a(d0a.begin(), d0a.end());
 281     polynomial<T> const aa(a);
 282     polynomial<T> const b(d0a1.begin(), d0a1.end());
 283     polynomial<T> const c(d0a5.begin(), d0a5.end());
 284     a <<= 0u;
 285     BOOST_CHECK_EQUAL(a, aa);
 286     a <<= 1u;
 287     BOOST_CHECK_EQUAL(a, b);
 288     a = a << 4u;
 289     BOOST_CHECK_EQUAL(a, c);
 290     polynomial<T> zero;
 291     // Multiplying zero by x should still be zero.
 292     zero <<= 1u;
 293     BOOST_CHECK_EQUAL(zero, zero_element(multiplies< polynomial<T> >()));
 294 }
 295
 296
 297 BOOST_AUTO_TEST_CASE_TEMPLATE(test_odd_even, T, all_test_types)
 298 {
 299     polynomial<T> const zero;
 300     BOOST_CHECK_EQUAL(odd(zero), false);
 301     BOOST_CHECK_EQUAL(even(zero), true);
 302     polynomial<T> const a(d0a.begin(), d0a.end());
 303     BOOST_CHECK_EQUAL(odd(a), true);
 304     BOOST_CHECK_EQUAL(even(a), false);
 305     polynomial<T> const b(d0a1.begin(), d0a1.end());
 306     BOOST_CHECK_EQUAL(odd(b), false);
 307     BOOST_CHECK_EQUAL(even(b), true);
 308 }
 309
 310
 311 BOOST_AUTO_TEST_CASE_TEMPLATE( test_pow, T, all_test_types )
 312 {
 313     polynomial<T> a(d3a.begin(), d3a.end());
 314     polynomial<T> const one(T(1));
 315     boost::array<double, 7> const d3a_sqr = {{100, -120, -44, 108, -20, -24, 9}};
 316     boost::array<double, 10> const d3a_cub =
 317         {{1000, -1800, -120, 2124, -1032, -684, 638, -18, -108, 27}};
 318     polynomial<T> const asqr(d3a_sqr.begin(), d3a_sqr.end());
 319     polynomial<T> const acub(d3a_cub.begin(), d3a_cub.end());
 320
 321     BOOST_CHECK_EQUAL(pow(a, 0), one);
 322     BOOST_CHECK_EQUAL(pow(a, 1), a);
 323     BOOST_CHECK_EQUAL(pow(a, 2), asqr);
 324     BOOST_CHECK_EQUAL(pow(a, 3), acub);
 325     BOOST_CHECK_EQUAL(pow(a, 4), pow(asqr, 2));
 326     BOOST_CHECK_EQUAL(pow(a, 5), asqr * acub);
 327     BOOST_CHECK_EQUAL(pow(a, 6), pow(acub, 2));
 328     BOOST_CHECK_EQUAL(pow(a, 7), acub * acub * a);
 329
 330     BOOST_CHECK_THROW(pow(a, -1), std::domain_error);
 331     BOOST_CHECK_EQUAL(pow(one, 137), one);
 332 }
 333
 334
 335 BOOST_AUTO_TEST_CASE_TEMPLATE(test_bool, T, all_test_types)
 336 {
 337     polynomial<T> const zero;
 338     polynomial<T> const a(d0a.begin(), d0a.end());
 339     BOOST_CHECK_EQUAL(bool(zero), false);
 340     BOOST_CHECK_EQUAL(bool(a), true);
 341 }
 342
 343
 344 BOOST_AUTO_TEST_CASE_TEMPLATE(test_set_zero, T, all_test_types)
 345 {
 346     polynomial<T> const zero;
 347     polynomial<T> a(d0a.begin(), d0a.end());
 348     a.set_zero();
 349     BOOST_CHECK_EQUAL(a, zero);
 350     a.set_zero(); // Ensure that setting zero to zero is a no-op.
 351     BOOST_CHECK_EQUAL(a, zero);
 352 }