real_trigonometric_kronecker_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_REAL_TRIGONOMETRIC_KRONECKER_MONOMIAL_HPP
#define PIRANHA_REAL_TRIGONOMETRIC_KRONECKER_MONOMIAL_HPP

#include <algorithm>
#include <array>
#include <cinttypes>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <cstdlib>
#include <functional>
#include <initializer_list>
#include <iostream>
#include <iterator>
#include <limits>
#include <stdexcept>
#include <string>
#include <type_traits>
#include <unordered_map>
#include <utility>
#include <vector>

#include <piranha/binomial.hpp>
#include <piranha/config.hpp>
#include <piranha/detail/cf_mult_impl.hpp>
#include <piranha/detail/km_commons.hpp>
#include <piranha/detail/prepare_for_print.hpp>
#include <piranha/detail/safe_integral_adder.hpp>
#include <piranha/exceptions.hpp>
#include <piranha/is_cf.hpp>
#include <piranha/is_key.hpp>
#include <piranha/kronecker_array.hpp>
#include <piranha/math.hpp>
#include <piranha/mp_integer.hpp>
#include <piranha/s11n.hpp>
#include <piranha/safe_cast.hpp>
#include <piranha/static_vector.hpp>
#include <piranha/symbol_utils.hpp>
#include <piranha/term.hpp>
#include <piranha/type_traits.hpp>

namespace piranha
{


// NOTES:
// - it might make sense, for canonicalisation and is_compatible(), to provide a method in kronecker_array to get only
//   the first element of the array. This should be quite fast, and it will provide enough information for the
//   canon/compatibility.
// - related to the above: we can embed the flavour as the first element of the kronecker array - at that point checking
//   the flavour is just determining if the int value is even or odd.
// - need to do some reasoning on the impact of the codification in a specialised fast poisson series multiplier: how do
//   we deal with the canonical form without going through code/decode? if we require the last multiplier to be always
//   positive (instead of the first), can we just check/flip the sign of the coded value?
template <typename T = std::make_signed<std::size_t>::type>
class real_trigonometric_kronecker_monomial
{
public:
    typedef T value_type;

private:
    typedef kronecker_array<value_type> ka;

public:
    static const std::size_t multiply_arity = 2u;

    typedef typename ka::size_type size_type;
    static const size_type max_size = 255u;
    typedef static_vector<value_type, max_size> v_type;

private:
    static_assert(max_size <= std::numeric_limits<static_vector<int, 1u>::size_type>::max(), "Invalid max size.");

public:

    real_trigonometric_kronecker_monomial() : m_value(0), m_flavour(true) {}
    real_trigonometric_kronecker_monomial(const real_trigonometric_kronecker_monomial &) = default;
    real_trigonometric_kronecker_monomial(real_trigonometric_kronecker_monomial &&) = default;

private:
    // Enabler for ctor from init list.
    template <typename U>
    using init_list_enabler = enable_if_t<has_safe_cast<value_type, U>::value, int>;

public:

    template <typename U, init_list_enabler<U> = 0>
    explicit real_trigonometric_kronecker_monomial(std::initializer_list<U> list, bool flavour = true)
        : real_trigonometric_kronecker_monomial(std::begin(list), std::end(list), flavour)
    {
    }

private:
    // Enabler for ctor from iterator.
    template <typename Iterator>
    using it_ctor_enabler
        = enable_if_t<conjunction<is_input_iterator<Iterator>,
                                  has_safe_cast<value_type, decltype(*std::declval<const Iterator &>())>>::value,
                      int>;

public:

    template <typename Iterator, it_ctor_enabler<Iterator> = 0>
    explicit real_trigonometric_kronecker_monomial(const Iterator &start, const Iterator &end, bool flavour = true)
        : m_value(0), m_flavour(flavour)
    {
        v_type tmp;
        std::transform(
            start, end, std::back_inserter(tmp),
            [](const typename std::iterator_traits<Iterator>::value_type &v) { return safe_cast<value_type>(v); });
        m_value = ka::encode(tmp);
    }

    explicit real_trigonometric_kronecker_monomial(const symbol_fset &) : real_trigonometric_kronecker_monomial() {}

    explicit real_trigonometric_kronecker_monomial(const value_type &n, bool f) : m_value(n), m_flavour(f) {}

