monomial.hpp

/* Copyright 2009-2017 Francesco Biscani (bluescarni@gmail.com)

This file is part of the Piranha library.

The Piranha library is free software; you can redistribute it and/or modify
it under the terms of either:

  * the GNU Lesser General Public License as published by the Free
    Software Foundation; either version 3 of the License, or (at your
    option) any later version.

or

  * the GNU General Public License as published by the Free Software
    Foundation; either version 3 of the License, or (at your option) any
    later version.

or both in parallel, as here.

The Piranha library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
for more details.

You should have received copies of the GNU General Public License and the
GNU Lesser General Public License along with the Piranha library.  If not,
see https://www.gnu.org/licenses/. */

#ifndef PIRANHA_MONOMIAL_HPP
#define PIRANHA_MONOMIAL_HPP

#include <algorithm>
#include <array>
#include <cstddef>
#include <cstdint>
#include <functional>
#include <initializer_list>
#include <iostream>
#include <iterator>
#include <limits>
#include <sstream>
#include <stdexcept>
#include <string>
#include <tuple>
#include <type_traits>
#include <utility>
#include <vector>

#include <piranha/array_key.hpp>
#include <piranha/config.hpp>
#include <piranha/detail/cf_mult_impl.hpp>
#include <piranha/detail/monomial_common.hpp>
#include <piranha/detail/prepare_for_print.hpp>
#include <piranha/detail/safe_integral_adder.hpp>
#include <piranha/exceptions.hpp>
#include <piranha/forwarding.hpp>
#include <piranha/is_cf.hpp>
#include <piranha/is_key.hpp>
#include <piranha/math.hpp>
#include <piranha/mp_integer.hpp>
#include <piranha/mp_rational.hpp>
#include <piranha/pow.hpp>
#include <piranha/s11n.hpp>
#include <piranha/safe_cast.hpp>
#include <piranha/symbol_utils.hpp>
#include <piranha/term.hpp>
#include <piranha/type_traits.hpp>

namespace piranha
{


template <typename T, typename S = std::integral_constant<std::size_t, 0u>>
class monomial : public array_key<T, monomial<T, S>, S>
{
    PIRANHA_TT_CHECK(has_is_unitary, T);
    PIRANHA_TT_CHECK(is_ostreamable, T);
    PIRANHA_TT_CHECK(has_negate, T);
    PIRANHA_TT_CHECK(std::is_copy_assignable, T);
    using base = array_key<T, monomial<T, S>, S>;
    // Less-than operator.
    template <typename U>
    using comparison_enabler = typename std::enable_if<is_less_than_comparable<U>::value, int>::type;

public:
    static const std::size_t multiply_arity = 1u;
    monomial() = default;
    monomial(const monomial &) = default;
    monomial(monomial &&) = default;

private:
    // Enabler for ctor from init list.
    template <typename U>
    using init_list_enabler = enable_if_t<std::is_constructible<base, std::initializer_list<U>>::value, int>;

public:

    template <typename U, init_list_enabler<U> = 0>
    explicit monomial(std::initializer_list<U> list) : base(list)
    {
    }

private:
    // Enabler for ctor from range.
    template <typename Iterator>
    using it_ctor_enabler
        = enable_if_t<conjunction<is_input_iterator<Iterator>,
                                  has_safe_cast<T, typename std::iterator_traits<Iterator>::value_type>>::value,
                      int>;

public:

    template <typename Iterator, it_ctor_enabler<Iterator> = 0>
    explicit monomial(Iterator begin, Iterator end)
    {
        for (; begin != end; ++begin) {
            this->push_back(safe_cast<T>(*begin));
        }
    }

    template <typename Iterator, it_ctor_enabler<Iterator> = 0>
    explicit monomial(Iterator begin, Iterator end, const symbol_fset &s) : monomial(begin, end)
    {
        if (unlikely(this->size() != s.size())) {
            piranha_throw(std::invalid_argument, "the monomial constructor from range and symbol set "
                                                 "yielded an invalid monomial: the final size is "
                                                     + std::to_string(this->size())
                                                     + ", while the size of the symbol set is "
                                                     + std::to_string(s.size()));
        }
    }
    PIRANHA_FORWARDING_CTOR(monomial, base)
    ~monomial()
    {
        PIRANHA_TT_CHECK(is_key, monomial);
    }

    monomial &operator=(const monomial &other) = default;

