]>
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 | ||
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 | } |