    explicit real_trigonometric_kronecker_monomial(const real_trigonometric_kronecker_monomial &other,
                                                   const symbol_fset &)
        : real_trigonometric_kronecker_monomial(other)
    {
    }
    ~real_trigonometric_kronecker_monomial()
    {
        PIRANHA_TT_CHECK(is_key, real_trigonometric_kronecker_monomial);
        PIRANHA_TT_CHECK(key_has_t_degree, real_trigonometric_kronecker_monomial);
        PIRANHA_TT_CHECK(key_has_t_ldegree, real_trigonometric_kronecker_monomial);
        PIRANHA_TT_CHECK(key_has_t_order, real_trigonometric_kronecker_monomial);
        PIRANHA_TT_CHECK(key_has_t_lorder, real_trigonometric_kronecker_monomial);
        PIRANHA_TT_CHECK(key_is_differentiable, real_trigonometric_kronecker_monomial);
    }

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

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

    void set_int(const value_type &n)
    {
        m_value = n;
    }

    value_type get_int() const
    {
        return m_value;
    }

    bool get_flavour() const
    {
        return m_flavour;
    }

    void set_flavour(bool f)
    {
        m_flavour = f;
    }

private:
    // Implementation of canonicalisation.
    static bool canonicalise_impl(v_type &unpacked)
    {
        const auto size = unpacked.size();
        bool sign_change = false;
        for (decltype(unpacked.size()) i = 0u; i < size; ++i) {
            if (sign_change || unpacked[i] < value_type(0)) {
                unpacked[i] = static_cast<value_type>(-unpacked[i]);
                sign_change = true;
            } else if (unpacked[i] > value_type(0)) {
                break;
            }
        }
        return sign_change;
    }

public:

    bool canonicalise(const symbol_fset &args)
    {
        auto unpacked = unpack(args);
        const bool retval = canonicalise_impl(unpacked);
        if (retval) {
            m_value = ka::encode(unpacked);
        }
        return retval;
    }

    bool is_compatible(const symbol_fset &args) const noexcept
    {
        const auto s = args.size();
        // No args means the value must also be zero.
        if (s == 0u) {
            return !m_value;
        }
        const auto &limits = ka::get_limits();
        // If we overflow the maximum size available, we cannot use this object as key in series.
        if (s >= limits.size()) {
            return false;
        }
        const auto &l = limits[s];
        // Value is compatible if it is within the bounds for the given size.
        if (m_value < std::get<1u>(l) || m_value > std::get<2u>(l)) {
            return false;
        }
        // Now check for the first multiplier.
        // NOTE: here we have already checked all the conditions that could lead to unpack() throwing, so
        // we do not need to put @throw specifications in the doc.
        const auto unpacked = unpack(args);
        // We know that s != 0.
        piranha_assert(unpacked.size() > 0u);
        for (typename v_type::size_type i = 0u; i < s; ++i) {
            if (unpacked[i] < value_type(0)) {
                return false;
            } else if (unpacked[i] > value_type(0)) {
                break;
            }
        }
        return true;
    }

    bool is_zero(const symbol_fset &) const noexcept
    {
        return m_value == value_type(0) && !m_flavour;
    }

    real_trigonometric_kronecker_monomial merge_symbols(const symbol_idx_fmap<symbol_fset> &ins_map,
                                                        const symbol_fset &args) const
    {
        return real_trigonometric_kronecker_monomial(detail::km_merge_symbols<v_type, ka>(ins_map, args, m_value),
                                                     m_flavour);
    }

    bool is_unitary(const symbol_fset &) const
    {
        return (!m_value && m_flavour);
    }

private:
    // Degree utils.
    using degree_type = add_t<T, T>;

public:

    degree_type t_degree(const symbol_fset &args) const
    {
        const auto tmp = unpack(args);
        // NOTE: this should be guaranteed by the unpack function.
        piranha_assert(tmp.size() == args.size());
        degree_type retval(0);
        for (const auto &x : tmp) {
            detail::safe_integral_adder(retval, static_cast<degree_type>(x));
        }
        return retval;
    }

    degree_type t_ldegree(const symbol_fset &args) const
    {
        return t_degree(args);
    }

    degree_type t_degree(const symbol_idx_fset &p, const symbol_fset &args) const
    {
        const auto tmp = unpack(args);
        piranha_assert(tmp.size() == args.size());
        if (unlikely(p.size() && *p.rbegin() >= tmp.size())) {
            piranha_throw(std::invalid_argument,
                          "the largest value in the positions set for the computation of the "
                          "partial trigonometric degree of a real trigonometric Kronecker monomial is "
                              + std::to_string(*p.rbegin()) + ", but the monomial has a size of only "
                              + std::to_string(tmp.size()));
        }
        degree_type retval(0);
        for (auto idx : p) {
            detail::safe_integral_adder(retval, static_cast<degree_type>(tmp[static_cast<decltype(tmp.size())>(idx)]));
        }
        return retval;
    }