    monomial &operator=(monomial &&other) = default;

    bool is_compatible(const symbol_fset &args) const noexcept
    {
        return this->size() == args.size();
    }

    bool is_zero(const symbol_fset &) const noexcept
    {
        return false;
    }

    bool is_unitary(const symbol_fset &args) const
    {
        const auto sbe = this->size_begin_end();
        if (unlikely(args.size() != std::get<0>(sbe))) {
            piranha_throw(std::invalid_argument,
                          "invalid sizes in the invocation of is_unitary() for a monomial: the monomial has a size of "
                              + std::to_string(std::get<0>(sbe)) + ", while the reference symbol set has a size of "
                              + std::to_string(args.size()));
        }
        return std::all_of(std::get<1>(sbe), std::get<2>(sbe), [](const T &element) { return math::is_zero(element); });
    }

private:
    // Machinery to determine the degree type.
    template <typename U>
    using degree_type = enable_if_t<conjunction<std::is_constructible<add_t<U, U>, int>,
                                                is_addable_in_place<add_t<U, U>, U>, is_returnable<add_t<U, U>>>::value,
                                    add_t<U, U>>;
    // Helpers to add exponents in the degree computation.
    template <typename U, enable_if_t<std::is_integral<U>::value, int> = 0>
    static void expo_add(degree_type<U> &retval, const U &n)
    {
        detail::safe_integral_adder(retval, static_cast<degree_type<U>>(n));
    }
    template <typename U, enable_if_t<!std::is_integral<U>::value, int> = 0>
    static void expo_add(degree_type<U> &retval, const U &x)
    {
        retval += x;
    }

public:

    template <typename U = T>
    degree_type<U> degree(const symbol_fset &args) const
    {
        auto sbe = this->size_begin_end();
        if (unlikely(args.size() != std::get<0>(sbe))) {
            piranha_throw(
                std::invalid_argument,
                "invalid symbol set for the computation of the degree of a monomial: the size of the symbol set ("
                    + std::to_string(args.size()) + ") differs from the size of the monomial ("
                    + std::to_string(std::get<0>(sbe)) + ")");
        }
        degree_type<U> retval(0);
        for (; std::get<1>(sbe) != std::get<2>(sbe); ++std::get<1>(sbe)) {
            expo_add(retval, *std::get<1>(sbe));
        }
        return retval;
    }

    template <typename U = T>
    degree_type<U> degree(const symbol_idx_fset &p, const symbol_fset &args) const
    {
        auto sbe = this->size_begin_end();
        if (unlikely(args.size() != std::get<0>(sbe))) {
            piranha_throw(std::invalid_argument, "invalid symbol set for the computation of the partial degree of a "
                                                 "monomial: the size of the symbol set ("
                                                     + std::to_string(args.size())
                                                     + ") differs from the size of the monomial ("
                                                     + std::to_string(std::get<0>(sbe)) + ")");
        }
        if (unlikely(p.size() && *p.rbegin() >= args.size())) {
            piranha_throw(std::invalid_argument, "the largest value in the positions set for the computation of the "
                                                 "partial degree of a monomial is "
                                                     + std::to_string(*p.rbegin())
                                                     + ", but the monomial has a size of only "
                                                     + std::to_string(args.size()));
        }
        degree_type<U> retval(0);
        for (const auto &i : p) {
            expo_add(retval, std::get<1>(sbe)[i]);
        }
        return retval;
    }

    template <typename U = T>
    degree_type<U> ldegree(const symbol_fset &args) const
    {
        return degree(args);
    }

    template <typename U = T>
    degree_type<U> ldegree(const symbol_idx_fset &p, const symbol_fset &args) const
    {
        return degree(p, args);
    }

