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