    degree_type t_ldegree(const symbol_idx_fset &p, const symbol_fset &args) const
    {
        return t_degree(p, args);
    }

private:
    // Order utils.
    using order_type = decltype(math::abs(std::declval<const T &>()) + math::abs(std::declval<const T &>()));

public:

    order_type t_order(const symbol_fset &args) const
    {
        const auto tmp = unpack(args);
        piranha_assert(tmp.size() == args.size());
        order_type retval(0);
        for (const auto &x : tmp) {
            // NOTE: here the k codification is symmetric, we can always take the negative
            // safely.
            detail::safe_integral_adder(retval, static_cast<order_type>(math::abs(x)));
        }
        return retval;
    }

    order_type t_lorder(const symbol_fset &args) const
    {
        return t_order(args);
    }

    order_type t_order(const symbol_idx_fset &p, const symbol_fset &args) const
    {
        const auto tmp = unpack(args);
        piranha_assert(tmp.size() == args.size());
        if (unlikely(p.size() && *p.rbegin() >= tmp.size())) {
            piranha_throw(std::invalid_argument,
                          "the largest value in the positions set for the computation of the "
                          "partial trigonometric order of a real trigonometric Kronecker monomial is "
                              + std::to_string(*p.rbegin()) + ", but the monomial has a size of only "
                              + std::to_string(tmp.size()));
        }
        order_type retval(0);
        for (auto idx : p) {
            detail::safe_integral_adder(
                retval, static_cast<degree_type>(math::abs(tmp[static_cast<decltype(tmp.size())>(idx)])));
        }
        return retval;
    }

    order_type t_lorder(const symbol_idx_fset &p, const symbol_fset &args) const
    {
        return t_order(p, args);
    }

private:
    // Enabler for multiplication.
    template <typename Cf>
    using multiply_enabler = enable_if_t<
        conjunction<is_divisible_in_place<Cf, int>, has_negate<Cf>, std::is_copy_assignable<Cf>, has_mul3<Cf>>::value,
        int>;

public:

    template <typename Cf, multiply_enabler<Cf> = 0>
    static void multiply(std::array<term<Cf, real_trigonometric_kronecker_monomial>, multiply_arity> &res,
                         const term<Cf, real_trigonometric_kronecker_monomial> &t1,
                         const term<Cf, real_trigonometric_kronecker_monomial> &t2, const symbol_fset &args)
    {
        // Coefficients first.
        cf_mult_impl(res[0u].m_cf, t1.m_cf, t2.m_cf);
        res[1u].m_cf = res[0u].m_cf;
        const bool f1 = t1.m_key.get_flavour(), f2 = t2.m_key.get_flavour();
        if (f1 && f2) {
            // cos, cos: no change.
        } else if (!f1 && !f2) {
            // sin, sin: negate the plus.
            math::negate(res[0u].m_cf);
        } else if (!f1 && f2) {
            // sin, cos: no change.
        } else {
            // cos, sin: negate the minus.
            math::negate(res[1u].m_cf);
        }
        // Now the keys.
        auto &retval_plus = res[0u].m_key;
        auto &retval_minus = res[1u].m_key;
        // Flags to signal if a sign change in the multipliers was needed as part of the canonicalization.
        bool sign_plus = false, sign_minus = false;
        const auto size = args.size();
        const auto tmp1 = t1.m_key.unpack(args), tmp2 = t2.m_key.unpack(args);
        v_type result_plus, result_minus;
        for (typename v_type::size_type i = 0u; i < size; ++i) {
            result_plus.push_back(tmp1[i]);
            detail::safe_integral_adder(result_plus[i], tmp2[i]);
            // NOTE: it is safe here to take the negative because in kronecker_array we are guaranteed
            // that the range of each element is symmetric, so if tmp2[i] is representable also -tmp2[i] is.
            // NOTE: the static cast here is because if value_type is narrower than int, the unary minus will promote
            // to int and safe_adder won't work as it expects identical types.
            result_minus.push_back(tmp1[i]);
            detail::safe_integral_adder(result_minus[i], static_cast<value_type>(-tmp2[i]));
        }
        // Handle sign changes.
        sign_plus = canonicalise_impl(result_plus);
        sign_minus = canonicalise_impl(result_minus);
        // Compute them before assigning, so in case of exceptions we do not touch the return values.
        const auto re_plus = ka::encode(result_plus), re_minus = ka::encode(result_minus);
        retval_plus.m_value = re_plus;
        retval_minus.m_value = re_minus;
        const bool f = (t1.m_key.get_flavour() == t2.m_key.get_flavour());
        retval_plus.m_flavour = f;
        retval_minus.m_flavour = f;
        // If multiplier sign was changed and the result is a sine, negate the coefficient.
        if (sign_plus && !res[0u].m_key.get_flavour()) {
            math::negate(res[0u].m_cf);
        }
        if (sign_minus && !res[1u].m_key.get_flavour()) {
            math::negate(res[1u].m_cf);
        }
    }