    std::pair<bool, symbol_idx> is_linear(const symbol_fset &args) const
    {
        auto sbe = this->size_begin_end();
        if (unlikely(args.size() != std::get<0>(sbe))) {
            piranha_throw(std::invalid_argument, "invalid symbol set for the identification of a linear "
                                                 "monomial: the size of the symbol set ("
                                                     + std::to_string(args.size())
                                                     + ") differs from the size of the monomial ("
                                                     + std::to_string(std::get<0>(sbe)) + ")");
        }
        using size_type = typename base::size_type;
        size_type n_linear = 0u, candidate = 0u;
        for (size_type i = 0u; i < std::get<0>(sbe); ++i, ++std::get<1>(sbe)) {
            // NOTE: is_zero()'s availability is guaranteed by array_key's reqs,
            // is_unitary() is required by the monomial reqs.
            if (math::is_zero(*std::get<1>(sbe))) {
                continue;
            }
            if (!math::is_unitary(*std::get<1>(sbe))) {
                return std::make_pair(false, symbol_idx{0});
            }
            candidate = i;
            ++n_linear;
        }
        piranha_assert(std::get<1>(sbe) == std::get<2>(sbe));
        if (n_linear != 1u) {
            return std::make_pair(false, symbol_idx{0});
        }
        return std::make_pair(true, symbol_idx{candidate});
    }

private:
    // Enabler for pow.
    template <typename U>
    using pow_enabler = monomial_pow_enabler<T, U>;

public:

    template <typename U, pow_enabler<U> = 0>
    monomial pow(const U &x, const symbol_fset &args) const
    {
        auto sbe = this->size_begin_end();
        if (unlikely(args.size() != std::get<0>(sbe))) {
            piranha_throw(std::invalid_argument, "invalid symbol set for the exponentiation of a "
                                                 "monomial: the size of the symbol set ("
                                                     + std::to_string(args.size())
                                                     + ") differs from the size of the monomial ("
                                                     + std::to_string(std::get<0>(sbe)) + ")");
        }
        // Init with zeroes.
        monomial retval(args);
        for (decltype(retval.size()) i = 0u; i < std::get<0>(sbe); ++i, ++std::get<1>(sbe)) {
            monomial_pow_mult_exp(retval[i], *std::get<1>(sbe), x, monomial_pow_dispatcher<T, U>{});
        }
        return retval;
    }

private:
    // Partial utils.
    // Detect decrement operator.
    template <typename U>
    using dec_t = decltype(--std::declval<U &>());
    // Main dispatcher.
    template <typename U>
    using monomial_partial_dispatcher = disjunction_idx<
        // Case 0: U is integral.
        std::is_integral<U>,
        // Case 1: U supports the decrement operator.
        is_detected<dec_t, U>>;
    // In-place decrement by one, checked for integral types.
    template <typename U>
    static void ip_dec(U &x, const std::integral_constant<std::size_t, 0u> &)
    {
        // NOTE: maybe replace with the safe integral subber, eventually.
        if (unlikely(x == std::numeric_limits<U>::min())) {
            piranha_throw(std::overflow_error, "negative overflow error in the calculation of the "
                                               "partial derivative of a monomial");
        }
        x = static_cast<U>(x - U(1));
    }
    template <typename U>
    static void ip_dec(U &x, const std::integral_constant<std::size_t, 1u> &)
    {
        --x;
    }
    // The enabler.
    template <typename U>
    using partial_enabler = enable_if_t<(monomial_partial_dispatcher<U>::value < 2u), int>;

public:

    template <typename U = T, partial_enabler<U> = 0>
    std::pair<T, monomial> partial(const symbol_idx &p, const symbol_fset &args) const
    {
        auto sbe = this->size_begin_end();
        if (unlikely(args.size() != std::get<0>(sbe))) {
            piranha_throw(std::invalid_argument,
                          "invalid symbol set for the computation of the partial derivative of a "
                          "monomial: the size of the symbol set ("
                              + std::to_string(args.size()) + ") differs from the size of the monomial ("
                              + std::to_string(std::get<0>(sbe)) + ")");
        }
        if (p >= std::get<0>(sbe) || math::is_zero(std::get<1>(sbe)[p])) {
            // Derivative wrt a variable not in the monomial: position is outside the bounds, or it refers to a
            // variable with zero exponent.
            return std::make_pair(T(0), monomial{args});
        }
        // Copy the original exponent.
        T expo(std::get<1>(sbe)[p]);
        // Copy the original monomial.
        monomial m(*this);
        // Decrement the affected exponent in m.
        ip_dec(m[static_cast<typename base::size_type>(p)], monomial_partial_dispatcher<U>{});
        // Return the result.
        return std::make_pair(std::move(expo), std::move(m));
    }

private:
    // Integrate utils.
    // Detect increment operator.
    template <typename U>
    using inc_t = decltype(++std::declval<U &>());
    // Main dispatcher.
    template <typename U>
    using monomial_int_dispatcher = disjunction_idx<
        // Case 0: U is integral.
        std::is_integral<U>,
        // Case 1: U supports the increment operator.
        is_detected<inc_t, U>>;
    // In-place increment by one, checked for integral types.
    template <typename U>
    static void ip_inc(U &x, const std::integral_constant<std::size_t, 0u> &)
    {
        // NOTE: maybe replace with the safe integral adder, eventually.
        if (unlikely(x == std::numeric_limits<U>::max())) {
            piranha_throw(std::overflow_error, "positive overflow error in the calculation of the "
                                               "antiderivative of a monomial");
        }
        x = static_cast<U>(x + U(1));
    }
    template <typename U>
    static void ip_inc(U &x, const std::integral_constant<std::size_t, 1u> &)
    {
        ++x;
    }
    // The enabler.
    template <typename U>
    using integrate_enabler = enable_if_t<(monomial_int_dispatcher<U>::value < 2u), int>;

public:

    template <typename U = T, integrate_enabler<U> = 0>
    std::pair<T, monomial> integrate(const std::string &s, const symbol_fset &args) const
    {
        auto sbe = this->size_begin_end();
        if (unlikely(args.size() != std::get<0>(sbe))) {
            piranha_throw(std::invalid_argument, "invalid symbol set for the computation of the antiderivative of a "
                                                 "monomial: the size of the symbol set ("
                                                     + std::to_string(args.size())
                                                     + ") differs from the size of the monomial ("
                                                     + std::to_string(std::get<0>(sbe)) + ")");
        }
        monomial retval;
        T expo(0);
        const PIRANHA_MAYBE_TLS T one(1);
        auto it_args = args.begin();
        for (decltype(retval.size()) i = 0; std::get<1>(sbe) != std::get<2>(sbe); ++i, ++std::get<1>(sbe), ++it_args) {
            const auto &cur_sym = *it_args;
            if (math::is_zero(expo) && s < cur_sym) {
                // If we went past the position of s in args and still we
                // have not performed the integration, it means that we need to add
                // a new exponent.
                retval.push_back(one);
                expo = one;
            }
            retval.push_back(*std::get<1>(sbe));
            if (cur_sym == s) {
                // NOTE: here using i is safe: if retval gained an extra exponent in the condition above,
                // we are never going to land here as cur_sym is at this point never going to be s.
                auto &ref = retval[i];
                // Do the addition and check for zero later, to detect -1 expo.
                ip_inc(ref, monomial_int_dispatcher<U>{});
                if (unlikely(math::is_zero(ref))) {
                    piranha_throw(std::invalid_argument,
                                  "unable to perform monomial integration: a negative "
                                  "unitary exponent was encountered in correspondence of the variable '"
                                      + cur_sym + "'");
                }
                expo = ref;
            }
        }
        // If expo is still zero, it means we need to add a new exponent at the end.
        if (math::is_zero(expo)) {
            retval.push_back(one);
            expo = one;
        }
        return std::make_pair(std::move(expo), std::move(retval));
    }

private:
    // Let's hard code the custom behaviour for rational exponents for the moment.
    // We can offer a customization point in the future.
    template <typename U, enable_if_t<is_mp_rational<U>::value, int> = 0>
    static void print_exponent(std::ostream &os, const U &e)
    {
        if (math::is_unitary(e.den())) {
            os << e;
        } else {
            os << '(' << e << ')';
        }
    }
    template <typename U, enable_if_t<!is_mp_rational<U>::value, int> = 0>
    static void print_exponent(std::ostream &os, const U &e)
    {
        os << detail::prepare_for_print(e);
    }

public:

    void print(std::ostream &os, const symbol_fset &args) const
    {
        auto sbe = this->size_begin_end();
        if (unlikely(args.size() != std::get<0>(sbe))) {
            piranha_throw(std::invalid_argument,
                          "cannot print monomial: the size of the symbol set (" + std::to_string(args.size())
                              + ") differs from the size of the monomial (" + std::to_string(std::get<0>(sbe)) + ")");
        }
        bool empty_output = true;
        auto it_args = args.begin();
        for (decltype(args.size()) i = 0; std::get<1>(sbe) != std::get<2>(sbe); ++i, ++std::get<1>(sbe), ++it_args) {
            if (!math::is_zero(*std::get<1>(sbe))) {
                // If we are going to print a symbol, and something has been printed before,
                // then we are going to place the multiplication sign.
                if (!empty_output) {
                    os << '*';
                }
                os << *it_args;
                empty_output = false;
                if (!math::is_unitary(*std::get<1>(sbe))) {
                    os << "**";
                    print_exponent(os, *std::get<1>(sbe));
                }
            }
        }
    }

