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 {
answeranswer36 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
46 boost::array<double, 2> const d1a = {{-2, 1}};
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}};
50
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}};
54
55
56
BOOST_AUTO_TEST_CASE(trivial)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
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
74
BOOST_AUTO_TEST_CASE(test_construction)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
82 #ifdef BOOST_MATH_HAS_IS_CONST_ITERABLE
83
84 #include <list>
85 #include <array>
86
BOOST_AUTO_TEST_CASE(test_range_construction)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
101
102 #if !defined(BOOST_NO_CXX11_HDR_INITIALIZER_LIST) && !BOOST_WORKAROUND(BOOST_GCC_VERSION, < 40500)
BOOST_AUTO_TEST_CASE(test_initializer_list_construction)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};
107 polynomial<double> c{10, -6, -4, 3};
108 polynomial<double> d{10, -6, -4, 3, 0, 0};
109 BOOST_CHECK_EQUAL(a, b);
110 BOOST_CHECK_EQUAL(b, c);
111 BOOST_CHECK_EQUAL(d.degree(), 3u);
112 }
113
BOOST_AUTO_TEST_CASE(test_initializer_list_assignment)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
BOOST_AUTO_TEST_CASE(test_degree)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
BOOST_AUTO_TEST_CASE(test_division_over_field)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.
152
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.
157
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
BOOST_AUTO_TEST_CASE(test_division_over_ufd)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());
180
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
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;
198
FM2GP_Ex_8_3__1FM2GP_Ex_8_3__1199 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;
216
FM2GP_Ex_8_3__2FM2GP_Ex_8_3__2217 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 };
227
228
229 template <typename T>
230 struct FM2GP_mixed
231 {
232 polynomial<T> x;
233 polynomial<T> y;
234 polynomial<T> z;
235
FM2GP_mixedFM2GP_mixed236 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;
254
FM2GP_trivialFM2GP_trivial255 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 };
265
266 // Sanity checks to make sure I didn't break it.
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<
273 #if !BOOST_WORKAROUND(BOOST_MSVC, <= 1500)
274 boost::multiprecision::cpp_int
275 #endif
276 > integral_test_types;
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<
289 #if !BOOST_WORKAROUND(BOOST_MSVC, <= 1500)
290 boost::multiprecision::cpp_rational
291 #endif
292 > non_integral_test_types;
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
301 typedef boost::mpl::joint_view<integral_test_types, non_integral_test_types> all_test_types;
302
303
304 template <typename T>
normalize(polynomial<T> & p)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.
BOOST_FIXTURE_TEST_CASE_TEMPLATE(gcd_interface,T,mp_integral_test_types,FM2GP_Ex_8_3__1<T>)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.
BOOST_FIXTURE_TEST_CASE_TEMPLATE(gcd_float_interface,T,non_integral_test_types,FM2GP_Ex_8_3__1<T>)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.
BOOST_FIXTURE_TEST_CASE_TEMPLATE(Ex_8_3__1,T,large_integral_test_types,FM2GP_Ex_8_3__1<T>)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
BOOST_FIXTURE_TEST_CASE_TEMPLATE(Ex_8_3__2,T,large_integral_test_types,FM2GP_Ex_8_3__2<T>)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
BOOST_FIXTURE_TEST_CASE_TEMPLATE(trivial_int,T,large_integral_test_types,FM2GP_trivial<T>)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
BOOST_AUTO_TEST_CASE_TEMPLATE(test_addition,T,all_test_types)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;
390
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
BOOST_AUTO_TEST_CASE_TEMPLATE(test_subtraction,T,all_test_types)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
BOOST_AUTO_TEST_CASE_TEMPLATE(test_multiplication,T,all_test_types)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());
418
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
BOOST_AUTO_TEST_CASE_TEMPLATE(test_arithmetic_relations,T,all_test_types)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
BOOST_AUTO_TEST_CASE_TEMPLATE(test_non_integral_arithmetic_relations,T,non_integral_test_types)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());
446
447 BOOST_CHECK_EQUAL(a * T(0.5), a / T(2));
448 }
449
BOOST_AUTO_TEST_CASE_TEMPLATE(test_cont_and_pp,T,integral_test_types)450 BOOST_AUTO_TEST_CASE_TEMPLATE(test_cont_and_pp, T, integral_test_types)
451 {
452 boost::array<polynomial<T>, 4> const q={{
453 polynomial<T>(d8.begin(), d8.end()),
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
BOOST_AUTO_TEST_CASE_TEMPLATE(test_self_multiply_assign,T,all_test_types)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, b*b*b*b);
484 }
485
486
BOOST_AUTO_TEST_CASE_TEMPLATE(test_right_shift,T,all_test_types)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
BOOST_AUTO_TEST_CASE_TEMPLATE(test_left_shift,T,all_test_types)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;
509 BOOST_CHECK_EQUAL(a, aa);
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
BOOST_AUTO_TEST_CASE_TEMPLATE(test_odd_even,T,all_test_types)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
534 // NOTE: Slightly unexpected: this unit test passes even when T = char.
BOOST_AUTO_TEST_CASE_TEMPLATE(test_pow,T,all_test_types)535 BOOST_AUTO_TEST_CASE_TEMPLATE( test_pow, T, all_test_types )
536 {
537 if (std::numeric_limits<T>::digits < 32)
538 return; // Invokes undefined behaviour
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
BOOST_AUTO_TEST_CASE_TEMPLATE(test_bool,T,all_test_types)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
BOOST_AUTO_TEST_CASE_TEMPLATE(test_set_zero,T,all_test_types)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 }
579
580
BOOST_AUTO_TEST_CASE_TEMPLATE(test_leading_coefficient,T,all_test_types)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 }
588
589 #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)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
615