    std::size_t hash() const
    {
        return static_cast<std::size_t>(m_value);
    }

    bool operator==(const real_trigonometric_kronecker_monomial &other) const
    {
        return (m_value == other.m_value && m_flavour == other.m_flavour);
    }

    bool operator!=(const real_trigonometric_kronecker_monomial &other) const
    {
        return (m_value != other.m_value || m_flavour != other.m_flavour);
    }

    v_type unpack(const symbol_fset &args) const
    {
        return detail::km_unpack<v_type, ka>(args, m_value);
    }

    void print(std::ostream &os, const symbol_fset &args) const
    {
        // Don't print anything in case all multipliers are zero.
        if (m_value == value_type(0)) {
            return;
        }
        if (m_flavour) {
            os << "cos(";
        } else {
            os << "sin(";
        }
        const auto tmp = unpack(args);
        piranha_assert(tmp.size() == args.size());
        const value_type zero(0), one(1), m_one(-1);
        bool empty_output = true;
        auto it_args = args.begin();
        for (decltype(tmp.size()) i = 0u; i < tmp.size(); ++i, ++it_args) {
            if (tmp[i] != zero) {
                // A positive multiplier, in case previous output exists, must be preceded
                // by a "+" sign.
                if (tmp[i] > zero && !empty_output) {
                    os << '+';
                }
                // Print the multiplier, unless it's "-1": in that case, just print the minus sign.
                if (tmp[i] == m_one) {
                    os << '-';
                } else if (tmp[i] != one) {
                    os << detail::prepare_for_print(tmp[i]) << '*';
                }
                // Finally, print name of variable.
                os << *it_args;
                empty_output = false;
            }
        }
        os << ")";
    }

    void print_tex(std::ostream &os, const symbol_fset &args) const
    {
        // Don't print anything in case all multipliers are zero.
        if (m_value == value_type(0)) {
            return;
        }
        if (m_flavour) {
            os << "\\cos{\\left(";
        } else {
            os << "\\sin{\\left(";
        }
        const auto tmp = unpack(args);
        piranha_assert(tmp.size() == args.size());
        const value_type zero(0), one(1), m_one(-1);
        bool empty_output = true;
        auto it_args = args.begin();
        for (decltype(tmp.size()) i = 0u; i < tmp.size(); ++i, ++it_args) {
            if (tmp[i] != zero) {
                // A positive multiplier, in case previous output exists, must be preceded
                // by a "+" sign.
                if (tmp[i] > zero && !empty_output) {
                    os << '+';
                }
                // Print the multiplier, unless it's "-1": in that case, just print the minus sign.
                if (tmp[i] == m_one) {
                    os << '-';
                } else if (tmp[i] != one) {
                    os << static_cast<long long>(tmp[i]);
                }
                // Finally, print name of variable.
                os << '{' << *it_args << '}';
                empty_output = false;
            }
        }
        os << "\\right)}";
    }

    std::pair<T, real_trigonometric_kronecker_monomial> partial(const symbol_idx &p, const symbol_fset &args) const
    {
        auto v = unpack(args);
        if (p >= args.size() || v[static_cast<decltype(v.size())>(p)] == T(0)) {
            // Derivative wrt a variable not in the monomial: the position is outside the bounds, or it refers to a
            // variable with zero multiplier.
            return std::make_pair(T(0), real_trigonometric_kronecker_monomial{args});
        }
        const auto v_b = v.begin();
        if (get_flavour()) {
            // cos(nx + b) -> -n*sin(nx + b)
            return std::make_pair(static_cast<T>(-v_b[p]), real_trigonometric_kronecker_monomial(m_value, false));
        }
        // sin(nx + b) -> n*cos(nx + b)
        return std::make_pair(v_b[p], real_trigonometric_kronecker_monomial(m_value, true));
    }