    void print_tex(std::ostream &os, const symbol_fset &args) const
    {
        auto sbe = this->size_begin_end();
        if (unlikely(args.size() != std::get<0>(sbe))) {
            piranha_throw(std::invalid_argument, "cannot print monomial in TeX mode: the size of the symbol set ("
                                                     + std::to_string(args.size())
                                                     + ") differs from the size of the monomial ("
                                                     + std::to_string(std::get<0>(sbe)) + ")");
        }
        std::ostringstream oss_num, oss_den, *cur_oss;
        const PIRANHA_MAYBE_TLS T zero(0);
        T cur_value;
        auto it_args = args.begin();
        for (decltype(args.size()) i = 0; std::get<1>(sbe) != std::get<2>(sbe); ++i, ++std::get<1>(sbe), ++it_args) {
            cur_value = *std::get<1>(sbe);
            if (!math::is_zero(cur_value)) {
                // NOTE: use this weird form for the test because the presence of operator<()
                // is already guaranteed and thus we don't need additional requirements on T.
                // Maybe in the future we want to do it with a math::sign() function.
                if (zero < cur_value) {
                    cur_oss = std::addressof(oss_num);
                } else {
                    math::negate(cur_value);
                    cur_oss = std::addressof(oss_den);
                }
                *cur_oss << "{" << *it_args << "}";
                if (!math::is_unitary(cur_value)) {
                    *cur_oss << "^{" << detail::prepare_for_print(cur_value) << "}";
                }
            }
        }
        const std::string num_str = oss_num.str(), den_str = oss_den.str();
        if (!num_str.empty() && !den_str.empty()) {
            os << "\\frac{" << num_str << "}{" << den_str << "}";
        } else if (!num_str.empty() && den_str.empty()) {
            os << num_str;
        } else if (num_str.empty() && !den_str.empty()) {
            os << "\\frac{1}{" << den_str << "}";
        }
    }

private:
    // Eval type definition.
    template <typename U>
    using eval_type
        = enable_if_t<conjunction<is_multipliable_in_place<pow_t<U, T>>, std::is_constructible<pow_t<U, T>, int>,
                                  is_returnable<pow_t<U, T>>>::value,
                      pow_t<U, T>>;

public:

    template <typename U>
    eval_type<U> evaluate(const std::vector<U> &values, const symbol_fset &args) const
    {
        auto sbe = this->size_begin_end();
        if (unlikely(args.size() != std::get<0>(sbe))) {
            piranha_throw(std::invalid_argument,
                          "cannot evaluate monomial: the size of the symbol set (" + std::to_string(args.size())
                              + ") differs from the size of the monomial (" + std::to_string(std::get<0>(sbe)) + ")");
        }
        if (unlikely(values.size() != std::get<0>(sbe))) {
            piranha_throw(std::invalid_argument,
                          "cannot evaluate monomial: the size of the vector of values (" + std::to_string(values.size())
                              + ") differs from the size of the monomial (" + std::to_string(std::get<0>(sbe)) + ")");
        }
        if (args.size()) {
            eval_type<U> retval(math::pow(values[0], *std::get<1>(sbe)++));
            for (decltype(values.size()) i = 1; i < values.size(); ++i, ++std::get<1>(sbe)) {
                // NOTE: here maybe we could use mul3() and pow3() (to be implemented?).
                // NOTE: math::pow() for C++ integrals produces an integer result, no need
                // to worry about overflows.
                retval *= math::pow(values[i], *std::get<1>(sbe));
            }
            piranha_assert(std::get<1>(sbe) == std::get<2>(sbe));
            return retval;
        }
        return eval_type<U>(1);
    }

private:
    // Subs type is the same as eval type.
    template <typename U>
    using subs_type = eval_type<U>;

public:

    template <typename U>
    std::vector<std::pair<subs_type<U>, monomial>> subs(const symbol_idx_fmap<U> &smap, const symbol_fset &args) const
    {
        auto sbe = this->size_begin_end();
        if (unlikely(args.size() != std::get<0>(sbe))) {
            piranha_throw(std::invalid_argument,
                          "cannot perform substitution in a monomial: the size of the symbol set ("
                              + std::to_string(args.size()) + ") differs from the size of the monomial ("
                              + std::to_string(std::get<0>(sbe)) + ")");
        }
        if (unlikely(smap.size() && smap.rbegin()->first >= std::get<0>(sbe))) {
            // The last element of the substitution map must be a valid index.
            piranha_throw(std::invalid_argument,
                          "invalid argument(s) for substitution in a monomial: the last index of the substitution map ("
                              + std::to_string(smap.rbegin()->first) + ") must be smaller than the monomial's size ("
                              + std::to_string(std::get<0>(sbe)) + ")");
        }
        std::vector<std::pair<subs_type<U>, monomial>> retval;
        if (smap.size()) {
            // The substitution map contains something, proceed to the substitution.
            // Cache-init the zero constant.
            const PIRANHA_MAYBE_TLS T zero(0);
            // Init the subs return value from the exponentiation of the first value in the map.
            auto it = smap.begin();
            auto ret(math::pow(it->second, std::get<1>(sbe)[it->first]));
            // Init the monomial return value with a copy of this.
            auto mon_ret(*this);
            // Zero out the corresponding exponent.
            mon_ret[static_cast<decltype(mon_ret.size())>(it->first)] = zero;
            //  NOTE: move to the next element in the init statement of the for loop.
            for (++it; it != smap.end(); ++it) {
                ret *= math::pow(it->second, std::get<1>(sbe)[it->first]);
                mon_ret[static_cast<decltype(mon_ret.size())>(it->first)] = zero;
            }
            // NOTE: the is_returnable requirement ensures we can emplace back a pair
            // containing the subs type.
            retval.emplace_back(std::move(ret), std::move(mon_ret));
        } else {
            // Otherwise, the substitution yields 1 and the monomial is the original one.
            retval.emplace_back(subs_type<U>(1), *this);
        }
        return retval;
    }

private:
    // ipow subs utilities.
    // Dispatcher for the assignment of an exponent to an integer.
    template <typename U>
    using ipow_subs_d_assign_dispatcher = disjunction_idx<
        // Case 0: U is integer or a C++ integral.
        disjunction<std::is_integral<U>, std::is_same<integer, U>>,
        // Case 1: U supports safe cast to integer.
        has_safe_cast<integer, U>>;
    template <typename U>
    static void ipow_subs_d_assign(integer &d, const U &expo, const std::integral_constant<std::size_t, 0> &)
    {
        d = expo;
    }
    template <typename U>
    static void ipow_subs_d_assign(integer &d, const U &expo, const std::integral_constant<std::size_t, 1> &)
    {
        d = safe_cast<integer>(expo);
    }
    // Dispatcher for the assignment of an integer to an exponent.
    template <typename U>
    using ipow_subs_expo_assign_dispatcher = disjunction_idx<
        // Case 0: U is integer.
        std::is_same<integer, U>,
        // Case 1: integer supports safe cast to U.
        has_safe_cast<U, integer>>;
    template <typename U>
    static void ipow_subs_expo_assign(U &expo, const integer &r, const std::integral_constant<std::size_t, 0> &)
    {
        expo = r;
    }
    template <typename U>
    static void ipow_subs_expo_assign(U &expo, const integer &r, const std::integral_constant<std::size_t, 1> &)
    {
        expo = safe_cast<U>(r);
    }
    // Definition of the return type.
    template <typename U>
    using ips_type = decltype(math::pow(std::declval<const U &>(), std::declval<const integer &>()));
    // Final enabler.
    template <typename U>
    using ipow_subs_type
        = enable_if_t<conjunction<std::is_constructible<ips_type<U>, int>, is_returnable<ips_type<U>>>::value
                          && (ipow_subs_d_assign_dispatcher<T>::value < 2u)
                          && (ipow_subs_expo_assign_dispatcher<T>::value < 2u),
                      ips_type<U>>;

public:

    template <typename U>
    std::vector<std::pair<ipow_subs_type<U>, monomial>> ipow_subs(const symbol_idx &p, const integer &n, const U &x,
                                                                  const symbol_fset &args) const
    {
        if (unlikely(this->size() != args.size())) {
            piranha_throw(std::invalid_argument,
                          "cannot perform integral power substitution in a monomial: the size of the symbol set ("
                              + std::to_string(args.size()) + ") differs from the size of the monomial ("
                              + std::to_string(this->size()) + ")");
        }
        if (unlikely(!n.sgn())) {
            piranha_throw(std::invalid_argument,
                          "invalid integral power for ipow_subs() in a monomial: the power must be nonzero");
        }
        // Initialise the return value.
        std::vector<std::pair<ipow_subs_type<U>, monomial>> retval;
        auto mon(*this);
        if (p < args.size()) {
            auto &expo = mon[static_cast<decltype(mon.size())>(p)];
            PIRANHA_MAYBE_TLS integer q, r, d;
            // Assign expo to d, possibly safely converting it.
            ipow_subs_d_assign(d, expo, ipow_subs_d_assign_dispatcher<T>{});
            // NOTE: regarding the sign of r: tdiv_qr() sets the sign of r to the sign of q.
            // The only two cases we are interested in here are where d and n have the same sign
            // (otherwise q will have negative sign and we never enter the 'if' below). With
            // d and n positive, everything is straightforward (r's sign will be positive).
            // If d and n are both negative, r will have negative sign, and it will satisfy:
            // q*n + r == d (with d < 0 and d < q*n)
            // This is the result we want: r is the number of steps towards -inf that q*n
            // must take to reach d.
            tdiv_qr(q, r, d, n);
            if (q.sgn() > 0) {
                // Assign back the remainder r to expo, possibly with a safe
                // conversion involved.
                ipow_subs_expo_assign(expo, r, ipow_subs_expo_assign_dispatcher<T>{});
                retval.emplace_back(math::pow(x, q), std::move(mon));
                return retval;
            }
        }
        // Otherwise, the substitution yields 1 and the monomial is the original one.
        retval.emplace_back(ipow_subs_type<U>(1), std::move(mon));
        return retval;
    }

private:
    // Multiplication.
    template <typename Cf>
    using multiply_enabler = enable_if_t<conjunction<has_add3<T>, has_mul3<Cf>>::value, int>;

public:

    // NOTE: it should be ok to use this method (and the one below) with overlapping arguments, as this is allowed
    // in small_vector::add().
    template <typename Cf, multiply_enabler<Cf> = 0>
    static void multiply(std::array<term<Cf, monomial>, multiply_arity> &res, const term<Cf, monomial> &t1,
                         const term<Cf, monomial> &t2, const symbol_fset &args)
    {
        term<Cf, monomial> &t = res[0u];
        // NOTE: the check on the monomials' size is in vector_add().
        if (unlikely(t1.m_key.size() != args.size())) {
            piranha_throw(std::invalid_argument,
                          "cannot multiply terms with monomial keys: the size of the symbol set ("
                              + std::to_string(args.size()) + ") differs from the size of the first monomial ("
                              + std::to_string(t1.m_key.size()) + ")");
        }
        // Coefficient.
        cf_mult_impl(t.m_cf, t1.m_cf, t2.m_cf);
        // Now deal with the key.
        t1.m_key.vector_add(t.m_key, t2.m_key);
    }

    bool operator<(const monomial &other) const
    {
        const auto sbe1 = this->size_begin_end();
        const auto sbe2 = other.size_begin_end();
        if (unlikely(std::get<0u>(sbe1) != std::get<0u>(sbe2))) {
            piranha_throw(std::invalid_argument,
                          "mismatched sizes in a monomial comparison: the first monomial has a size of "
                              + std::to_string(std::get<0u>(sbe1)) + ", the second monomial has a size of "
                              + std::to_string(std::get<0u>(sbe2)));
        }
        return std::lexicographical_compare(std::get<1u>(sbe1), std::get<2u>(sbe1), std::get<1u>(sbe2),
                                            std::get<2u>(sbe2));
    }

private:
#if !defined(PIRANHA_DOXYGEN_INVOKED)
    // Make friend with the s11n functions.
    template <typename Archive, typename T1, typename S1>
    friend void
    boost::serialization::save(Archive &, const piranha::boost_s11n_key_wrapper<piranha::monomial<T1, S1>> &, unsigned);
    template <typename Archive, typename T1, typename S1>
    friend void boost::serialization::load(Archive &, piranha::boost_s11n_key_wrapper<piranha::monomial<T1, S1>> &,
                                           unsigned);
#endif
#if defined(PIRANHA_WITH_MSGPACK)
    // Enablers for msgpack serialization.
    template <typename Stream>
    using msgpack_pack_enabler = enable_if_t<
        conjunction<is_msgpack_stream<Stream>, has_msgpack_pack<Stream, typename base::container_type>>::value, int>;
    template <typename U>
    using msgpack_convert_enabler = enable_if_t<has_msgpack_convert<typename U::container_type>::value, int>;

public:

    template <typename Stream, msgpack_pack_enabler<Stream> = 0>
    void msgpack_pack(msgpack::packer<Stream> &packer, msgpack_format f, const symbol_fset &s) const
    {
        if (unlikely(this->size() != s.size())) {
            piranha_throw(std::invalid_argument, "incompatible symbol set in monomial serialization: the reference "
                                                 "symbol set has a size of "
                                                     + std::to_string(s.size())
                                                     + ", while the monomial being serialized has a size of "
                                                     + std::to_string(this->size()));
        }
        piranha::msgpack_pack(packer, this->m_container, f);
    }

    template <typename U = monomial, msgpack_convert_enabler<U> = 0>
    void msgpack_convert(const msgpack::object &o, msgpack_format f, const symbol_fset &s)
    {
        piranha::msgpack_convert(this->m_container, o, f);
        if (unlikely(this->size() != s.size())) {
            piranha_throw(std::invalid_argument, "incompatible symbol set in monomial serialization: the reference "
                                                 "symbol set has a size of "
                                                     + std::to_string(s.size())
                                                     + ", while the monomial being deserialized has a size of "
                                                     + std::to_string(this->size()));
        }
    }
#endif
};

template <typename T, typename S>
const std::size_t monomial<T, S>::multiply_arity;
}

// Implementation of the Boost s11n api.
namespace boost
{
namespace serialization
{

template <typename Archive, typename T, typename S>
inline void save(Archive &ar, const piranha::boost_s11n_key_wrapper<piranha::monomial<T, S>> &k, unsigned)
{
    if (unlikely(k.key().size() != k.ss().size())) {
        piranha_throw(std::invalid_argument, "incompatible symbol set in monomial serialization: the reference "
                                             "symbol set has a size of "
                                                 + std::to_string(k.ss().size())
                                                 + ", while the monomial being serialized has a size of "
                                                 + std::to_string(k.key().size()));
    }
    piranha::boost_save(ar, k.key().m_container);
}

template <typename Archive, typename T, typename S>
inline void load(Archive &ar, piranha::boost_s11n_key_wrapper<piranha::monomial<T, S>> &k, unsigned)
{
    piranha::boost_load(ar, k.key().m_container);
    if (unlikely(k.key().size() != k.ss().size())) {
        piranha_throw(std::invalid_argument, "incompatible symbol set in monomial serialization: the reference "
                                             "symbol set has a size of "
                                                 + std::to_string(k.ss().size())
                                                 + ", while the monomial being deserialized has a size of "
                                                 + std::to_string(k.key().size()));
    }
}

template <typename Archive, typename T, typename S>
inline void serialize(Archive &ar, piranha::boost_s11n_key_wrapper<piranha::monomial<T, S>> &k, unsigned version)
{
    split_free(ar, k, version);
}
}
}

namespace piranha
{

inline namespace impl
{

// Enablers for the boost s11n methods.
template <typename Archive, typename T, typename S>
using monomial_boost_save_enabler
    = enable_if_t<has_boost_save<Archive, typename monomial<T, S>::container_type>::value>;

template <typename Archive, typename T, typename S>
using monomial_boost_load_enabler
    = enable_if_t<has_boost_load<Archive, typename monomial<T, S>::container_type>::value>;
}


template <typename Archive, typename T, typename S>
struct boost_save_impl<Archive, boost_s11n_key_wrapper<monomial<T, S>>, monomial_boost_save_enabler<Archive, T, S>>
    : boost_save_via_boost_api<Archive, boost_s11n_key_wrapper<monomial<T, S>>> {
};


template <typename Archive, typename T, typename S>
struct boost_load_impl<Archive, boost_s11n_key_wrapper<monomial<T, S>>, monomial_boost_load_enabler<Archive, T, S>>
    : boost_load_via_boost_api<Archive, boost_s11n_key_wrapper<monomial<T, S>>> {
};
}

namespace std
{


template <typename T, typename S>
struct hash<piranha::monomial<T, S>> : public hash<piranha::array_key<T, piranha::monomial<T, S>, S>> {
};
}

#endif