1 /* Copyright 2009-2016 Francesco Biscani (bluescarni@gmail.com)
2
3 This file is part of the Piranha library.
4
5 The Piranha library is free software; you can redistribute it and/or modify
6 it under the terms of either:
7
8 * the GNU Lesser General Public License as published by the Free
9 Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11
12 or
13
14 * the GNU General Public License as published by the Free Software
15 Foundation; either version 3 of the License, or (at your option) any
16 later version.
17
18 or both in parallel, as here.
19
20 The Piranha library is distributed in the hope that it will be useful, but
21 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
22 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
23 for more details.
24
25 You should have received copies of the GNU General Public License and the
26 GNU Lesser General Public License along with the Piranha library. If not,
27 see https://www.gnu.org/licenses/. */
28
29 #include "../src/divisor_series.hpp"
30
31 #define BOOST_TEST_MODULE divisor_series_01_test
32 #include <boost/test/included/unit_test.hpp>
33
34 #include <boost/lexical_cast.hpp>
35 #include <boost/mpl/for_each.hpp>
36 #include <boost/mpl/vector.hpp>
37 #include <limits>
38 #include <stdexcept>
39 #include <string>
40 #include <type_traits>
41
42 #include "../src/divisor.hpp"
43 #include "../src/exceptions.hpp"
44 #include "../src/init.hpp"
45 #include "../src/invert.hpp"
46 #include "../src/kronecker_monomial.hpp"
47 #include "../src/math.hpp"
48 #include "../src/monomial.hpp"
49 #include "../src/mp_integer.hpp"
50 #include "../src/mp_rational.hpp"
51 #include "../src/poisson_series.hpp"
52 #include "../src/polynomial.hpp"
53 #include "../src/pow.hpp"
54 #include "../src/real.hpp"
55 #include "../src/type_traits.hpp"
56
57 using namespace piranha;
58
59 using cf_types = boost::mpl::vector<double, integer, real, rational, polynomial<rational, monomial<int>>>;
60 using expo_types = boost::mpl::vector<short, int, long, integer>;
61
62 struct test_00_tester {
63 template <typename T>
operator ()test_00_tester64 void operator()(const T &)
65 {
66 using s_type = divisor_series<T, divisor<short>>;
67 s_type s0{3};
68 // Just test some math operations and common functionalities.
69 BOOST_CHECK_EQUAL(s0 + s0, 6);
70 BOOST_CHECK_EQUAL(s0 * s0, 9);
71 BOOST_CHECK_EQUAL(s0 * 4, 12);
72 BOOST_CHECK_EQUAL(4 * s0, 12);
73 BOOST_CHECK_EQUAL(math::pow(s0, 3), 27);
74 BOOST_CHECK_EQUAL(math::cos(s_type{0}), 1);
75 BOOST_CHECK_EQUAL(math::sin(s_type{0}), 0);
76 BOOST_CHECK_EQUAL(math::evaluate<int>(math::pow(s0, 3), {{"x", 4}}), 27);
77 BOOST_CHECK(is_differentiable<s_type>::value);
78 BOOST_CHECK_EQUAL(s_type{1}.partial("x"), 0);
79 if (std::is_base_of<detail::polynomial_tag, T>::value) {
80 BOOST_CHECK((has_subs<s_type, s_type>::value));
81 BOOST_CHECK((has_subs<s_type, int>::value));
82 BOOST_CHECK((has_subs<s_type, integer>::value));
83 }
84 BOOST_CHECK((!has_subs<s_type, std::string>::value));
85 if (std::is_base_of<detail::polynomial_tag, T>::value) {
86 BOOST_CHECK((has_ipow_subs<s_type, s_type>::value));
87 BOOST_CHECK((has_ipow_subs<s_type, int>::value));
88 BOOST_CHECK((has_ipow_subs<s_type, integer>::value));
89 }
90 BOOST_CHECK((!has_ipow_subs<s_type, std::string>::value));
91 }
92 };
93
BOOST_AUTO_TEST_CASE(divisor_series_test_00)94 BOOST_AUTO_TEST_CASE(divisor_series_test_00)
95 {
96 init();
97 boost::mpl::for_each<cf_types>(test_00_tester());
98 }
99
100 struct partial_tester {
101 template <typename T>
operator ()partial_tester102 void operator()(const T &)
103 {
104 using p_type = polynomial<rational, monomial<int>>;
105 using s_type = divisor_series<p_type, divisor<T>>;
106 s_type x{"x"}, y{"y"}, z{"z"};
107 // First with variables only in the divisors.
108 auto s0 = math::invert(x + y - 2 * z);
109 BOOST_CHECK((std::is_same<s_type, decltype(s0.partial("x"))>::value));
110 BOOST_CHECK((std::is_same<s_type, decltype(math::partial(s0, "x"))>::value));
111 BOOST_CHECK_EQUAL(s0.partial("x"), -s0 * s0);
112 BOOST_CHECK_EQUAL(math::partial(s0, "x"), -s0 * s0);
113 BOOST_CHECK_EQUAL(s0.partial("z"), 2 * s0 * s0);
114 auto s1 = s0 * s0;
115 BOOST_CHECK_EQUAL(s1.partial("x"), -2 * s0 * s1);
116 BOOST_CHECK_EQUAL(s1.partial("z"), 4 * s0 * s1);
117 auto s2 = math::invert(x - y);
118 auto s3 = s0 * s2;
119 BOOST_CHECK_EQUAL(s3.partial("x"), -s0 * s0 * s2 - s0 * s2 * s2);
120 auto s4 = math::invert(x);
121 auto s5 = s0 * s2 * s4;
122 BOOST_CHECK_EQUAL(s5.partial("x"), -s0 * s0 * s2 * s4 - s0 * s2 * s2 * s4 - s0 * s2 * s4 * s4);
123 BOOST_CHECK_EQUAL(s5.partial("z"), 2 * s0 * s0 * s2 * s4);
124 auto s6 = s0 * s0 * s2 * s4;
125 BOOST_CHECK_EQUAL(s6.partial("x"),
126 -2 * s0 * s0 * s0 * s2 * s4 - s0 * s0 * s2 * s2 * s4 - s0 * s0 * s2 * s4 * s4);
127 // Variables only in the coefficients.
128 auto s7 = s2 * s4 * (x * x / 5 + y - 3 * z);
129 BOOST_CHECK_EQUAL(s7.partial("z"), s2 * s4 * -3);
130 auto s8 = s2 * s4 * (x * x / 5 + y - 3 * z) + z * s2 * s4 * y;
131 BOOST_CHECK_EQUAL(s8.partial("z"), s2 * s4 * -3 + s2 * s4 * y);
132 BOOST_CHECK_EQUAL((x * x * math::invert(z)).partial("x"), 2 * x * math::invert(z));
133 // This excercises the presense of an additional divisor variable with a zero multiplier.
134 BOOST_CHECK_EQUAL((x * x * math::invert(z) + s4 - s4).partial("x"), 2 * x * math::invert(z));
135 // Variables both in the coefficients and in the divisors.
136 auto s9 = x * s2;
137 BOOST_CHECK_EQUAL(s9.partial("x"), s2 - x * s2 * s2);
138 BOOST_CHECK_EQUAL(math::partial(s9, "x"), s2 - x * s2 * s2);
139 auto s10 = x * s2 * s4;
140 BOOST_CHECK_EQUAL(s10.partial("x"), s2 * s4 + x * (-s2 * s2 * s4 - s2 * s4 * s4));
141 auto s11 = math::invert(-3 * x - y);
142 auto s12 = math::invert(z);
143 auto s13 = x * s11 * s4 + x * y * z * s2 * s2 * s2 * s12;
144 BOOST_CHECK_EQUAL(s13.partial("x"), s11 * s4 + x * (3 * s11 * s11 * s4 - s11 * s4 * s4)
145 + y * z * s2 * s2 * s2 * s12
146 + x * y * z * (-3 * s2 * s2 * s2 * s2 * s12));
147 BOOST_CHECK_EQUAL(math::partial(s13, "x"), s11 * s4 + x * (3 * s11 * s11 * s4 - s11 * s4 * s4)
148 + y * z * s2 * s2 * s2 * s12
149 + x * y * z * (-3 * s2 * s2 * s2 * s2 * s12));
150 auto s15 = x * s11 * s4 + x * y * z * s2 * s2 * s2 * s12 + s4 * s12;
151 BOOST_CHECK_EQUAL(s15.partial("x"), s11 * s4 + x * (3 * s11 * s11 * s4 - s11 * s4 * s4)
152 + y * z * s2 * s2 * s2 * s12
153 + x * y * z * (-3 * s2 * s2 * s2 * s2 * s12) - s4 * s4 * s12);
154 // Overflow in an exponent.
155 overflow_check<T>();
156 auto s16 = math::invert(x - 4 * y);
157 auto s17 = s2 * s2 * s2 * s2 * s2 * s16 * s16 * s16 * s12;
158 BOOST_CHECK_EQUAL(s17.partial("x"), -5 * s2 * s2 * s2 * s2 * s2 * s2 * s16 * s16 * s16 * s12
159 - 3 * s2 * s2 * s2 * s2 * s2 * s16 * s16 * s16 * s16 * s12);
160 // Excercise the chain rule.
161 auto s18 = x * x * 3 / 4_q * y * z * z;
162 auto s19 = -y * y * x * z * z;
163 auto s20 = y * x * x * 4;
164 auto s21 = s18 * s17 + s19 * s2 * s11 * s12 + s20 * s16 * s2 * s3;
165 BOOST_CHECK_EQUAL(s21.partial("x"), s18.partial("x") * s17 + s18 * s17.partial("x")
166 + s19.partial("x") * s2 * s11 * s12
167 + s19 * (s2 * s11 * s12).partial("x") + s20.partial("x") * s16 * s2 * s3
168 + s20 * (s16 * s2 * s3).partial("x"));
169 BOOST_CHECK_EQUAL(s21.partial("y"), s18.partial("y") * s17 + s18 * s17.partial("y")
170 + s19.partial("y") * s2 * s11 * s12
171 + s19 * (s2 * s11 * s12).partial("y") + s20.partial("y") * s16 * s2 * s3
172 + s20 * (s16 * s2 * s3).partial("y"));
173 BOOST_CHECK_EQUAL(s21.partial("z"), s18.partial("z") * s17 + s18 * s17.partial("z")
174 + s19.partial("z") * s2 * s11 * s12
175 + s19 * (s2 * s11 * s12).partial("z") + s20.partial("z") * s16 * s2 * s3
176 + s20 * (s16 * s2 * s3).partial("z"));
177 BOOST_CHECK_EQUAL(s21.partial("v"), 0);
178 BOOST_CHECK_EQUAL(s_type{1}.partial("x"), 0);
179 }
180 template <typename T, typename U = T, typename std::enable_if<std::is_integral<U>::value, int>::type = 0>
overflow_checkpartial_tester181 static void overflow_check()
182 {
183 using p_type = polynomial<rational, monomial<int>>;
184 using s_type = divisor_series<p_type, divisor<T>>;
185 s_type s14;
186 symbol_set ss;
187 ss.add("x");
188 s14.set_symbol_set(ss);
189 typename s_type::term_type::key_type k0;
190 std::vector<T> vs;
191 vs.push_back(T(1));
192 T expo = std::numeric_limits<T>::max();
193 k0.insert(vs.begin(), vs.end(), expo);
194 s14.insert(typename s_type::term_type(typename s_type::term_type::cf_type(1), k0));
195 BOOST_CHECK_THROW(s14.partial("x"), std::overflow_error);
196 // Skip this overflow test if T is short, as short * short will become int * int and it will
197 // not overflow.
198 if (std::is_same<T, short>::value) {
199 return;
200 }
201 s_type s15;
202 ss.add("y");
203 s15.set_symbol_set(ss);
204 vs[0] = static_cast<T>(std::numeric_limits<T>::max() / T(4));
205 vs.push_back(T(1));
206 expo = static_cast<T>(std::numeric_limits<T>::max() - 1);
207 typename s_type::term_type::key_type k1;
208 k1.insert(vs.begin(), vs.end(), expo);
209 s15.insert(typename s_type::term_type(typename s_type::term_type::cf_type(1), k1));
210 BOOST_CHECK_THROW(s15.partial("x"), std::overflow_error);
211 }
212 template <typename T, typename U = T, typename std::enable_if<!std::is_integral<U>::value, int>::type = 0>
overflow_checkpartial_tester213 static void overflow_check()
214 {
215 }
216 };
217
BOOST_AUTO_TEST_CASE(divisor_series_partial_test)218 BOOST_AUTO_TEST_CASE(divisor_series_partial_test)
219 {
220 // A couple of general tests to start.
221 using p_type = polynomial<rational, monomial<int>>;
222 using s_type = divisor_series<p_type, divisor<short>>;
223 {
224 BOOST_CHECK_EQUAL(s_type{}.partial("x"), 0);
225 s_type s0{3};
226 BOOST_CHECK_EQUAL(s0.partial("x"), 0);
227 s_type x{"x"};
228 BOOST_CHECK_EQUAL((x * 3).partial("x"), 3);
229 BOOST_CHECK_EQUAL((x * 3).partial("y"), 0);
230 // Define an EPS.
231 using ps_type = poisson_series<s_type>;
232 ps_type a{"a"}, b{"b"}, c{"c"};
233 auto p1 = 3 * a * b * math::cos(3 * c);
234 BOOST_CHECK_EQUAL(boost::lexical_cast<std::string>(p1.t_integrate()), "a*b*1/[(\\nu_{c})]*sin(3*c)");
235 BOOST_CHECK_EQUAL(boost::lexical_cast<std::string>(p1.t_integrate().partial("a")), "b*1/[(\\nu_{c})]*sin(3*c)");
236 BOOST_CHECK_EQUAL(boost::lexical_cast<std::string>(p1.t_integrate().partial("b")), "a*1/[(\\nu_{c})]*sin(3*c)");
237 BOOST_CHECK_EQUAL(boost::lexical_cast<std::string>(p1.t_integrate().partial("c")),
238 "3*a*b*1/[(\\nu_{c})]*cos(3*c)");
239 BOOST_CHECK_EQUAL(boost::lexical_cast<std::string>(p1.t_integrate().partial("\\nu_{c}")),
240 "-a*b*1/[(\\nu_{c})**2]*sin(3*c)");
241 }
242 // Test with various exponent types.
243 boost::mpl::for_each<expo_types>(partial_tester());
244 // Test custom derivatives.
245 s_type x{"x"}, y{"y"};
246 s_type::register_custom_derivative(
247 "x", [x](const s_type &s) -> s_type { return s.partial("x") + math::partial(s, "y") * 2 * x; });
248 BOOST_CHECK_EQUAL(math::partial(math::invert(x + y), "x"), (-1 - 2 * x) * math::invert(x + y).pow(2));
249 s_type::register_custom_derivative(
250 "x", [y](const s_type &s) -> s_type { return s.partial("x") + math::partial(s, "y") * math::invert(y) / 2; });
251 BOOST_CHECK_EQUAL(math::partial(math::invert(x + 2 * y), "x"),
252 (-1 - 1 * y.invert()) * math::invert(x + 2 * y).pow(2));
253 s_type::register_custom_derivative(
254 "x", [y](const s_type &s) -> s_type { return s.partial("x") + math::partial(s, "y") * math::invert(y) / 2; });
255 BOOST_CHECK_EQUAL(math::partial(math::invert(x + y), "x"),
256 -math::invert(x + y).pow(2) - 1 / 2_q * math::invert(x + y).pow(2) * math::invert(y));
257 // Implicit variable dependency both in the poly and in the divisor.
258 s_type::register_custom_derivative(
259 "x", [x](const s_type &s) -> s_type { return s.partial("x") + math::partial(s, "y") * 2 * x; });
260 BOOST_CHECK_EQUAL(math::partial(y * math::invert(x + y), "x"),
261 2 * x * math::invert(x + y) - y * (2 * x + 1) * math::invert(x + y).pow(2));
262 }
263
BOOST_AUTO_TEST_CASE(divisor_series_integrate_test)264 BOOST_AUTO_TEST_CASE(divisor_series_integrate_test)
265 {
266 using s_type = divisor_series<polynomial<rational, monomial<short>>, divisor<short>>;
267 s_type x{"x"}, y{"y"}, z{"z"};
268 BOOST_CHECK((is_integrable<s_type>::value));
269 // A few cases with the variables only in the polynomial part.
270 BOOST_CHECK_EQUAL(x.integrate("x"), x * x / 2);
271 BOOST_CHECK_EQUAL(math::integrate(x, "x"), x * x / 2);
272 BOOST_CHECK((std::is_same<s_type, decltype(math::integrate(x, "x"))>::value));
273 BOOST_CHECK_EQUAL(math::integrate(x, "y"), x * y);
274 BOOST_CHECK_EQUAL(math::integrate(x + y, "x"), x * y + x * x / 2);
275 BOOST_CHECK_EQUAL(math::integrate(x + y, "y"), x * y + y * y / 2);
276 BOOST_CHECK_EQUAL(math::integrate(s_type{1}, "y"), y);
277 BOOST_CHECK_EQUAL(math::integrate(s_type{1}, "x"), x);
278 BOOST_CHECK_EQUAL(math::integrate(s_type{0}, "x"), 0);
279 // Put variables in the divisors as well.
280 BOOST_CHECK_EQUAL(math::integrate(x + y.invert(), "x"), x * x / 2 + x * y.invert());
281 BOOST_CHECK_THROW(math::integrate(x + y.invert() + x.invert(), "x"), std::invalid_argument);
282 BOOST_CHECK_EQUAL(math::integrate(x + y.invert() + x.invert() - x.invert(), "x"), x * x / 2 + x * y.invert());
283 }
284
BOOST_AUTO_TEST_CASE(divisor_series_invert_test)285 BOOST_AUTO_TEST_CASE(divisor_series_invert_test)
286 {
287 using s_type0 = divisor_series<int, divisor<short>>;
288 BOOST_CHECK_EQUAL(math::invert(s_type0{2}), 0);
289 using s_type1 = divisor_series<rational, divisor<short>>;
290 BOOST_CHECK_EQUAL(math::invert(s_type1{2}), 1 / 2_q);
291 BOOST_CHECK_EQUAL(math::invert(s_type1{2 / 3_q}), 3 / 2_q);
292 {
293 using s_type = divisor_series<polynomial<rational, monomial<short>>, divisor<short>>;
294 s_type x{"x"}, y{"y"}, z{"z"}, null;
295 BOOST_CHECK(is_invertible<s_type>::value);
296 BOOST_CHECK((std::is_same<s_type, decltype(math::invert(s_type{}))>::value));
297 BOOST_CHECK_EQUAL(boost::lexical_cast<std::string>(math::invert(x)), "1/[(x)]");
298 BOOST_CHECK_EQUAL(math::invert(s_type{2}), 1 / 2_q);
299 BOOST_CHECK_EQUAL(boost::lexical_cast<std::string>(math::pow(x, -1)), "x**-1");
300 BOOST_CHECK_THROW(math::invert(null), zero_division_error);
301 BOOST_CHECK((std::is_same<decltype(x.invert()), s_type>::value));
302 BOOST_CHECK((std::is_same<decltype(math::invert(x)), s_type>::value));
303 BOOST_CHECK_EQUAL(boost::lexical_cast<std::string>(math::invert(x - y)), "1/[(x-y)]");
304 BOOST_CHECK_EQUAL(boost::lexical_cast<std::string>(math::invert(2 * x - 4 * y)), "1/2*1/[(x-2*y)]");
305 BOOST_CHECK_EQUAL(boost::lexical_cast<std::string>(math::invert(-2 * x + 4 * y)), "-1/2*1/[(x-2*y)]");
306 BOOST_CHECK_EQUAL(boost::lexical_cast<std::string>(math::invert(x + y + z)), "1/[(x+y+z)]");
307 BOOST_CHECK_EQUAL(boost::lexical_cast<std::string>(math::invert(x + y + z - z)), "1/[(x+y)]");
308 BOOST_CHECK_THROW(math::invert(x - 1), std::invalid_argument);
309 BOOST_CHECK_THROW(math::invert(x - y / 2), std::invalid_argument);
310 BOOST_CHECK_THROW(math::invert(x - x), zero_division_error);
311 // Out of bounds for short.
312 BOOST_CHECK_THROW(math::invert((std::numeric_limits<short>::max() + 1_q) * x + y), std::invalid_argument);
313 // Check, if appropriate, construction from outside the bounds defined in divisor.
314 if (detail::safe_abs_sint<short>::value < std::numeric_limits<short>::max()) {
315 BOOST_CHECK_THROW(math::invert((detail::safe_abs_sint<short>::value + 1_q) * x + y), std::invalid_argument);
316 }
317 if (-detail::safe_abs_sint<short>::value > std::numeric_limits<short>::min()) {
318 BOOST_CHECK_THROW(math::invert((-detail::safe_abs_sint<short>::value - 1_q) * x + y),
319 std::invalid_argument);
320 }
321 }
322 {
323 using s_type = divisor_series<polynomial<rational, k_monomial>, divisor<short>>;
324 s_type x{"x"}, y{"y"}, z{"z"}, null;
325 BOOST_CHECK(is_invertible<s_type>::value);
326 BOOST_CHECK((std::is_same<s_type, decltype(math::invert(s_type{}))>::value));
327 BOOST_CHECK_EQUAL(boost::lexical_cast<std::string>(math::invert(x)), "1/[(x)]");
328 BOOST_CHECK_EQUAL(math::invert(s_type{2}), 1 / 2_q);
329 BOOST_CHECK_EQUAL(boost::lexical_cast<std::string>(math::pow(x, -1)), "x**-1");
330 BOOST_CHECK_THROW(math::invert(null), zero_division_error);
331 BOOST_CHECK((std::is_same<decltype(x.invert()), s_type>::value));
332 BOOST_CHECK((std::is_same<decltype(math::invert(x)), s_type>::value));
333 BOOST_CHECK_EQUAL(boost::lexical_cast<std::string>(math::invert(x - y)), "1/[(x-y)]");
334 BOOST_CHECK_EQUAL(boost::lexical_cast<std::string>(math::invert(2 * x - 4 * y)), "1/2*1/[(x-2*y)]");
335 BOOST_CHECK_EQUAL(boost::lexical_cast<std::string>(math::invert(-2 * x + 4 * y)), "-1/2*1/[(x-2*y)]");
336 BOOST_CHECK_EQUAL(boost::lexical_cast<std::string>(math::invert(x + y + z)), "1/[(x+y+z)]");
337 BOOST_CHECK_EQUAL(boost::lexical_cast<std::string>(math::invert(x + y + z - z)), "1/[(x+y)]");
338 BOOST_CHECK_THROW(math::invert(x - 1), std::invalid_argument);
339 BOOST_CHECK_THROW(math::invert(x - y / 2), std::invalid_argument);
340 BOOST_CHECK_THROW(math::invert(x - x), zero_division_error);
341 // Out of bounds for short.
342 BOOST_CHECK_THROW(math::invert((std::numeric_limits<short>::max() + 1_q) * x + y), std::invalid_argument);
343 // Check, if appropriate, construction from outside the bounds defined in divisor.
344 if (detail::safe_abs_sint<short>::value < std::numeric_limits<short>::max()) {
345 BOOST_CHECK_THROW(math::invert((detail::safe_abs_sint<short>::value + 1_q) * x + y), std::invalid_argument);
346 }
347 if (-detail::safe_abs_sint<short>::value > std::numeric_limits<short>::min()) {
348 BOOST_CHECK_THROW(math::invert((-detail::safe_abs_sint<short>::value - 1_q) * x + y),
349 std::invalid_argument);
350 }
351 }
352 {
353 using s_type = divisor_series<polynomial<rational, monomial<rational>>, divisor<short>>;
354 s_type x{"x"}, y{"y"}, z{"z"}, null;
355 BOOST_CHECK(is_invertible<s_type>::value);
356 BOOST_CHECK((std::is_same<s_type, decltype(math::invert(s_type{}))>::value));
357 BOOST_CHECK_EQUAL(boost::lexical_cast<std::string>(math::invert(x)), "1/[(x)]");
358 BOOST_CHECK_EQUAL(math::invert(s_type{2}), 1 / 2_q);
359 BOOST_CHECK_EQUAL(boost::lexical_cast<std::string>(math::pow(x, -1)), "x**-1");
360 BOOST_CHECK_THROW(math::invert(null), zero_division_error);
361 BOOST_CHECK((std::is_same<decltype(x.invert()), s_type>::value));
362 BOOST_CHECK((std::is_same<decltype(math::invert(x)), s_type>::value));
363 BOOST_CHECK_EQUAL(boost::lexical_cast<std::string>(math::invert(x - y)), "1/[(x-y)]");
364 BOOST_CHECK_EQUAL(boost::lexical_cast<std::string>(math::invert(2 * x - 4 * y)), "1/2*1/[(x-2*y)]");
365 BOOST_CHECK_EQUAL(boost::lexical_cast<std::string>(math::invert(-2 * x + 4 * y)), "-1/2*1/[(x-2*y)]");
366 BOOST_CHECK_EQUAL(boost::lexical_cast<std::string>(math::invert(x + y + z)), "1/[(x+y+z)]");
367 BOOST_CHECK_EQUAL(boost::lexical_cast<std::string>(math::invert(x + y + z - z)), "1/[(x+y)]");
368 BOOST_CHECK_THROW(math::invert(x - 1), std::invalid_argument);
369 BOOST_CHECK_THROW(math::invert(x - y / 2), std::invalid_argument);
370 BOOST_CHECK_THROW(math::invert(x - x), zero_division_error);
371 // Out of bounds for short.
372 BOOST_CHECK_THROW(math::invert((std::numeric_limits<short>::max() + 1_q) * x + y), std::invalid_argument);
373 // Check, if appropriate, construction from outside the bounds defined in divisor.
374 if (detail::safe_abs_sint<short>::value < std::numeric_limits<short>::max()) {
375 BOOST_CHECK_THROW(math::invert((detail::safe_abs_sint<short>::value + 1_q) * x + y), std::invalid_argument);
376 }
377 if (-detail::safe_abs_sint<short>::value > std::numeric_limits<short>::min()) {
378 BOOST_CHECK_THROW(math::invert((-detail::safe_abs_sint<short>::value - 1_q) * x + y),
379 std::invalid_argument);
380 }
381 }
382 {
383 // Try with something else between poly and divisor.
384 using s_type = divisor_series<poisson_series<polynomial<rational, monomial<short>>>, divisor<short>>;
385 BOOST_CHECK(is_invertible<s_type>::value);
386 BOOST_CHECK((std::is_same<s_type, decltype(math::invert(s_type{}))>::value));
387 s_type x{"x"}, y{"y"}, null;
388 BOOST_CHECK_EQUAL(boost::lexical_cast<std::string>(math::pow(2 * x, -1)), "1/2*x**-1");
389 BOOST_CHECK_EQUAL(boost::lexical_cast<std::string>(math::invert(2 * x)), "1/2*1/[(x)]");
390 BOOST_CHECK_THROW(math::invert(math::cos(2 * x)), std::invalid_argument);
391 BOOST_CHECK_THROW(math::pow(x + y, -1), std::invalid_argument);
392 BOOST_CHECK_EQUAL(boost::lexical_cast<std::string>(math::invert(-2 * x + 4 * y)), "-1/2*1/[(x-2*y)]");
393 BOOST_CHECK_THROW(math::invert(null), zero_division_error);
394 BOOST_CHECK_THROW(math::pow(null, -1), zero_division_error);
395 }
396 }
397
BOOST_AUTO_TEST_CASE(divisor_series_rational_multiplication_test)398 BOOST_AUTO_TEST_CASE(divisor_series_rational_multiplication_test)
399 {
400 // Test that we handle correctly rational coefficients wrt the lcm computation in the multiplier.
401 using s_type = divisor_series<rational, divisor<short>>;
402 s_type s1{1 / 2_q}, s2{2 / 3_q};
403 BOOST_CHECK_EQUAL(s1 * s2, 1 / 3_q);
404 }
405