    std::pair<T, real_trigonometric_kronecker_monomial> integrate(const std::string &s, const symbol_fset &args) const
    {
        auto v = unpack(args);
        auto it_args = args.begin();
        for (min_int<decltype(args.size()), typename v_type::size_type> i = 0u; i < args.size(); ++i, ++it_args) {
            if (*it_args == s && v[i] != value_type(0)) {
                if (get_flavour()) {
                    // cos(nx + b) -> sin(nx + b)
                    return std::make_pair(v[i], real_trigonometric_kronecker_monomial(m_value, false));
                }
                // sin(nx + b) -> -cos(nx + b)
                return std::make_pair(static_cast<value_type>(-v[i]),
                                      real_trigonometric_kronecker_monomial(m_value, true));
            }
            if (*it_args > s) {
                // The current symbol in args is lexicographically larger than s,
                // we won't find s any more in args.
                break;
            }
        }
        return std::make_pair(value_type(0), real_trigonometric_kronecker_monomial{args});
    }

private:
    // The candidate type, resulting from math::cos()/math::sin() on T * U.
    template <typename U>
    using eval_t_cos = decltype(math::cos(std::declval<const mul_t<T, U> &>()));
    template <typename U>
    using eval_t_sin = decltype(math::sin(std::declval<const mul_t<T, U> &>()));
    // Definition of the evaluation type.
    template <typename U>
    using eval_type
        = enable_if_t<conjunction<
                          // sin/cos eval types must be the same.
                          std::is_same<eval_t_cos<U>, eval_t_sin<U>>,
                          // Eval type must be ctible form int.
                          std::is_constructible<eval_t_cos<U>, const int &>,
                          // Eval type must be returnable.
                          is_returnable<eval_t_cos<U>>,
                          // T * U must be addable in-place.
                          is_addable_in_place<mul_t<T, U>>,
                          // T * U must be ctible from int.
                          std::is_constructible<mul_t<T, U>, const int &>,
                          // T * U must be move ctible and dtible
                          std::is_move_constructible<mul_t<T, U>>, std::is_destructible<mul_t<T, U>>>::value,
                      eval_t_cos<U>>;

public:

    template <typename U>
    eval_type<U> evaluate(const std::vector<U> &values, const symbol_fset &args) const
    {
        // NOTE: here we can check the values size only against args.
        if (unlikely(values.size() != args.size())) {
            piranha_throw(std::invalid_argument, "invalid vector of values for real trigonometric Kronecker monomial "
                                                 "evaluation: the size of the vector of values ("
                                                     + std::to_string(values.size())
                                                     + ") differs from the size of the reference set of symbols ("
                                                     + std::to_string(args.size()) + ")");
        }
        // Run the unpack before the checks below in order to check the suitability of args.
        const auto v = unpack(args);
        // Special casing if the monomial is empty, just return 0 or 1.
        if (!args.size()) {
            if (get_flavour()) {
                return eval_type<U>(1);
            }
            return eval_type<U>(0);
        }
        // Init the accumulator with the first element of the linear
        // combination.
        auto tmp(v[0] * values[0]);
        // Accumulate the sin/cos argument.
        for (decltype(values.size()) i = 1; i < values.size(); ++i) {
            // NOTE: this could be optimised with an FMA eventually.
            tmp += v[static_cast<decltype(v.size())>(i)] * values[i];
        }
        if (get_flavour()) {
            // NOTE: move it, in the future math::cos() may exploit this.
            return math::cos(std::move(tmp));
        }
        return math::sin(std::move(tmp));
    }

private:
    // Substitution utils.
    // NOTE: the only additional requirement here is to be able to negate.
    template <typename U>
    using subs_type = enable_if_t<has_negate<eval_type<U>>::value, eval_type<U>>;

public:

    template <typename U>
    std::vector<std::pair<subs_type<U>, real_trigonometric_kronecker_monomial>> subs(const symbol_idx_fmap<U> &smap,
                                                                                     const symbol_fset &args) const
    {
        if (unlikely(smap.size() && smap.rbegin()->first >= args.size())) {
            // The last element of the substitution map must be a valid index into args.
            piranha_throw(std::invalid_argument, "invalid argument(s) for substitution in a real trigonometric "
                                                 "Kronecker 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(args.size()) + ")");
        }
        std::vector<std::pair<subs_type<U>, real_trigonometric_kronecker_monomial>> retval;
        retval.reserve(2u);
        const auto f = get_flavour();
        if (smap.size()) {
            // The substitution map contains something, proceed to the substitution.
            auto v = unpack(args);
            // Init a tmp value from the linear combination of the first value
            // of the map with the first multiplier.
            auto it = smap.begin();
            auto tmp(v[static_cast<decltype(v.size())>(it->first)] * it->second);
            // Zero out the corresponding multiplier.
            v[static_cast<decltype(v.size())>(it->first)] = T(0);
            // Finish computing the linear combination of the values with
            // the corresponding multipliers.
            // NOTE: move to the next element in the init statement of the for loop.
            for (++it; it != smap.end(); ++it) {
                // NOTE: FMA in the future, maybe.
                tmp += v[static_cast<decltype(v.size())>(it->first)] * it->second;
                v[static_cast<decltype(v.size())>(it->first)] = T(0);
            }
            // Check if we need to canonicalise the vector, after zeroing out
            // one or more multipliers.
            const bool sign_changed = canonicalise_impl(v);
            // Encode the modified vector of multipliers.
            const auto new_value = ka::encode(v);
            // Pre-compute the sin/cos of tmp.
            auto s_tmp(math::sin(tmp)), c_tmp(math::cos(tmp));
            if (f) {
                // cos(tmp+x) -> cos(tmp)*cos(x) - sin(tmp)*sin(x)
                retval.emplace_back(std::move(c_tmp), real_trigonometric_kronecker_monomial(new_value, true));
                if (!sign_changed) {
                    // NOTE: we need to negate because of trig formulas, but if
                    // we switched the sign of x we have a double negation.
                    math::negate(s_tmp);
                }
                retval.emplace_back(std::move(s_tmp), real_trigonometric_kronecker_monomial(new_value, false));
            } else {
                // sin(tmp+x) -> sin(tmp)*cos(x) + cos(tmp)*sin(x)
                retval.emplace_back(std::move(s_tmp), real_trigonometric_kronecker_monomial(new_value, true));
                if (sign_changed) {
                    // NOTE: opposite of the above, if we switched the signs in x we need to negate
                    // cos(tmp) to restore the signs.
                    math::negate(c_tmp);
                }
                retval.emplace_back(std::move(c_tmp), real_trigonometric_kronecker_monomial(new_value, false));
            }
        } else {
            // The subs map is empty, return the original values.
            // NOTE: can we just return a 1-element vector here?
            if (f) {
                // cos(a) -> 1*cos(a) + 0*sin(a)
                retval.emplace_back(subs_type<U>(1), *this);
                retval.emplace_back(subs_type<U>(0), real_trigonometric_kronecker_monomial(m_value, false));
            } else {
                // sin(a) -> 0*cos(a) + 1*sin(a)
                retval.emplace_back(subs_type<U>(0), real_trigonometric_kronecker_monomial(m_value, true));
                retval.emplace_back(subs_type<U>(1), *this);
            }
        }
        return retval;
    }

private:
    // Trig subs bits.
    template <typename U>
    using t_subs_t = decltype(std::declval<const integer &>() * std::declval<const mul_t<U, U> &>());
    template <typename U>
    using t_subs_type = enable_if_t<
        conjunction<std::is_default_constructible<U>, std::is_constructible<U, const int &>, std::is_move_assignable<U>,
                    std::is_destructible<U>, std::is_assignable<U &, mul_t<U, U> &&>, std::is_destructible<mul_t<U, U>>,
                    std::is_move_constructible<t_subs_t<U>>, std::is_destructible<t_subs_t<U>>,
                    is_addable_in_place<t_subs_t<U>>, has_negate<t_subs_t<U>>, std::is_copy_constructible<t_subs_t<U>>,
                    std::is_move_constructible<mul_t<U, U>>>::value,
        t_subs_t<U>>;
    // Couple of helper functions for Vieta's formulae.
    static value_type cos_phase(const value_type &n)
    {
        piranha_assert(n >= value_type(0));
        constexpr value_type v[4] = {1, 0, -1, 0};
        return v[n % value_type(4)];
    }
    static value_type sin_phase(const value_type &n)
    {
        piranha_assert(n >= value_type(0));
        constexpr value_type v[4] = {0, 1, 0, -1};
        return v[n % value_type(4)];
    }

public:

    template <typename U>
    std::vector<std::pair<t_subs_type<U>, real_trigonometric_kronecker_monomial>>
    t_subs(const symbol_idx &idx, const U &c, const U &s, const symbol_fset &args) const
    {
        auto v = unpack(args);
        T n(0);
        if (idx < args.size()) {
            // NOTE: this will extract the affected multiplier
            // into n and zero out the multiplier in v.
            std::swap(n, v[static_cast<decltype(v.size())>(idx)]);
        }
        const auto abs_n = static_cast<T>(std::abs(n));
        // Prepare the powers of c and s to be used in the multiple angles formulae.
        // NOTE: probably a vector would be better here.
        std::unordered_map<T, U> c_map, s_map;
        // TREQ: U def ctible, ctible from int, move-assignable and dtible.
        c_map[0] = U(1);
        s_map[0] = U(1);
        for (T k = 0; k < abs_n; ++k) {
            // TREQ: U assignable from mult_t<U,U> &&, mult_t<U,U> dtible.
            c_map[static_cast<T>(k + 1)] = c_map[k] * c;
            s_map[static_cast<T>(k + 1)] = s_map[k] * s;
        }
        // Init with the first element in the summation.
        // NOTE: we promote abs_n to integer in these formulae in order to force
        // the use of the binomial() overload for integer. In the future we might want
        // to disable the binomial() overload for integral C++ types.
        // TREQ: t_subs_type<U> move-ctible, dtible.
        t_subs_type<U> cos_nx(cos_phase(abs_n) * math::binomial(integer(abs_n), T(0)) * (c_map[T(0)] * s_map[abs_n])),
            sin_nx(sin_phase(abs_n) * math::binomial(integer(abs_n), T(0)) * (c_map[T(0)] * s_map[abs_n]));
        // Run the main iteration.
        for (T k = 0; k < abs_n; ++k) {
            const auto p = static_cast<T>(abs_n - (k + 1));
            piranha_assert(p >= T(0));
            // TREQ: mult_t<U,U> move ctible.
            const auto tmp = c_map[static_cast<T>(k + 1)] * s_map[p];
            const auto tmp_bin = math::binomial(integer(abs_n), k + T(1));
            // TREQ: t_subs_type<U> addable in-place.
            cos_nx += cos_phase(p) * tmp_bin * tmp;
            sin_nx += sin_phase(p) * tmp_bin * tmp;
        }
        // Change sign as necessary.
        if (abs_n != n) {
            // TREQ: t_subs_type<U> has_negate.
            math::negate(sin_nx);
        }
        // Buld the new keys and canonicalise as needed.
        const bool sign_changed = canonicalise_impl(v);
        // Encode the modified vector.
        const auto new_value = ka::encode(v);
        // Init the retval.
        std::vector<std::pair<t_subs_type<U>, real_trigonometric_kronecker_monomial>> retval;
        retval.reserve(2);
        if (get_flavour()) {
            retval.emplace_back(std::move(cos_nx), real_trigonometric_kronecker_monomial(new_value, true));
            retval.emplace_back(std::move(sin_nx), real_trigonometric_kronecker_monomial(new_value, false));
            // Need to flip the sign on the sin * sin product if sign was not changed.
            if (!sign_changed) {
                math::negate(retval[1u].first);
            }
        } else {
            retval.emplace_back(std::move(sin_nx), real_trigonometric_kronecker_monomial(new_value, true));
            retval.emplace_back(std::move(cos_nx), real_trigonometric_kronecker_monomial(new_value, false));
            // Need to flip the sign on the cos * sin product if sign was changed.
            if (sign_changed) {
                math::negate(retval[1u].first);
            }
        }
        return retval;
    }

    void trim_identify(std::vector<char> &trim_mask, const symbol_fset &args) const
    {
        detail::km_trim_identify<v_type, ka>(trim_mask, args, m_value);
    }

    real_trigonometric_kronecker_monomial trim(const std::vector<char> &trim_mask, const symbol_fset &args) const
    {
        return real_trigonometric_kronecker_monomial(detail::km_trim<v_type, ka>(trim_mask, args, m_value), m_flavour);
    }

    bool operator<(const real_trigonometric_kronecker_monomial &other) const
    {
        if (m_value == other.m_value) {
            return m_flavour < other.m_flavour;
        }
        return m_value < other.m_value;
    }

#if defined(PIRANHA_WITH_MSGPACK)
private:
    template <typename Stream>
    using msgpack_pack_enabler
        = enable_if_t<conjunction<is_msgpack_stream<Stream>, has_msgpack_pack<Stream, T>,
                                  has_msgpack_pack<Stream, bool>, has_msgpack_pack<Stream, v_type>>::value,
                      int>;
    template <typename U>
    using msgpack_convert_enabler
        = enable_if_t<conjunction<has_msgpack_convert<typename U::value_type>, has_msgpack_convert<typename U::v_type>,
                                  has_msgpack_convert<bool>>::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
    {
        packer.pack_array(2);
        if (f == msgpack_format::binary) {
            piranha::msgpack_pack(packer, m_value, f);
            piranha::msgpack_pack(packer, m_flavour, f);
        } else {
            auto tmp = unpack(s);
            piranha::msgpack_pack(packer, tmp, f);
            piranha::msgpack_pack(packer, m_flavour, f);
        }
    }

    template <typename U = real_trigonometric_kronecker_monomial, msgpack_convert_enabler<U> = 0>
    void msgpack_convert(const msgpack::object &o, msgpack_format f, const symbol_fset &s)
    {
        std::array<msgpack::object, 2> tmp;
        o.convert(tmp);
        if (f == msgpack_format::binary) {
            piranha::msgpack_convert(m_value, tmp[0], f);
            piranha::msgpack_convert(m_flavour, tmp[1], f);
        } else {
            v_type tmp_v;
            piranha::msgpack_convert(tmp_v, tmp[0], f);
            if (unlikely(tmp_v.size() != s.size())) {
                piranha_throw(std::invalid_argument,
                              "incompatible symbol set in trigonometric monomial serialization: the reference "
                              "symbol set has a size of "
                                  + std::to_string(s.size())
                                  + ", while the trigonometric monomial being deserialized has a size of "
                                  + std::to_string(tmp_v.size()));
            }
            *this = real_trigonometric_kronecker_monomial(tmp_v.begin(), tmp_v.end());
            piranha::msgpack_convert(m_flavour, tmp[1], f);
        }
    }
#endif

private:
    value_type m_value;
    bool m_flavour;
};

template <typename T>
const std::size_t real_trigonometric_kronecker_monomial<T>::multiply_arity;

using rtk_monomial = real_trigonometric_kronecker_monomial<>;
}

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

template <typename Archive, typename T>
inline void save(Archive &ar,
                 const piranha::boost_s11n_key_wrapper<piranha::real_trigonometric_kronecker_monomial<T>> &k, unsigned)
{
    if (std::is_same<Archive, boost::archive::binary_oarchive>::value) {
        piranha::boost_save(ar, k.key().get_int());
    } else {
        auto tmp = k.key().unpack(k.ss());
        piranha::boost_save(ar, tmp);
    }
    piranha::boost_save(ar, k.key().get_flavour());
}

template <typename Archive, typename T>
inline void load(Archive &ar, piranha::boost_s11n_key_wrapper<piranha::real_trigonometric_kronecker_monomial<T>> &k,
                 unsigned)
{
    if (std::is_same<Archive, boost::archive::binary_iarchive>::value) {
        T value;
        piranha::boost_load(ar, value);
        k.key().set_int(value);
    } else {
        typename piranha::real_trigonometric_kronecker_monomial<T>::v_type tmp;
        piranha::boost_load(ar, tmp);
        if (unlikely(tmp.size() != k.ss().size())) {
            piranha_throw(std::invalid_argument, "invalid size detected in the deserialization of a real Kronercker "
                                                 "trigonometric monomial: the deserialized size is "
                                                     + std::to_string(tmp.size())
                                                     + " but the reference symbol set has a size of "
                                                     + std::to_string(k.ss().size()));
        }
        // NOTE: here the exception safety is basic, as the last boost_load() could fail in principle.
        // It does not really matter much, as there's no real dependency between the multipliers and the flavour,
        // any combination is valid.
        k.key() = piranha::real_trigonometric_kronecker_monomial<T>(tmp.begin(), tmp.end());
    }
    // The flavour loading is common.
    bool f;
    piranha::boost_load(ar, f);
    k.key().set_flavour(f);
}

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

namespace piranha
{

inline namespace impl
{

template <typename Archive, typename T>
using rtk_monomial_boost_save_enabler = enable_if_t<
    conjunction<has_boost_save<Archive, T>, has_boost_save<Archive, bool>,
                has_boost_save<Archive, typename real_trigonometric_kronecker_monomial<T>::v_type>>::value>;

template <typename Archive, typename T>
using rtk_monomial_boost_load_enabler = enable_if_t<
    conjunction<has_boost_load<Archive, T>, has_boost_load<Archive, bool>,
                has_boost_load<Archive, typename real_trigonometric_kronecker_monomial<T>::v_type>>::value>;
}


template <typename Archive, typename T>
struct boost_save_impl<Archive, boost_s11n_key_wrapper<real_trigonometric_kronecker_monomial<T>>,
                       rtk_monomial_boost_save_enabler<Archive, T>>
    : boost_save_via_boost_api<Archive, boost_s11n_key_wrapper<real_trigonometric_kronecker_monomial<T>>> {
};


template <typename Archive, typename T>
struct boost_load_impl<Archive, boost_s11n_key_wrapper<real_trigonometric_kronecker_monomial<T>>,
                       rtk_monomial_boost_load_enabler<Archive, T>>
    : boost_load_via_boost_api<Archive, boost_s11n_key_wrapper<real_trigonometric_kronecker_monomial<T>>> {
};
}

namespace std
{

template <typename T>
struct hash<piranha::real_trigonometric_kronecker_monomial<T>> {
    typedef size_t result_type;
    typedef piranha::real_trigonometric_kronecker_monomial<T> argument_type;

    result_type operator()(const argument_type &a) const
    {
        return a.hash();
    }
};
}

#endif