math.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_MATH_HPP
#define PIRANHA_MATH_HPP

#include <algorithm>
#include <boost/numeric/conversion/cast.hpp>
#include <cmath>
#include <complex>
#include <cstdarg>
#include <initializer_list>
#include <iterator>
#include <stdexcept>
#include <string>
#include <type_traits>
#include <unordered_set>
#include <utility>
#include <vector>

#include <piranha/detail/sfinae_types.hpp>
#include <piranha/exceptions.hpp>
#include <piranha/is_key.hpp>
#include <piranha/symbol_utils.hpp>
#include <piranha/type_traits.hpp>

namespace piranha
{


namespace math
{


template <typename T, typename = void>
struct is_zero_impl {
private:
    // NOTE: the equality comparable requirement already implies that the return type of
    // the comparison must be convertible to bool.
    template <typename U>
    using enabler =
        typename std::enable_if<std::is_constructible<U, int>::value && is_equality_comparable<U>::value, int>::type;

public:

    template <typename U, enabler<U> = 0>
    bool operator()(const U &x) const
    {
        return x == U(0);
    }
};
}

namespace detail
{

// Enabler for math::is_zero().
template <typename T>
using math_is_zero_enabler = typename std::enable_if<
    std::is_convertible<decltype(math::is_zero_impl<T>{}(std::declval<const T &>())), bool>::value, int>::type;
}

namespace math
{


template <typename T, detail::math_is_zero_enabler<T> = 0>
inline bool is_zero(const T &x)
{
    return is_zero_impl<T>{}(x);
}
}

namespace detail
{

// Enabler for the std complex specialisation of is_zero.
template <typename T>
using math_is_zero_std_complex_enabler =
    typename std::enable_if<std::is_same<T, std::complex<float>>::value || std::is_same<T, std::complex<double>>::value
                            || std::is_same<T, std::complex<long double>>::value>::type;
}

namespace math
{


template <typename T>
struct is_zero_impl<T, detail::math_is_zero_std_complex_enabler<T>> {

    bool operator()(const T &c) const
    {
        return is_zero(c.real()) && is_zero(c.imag());
    }
};


template <typename T, typename = void>
struct is_unitary_impl {
private:
    template <typename U>
    using enabler =
        typename std::enable_if<std::is_constructible<U, int>::value && is_equality_comparable<U>::value, int>::type;

public:

    template <typename U, enabler<U> = 0>
    bool operator()(const U &x) const
    {
        return x == U(1);
    }
};
}

namespace detail
{

// Enabler for piranha::math::is_unitary().
template <typename T>
using math_is_unitary_enabler = typename std::enable_if<
    std::is_convertible<decltype(math::is_unitary_impl<T>{}(std::declval<const T &>())), bool>::value, int>::type;
}

namespace math
{


template <typename T, detail::math_is_unitary_enabler<T> = 0>
inline bool is_unitary(const T &x)
{
    return is_unitary_impl<T>{}(x);
}


template <typename T, typename = void>
struct negate_impl {
private:
    // NOTE: inside this type trait, U is always a non-const non-reference type.
    template <typename U>
    using generic_enabler =
        typename std::enable_if<!std::is_integral<U>::value
                                    && detail::true_tt<decltype(std::declval<U &>() = -std::declval<U &>())>::value,
                                int>::type;
    template <typename U>
    using integral_enabler = typename std::enable_if<std::is_integral<U>::value, int>::type;

public:

    template <typename U, generic_enabler<U> = 0>
    void operator()(U &x) const
    {
        x = -x;
    }

    template <typename U, integral_enabler<U> = 0>
    void operator()(U &x) const
    {
        // NOTE: here we use the explicit static_cast to cope with integral promotions
        // (e.g., in case of char).
        x = static_cast<U>(-x);
    }
};
}

namespace detail
{

// Enabler for math::negate().
template <typename T>
using math_negate_enabler = typename std::enable_if<
    !std::is_const<T>::value && true_tt<decltype(math::negate_impl<T>{}(std::declval<T &>()))>::value, int>::type;
}

namespace math
{


template <typename T, detail::math_negate_enabler<T> = 0>
inline void negate(T &x)
{
    negate_impl<T>{}(x);
}


template <typename T, typename = void>
struct multiply_accumulate_impl {
private:
    // NOTE: as usual, we check the expression against const ref arguments.
    template <typename U>
    using addmul_t = decltype(std::declval<U &>() += std::declval<const U &>() * std::declval<const U &>());
    template <typename U>
    using enabler = enable_if_t<is_detected<addmul_t, U>::value, int>;

public:

    template <typename U, enabler<U> = 0>
    void operator()(U &x, const U &y, const U &z) const
    {
        x += y * z;
    }
};

#if defined(FP_FAST_FMA) && defined(FP_FAST_FMAF) && defined(FP_FAST_FMAL)

inline namespace impl
{

// Enabler for the fast floating-point implementation of multiply_accumulate().
template <typename T>
using math_multiply_accumulate_float_enabler = enable_if_t<std::is_floating_point<T>::value>;
}


template <typename T>
struct multiply_accumulate_impl<T, math_multiply_accumulate_float_enabler<T>> {

    void operator()(T &x, const T &y, const T &z) const
    {
        x = std::fma(y, z, x);
    }
};

#endif
}

inline namespace impl
{

// Enabler for multiply_accumulate.
template <typename T>
using math_multiply_accumulate_t = decltype(
    math::multiply_accumulate_impl<T>{}(std::declval<T &>(), std::declval<const T &>(), std::declval<const T &>()));

template <typename T>
using math_multiply_accumulate_enabler = enable_if_t<is_detected<math_multiply_accumulate_t, T>::value, int>;
}

namespace math
{


template <typename T, math_multiply_accumulate_enabler<T> = 0>
inline void multiply_accumulate(T &x, const T &y, const T &z)
{
    multiply_accumulate_impl<T>{}(x, y, z);
}


template <typename T, typename Enable = void>
struct cos_impl {
};


template <typename T>
struct cos_impl<T, typename std::enable_if<std::is_floating_point<T>::value>::type> {

    T operator()(const T &x) const
    {
        return std::cos(x);
    }
};


template <typename T>
struct cos_impl<T, typename std::enable_if<std::is_integral<T>::value>::type> {

    T operator()(const T &x) const
    {
        if (x == T(0)) {
            return T(1);
        }
        piranha_throw(std::invalid_argument, "cannot compute the cosine of a non-zero integral");
    }
};
}

namespace detail
{

// Type for the result of math::cos().
template <typename T>
using math_cos_type_ = decltype(math::cos_impl<T>{}(std::declval<const T &>()));

template <typename T>
using math_cos_type = typename std::enable_if<is_returnable<math_cos_type_<T>>::value, math_cos_type_<T>>::type;
}

namespace math
{


template <typename T>
inline detail::math_cos_type<T> cos(const T &x)
{
    return cos_impl<T>{}(x);
}


template <typename T, typename Enable = void>
struct sin_impl {
};


template <typename T>
struct sin_impl<T, typename std::enable_if<std::is_floating_point<T>::value>::type> {

    T operator()(const T &x) const
    {
        return std::sin(x);
    }
};


template <typename T>
struct sin_impl<T, typename std::enable_if<std::is_integral<T>::value>::type> {

    T operator()(const T &x) const
    {
        if (x == T(0)) {
            return T(0);
        }
        piranha_throw(std::invalid_argument, "cannot compute the sine of a non-zero integral");
    }
};
}

namespace detail
{

// Type for the result of math::sin().
template <typename T>
using math_sin_type_ = decltype(math::sin_impl<T>{}(std::declval<const T &>()));

template <typename T>
using math_sin_type = typename std::enable_if<is_returnable<math_sin_type_<T>>::value, math_sin_type_<T>>::type;
}

namespace math
{


template <typename T>
inline detail::math_sin_type<T> sin(const T &x)
{
    return sin_impl<T>{}(x);
}


template <typename T, typename Enable = void>
struct partial_impl {
};


template <typename T>
struct partial_impl<T, typename std::enable_if<std::is_arithmetic<T>::value>::type> {

    T operator()(const T &, const std::string &) const
    {
        return T(0);
    }
};
}

namespace detail
{

// Return type for math::partial().
template <typename T>
using math_partial_type_
    = decltype(math::partial_impl<T>{}(std::declval<const T &>(), std::declval<const std::string &>()));

template <typename T>
using math_partial_type =
    typename std::enable_if<is_returnable<math_partial_type_<T>>::value, math_partial_type_<T>>::type;
}

namespace math
{


template <typename T>
inline detail::math_partial_type<T> partial(const T &x, const std::string &str)
{
    return partial_impl<T>{}(x, str);
}


template <typename T, typename Enable = void>
struct integrate_impl {
};
}

namespace detail
{

// Return type for math::integrate().
template <typename T>
using math_integrate_type_
    = decltype(math::integrate_impl<T>{}(std::declval<const T &>(), std::declval<const std::string &>()));

template <typename T>
using math_integrate_type =
    typename std::enable_if<is_returnable<math_integrate_type_<T>>::value, math_integrate_type_<T>>::type;
}

namespace math
{


template <typename T>
inline detail::math_integrate_type<T> integrate(const T &x, const std::string &str)
{
    return integrate_impl<T>{}(x, str);
}


template <typename T, typename U, typename = void>
class evaluate_impl
{
    template <typename T1, typename U1>
    using ret_t = typename std::conditional<is_detected<add_t, T1, U1>::value, detected_t<add_t, T1, U1>, T1>::type;
    template <typename T1, typename U1>
    using ret_type = enable_if_t<
        conjunction<is_returnable<ret_t<T1, U1>>, std::is_constructible<ret_t<T1, U1>, const T1 &>>::value,
        ret_t<T1, U1>>;

public:

    template <typename T1, typename U1>
    ret_type<T1, U1> operator()(const T1 &x, const symbol_fmap<U1> &) const
    {
        return ret_type<T1, U1>(x);
    }
};
}

inline namespace impl
{

// Return type for math::evaluate().
template <typename T, typename U>
using math_evaluate_t_
    = decltype(math::evaluate_impl<T, U>{}(std::declval<const T &>(), std::declval<const symbol_fmap<U> &>()));

template <typename T, typename U>
using math_evaluate_t = enable_if_t<is_returnable<math_evaluate_t_<T, U>>::value, math_evaluate_t_<T, U>>;
}

namespace math
{


template <typename U, typename T>
inline math_evaluate_t<T, U> evaluate(const T &x, const symbol_fmap<U> &dict)
{
    return evaluate_impl<T, U>{}(x, dict);
}


template <typename T, typename U, typename Enable = void>
struct subs_impl {
};
}

namespace detail
{

// Return type for math::subs().
template <typename T, typename U>
using math_subs_type_
    = decltype(math::subs_impl<T, U>{}(std::declval<const T &>(), std::declval<const symbol_fmap<U> &>()));

template <typename T, typename U>
using math_subs_type = enable_if_t<is_returnable<math_subs_type_<T, U>>::value, math_subs_type_<T, U>>;
}

namespace math
{


template <typename U, typename T>
inline detail::math_subs_type<T, U> subs(const T &x, const symbol_fmap<U> &dict)
{
    return subs_impl<T, U>{}(x, dict);
}


template <typename T, typename U, typename V, typename = void>
struct t_subs_impl {
};
}

namespace detail
{

// Return type for math::t_subs().
template <typename T, typename U, typename V>
using math_t_subs_type_
    = decltype(math::t_subs_impl<T, U, V>{}(std::declval<const T &>(), std::declval<const std::string &>(),
                                            std::declval<const U &>(), std::declval<const V &>()));

template <typename T, typename U, typename V>
using math_t_subs_type =
    typename std::enable_if<is_returnable<math_t_subs_type_<T, U, V>>::value, math_t_subs_type_<T, U, V>>::type;
}

namespace math
{


template <typename T, typename U, typename V>
inline detail::math_t_subs_type<T, U, V> t_subs(const T &x, const std::string &name, const U &c, const V &s)
{
    return t_subs_impl<T, U, V>{}(x, name, c, s);
}


template <typename T, typename Enable = void>
struct abs_impl {
};
}

namespace detail
{

template <typename T>
using abs_arith_enabler = typename std::enable_if<(std::is_signed<T>::value && std::is_integral<T>::value)
                                                  || (std::is_unsigned<T>::value && std::is_integral<T>::value)
                                                  || std::is_floating_point<T>::value>::type;
}

namespace math
{


template <typename T>
struct abs_impl<T, detail::abs_arith_enabler<T>> {
private:
    template <typename U>
    static U impl(const U &x, typename std::enable_if<std::is_floating_point<U>::value>::type * = nullptr)
    {
        return std::abs(x);
    }
    // NOTE: use decltype here so, in the remote case we are dealing with an extended integer types (for which
    // std::abs() will not exist and cast to long long might be lossy), the overload will be discarded.
    template <typename U>
    static auto impl(const U &x,
                     typename std::enable_if<std::is_integral<U>::value && std::is_signed<U>::value>::type * = nullptr)
        -> decltype(static_cast<U>(std::abs(static_cast<long long>(x))))
    {
        // Cast to long long in order to avoid conversion derps when dealing with chars.
        return static_cast<U>(std::abs(static_cast<long long>(x)));
    }
    template <typename U>
    static U impl(const U &x,
                  typename std::enable_if<std::is_integral<U>::value && std::is_unsigned<U>::value>::type * = nullptr)
    {
        return x;
    }

public:

    auto operator()(const T &x) const -> decltype(impl(x))
    {
        return impl(x);
    }
};


template <typename T>
inline auto abs(const T &x) -> decltype(abs_impl<T>()(x))
{
    return abs_impl<T>()(x);
}
}


template <typename T>
class has_abs
{
    template <typename U>
    using abs_t = decltype(math::abs(std::declval<const U &>()));
    static const bool implementation_defined = is_detected<abs_t, T>::value;

public:
    static const bool value = implementation_defined;
};

// Static init.
template <typename T>
const bool has_abs<T>::value;


template <typename T>
class has_is_zero
{
    template <typename U>
    using is_zero_t = decltype(math::is_zero(std::declval<const U &>()));
    static const bool implementation_defined = is_detected<is_zero_t, T>::value;

public:
    static const bool value = implementation_defined;
};

// Static init.
template <typename T>
const bool has_is_zero<T>::value;


template <typename T>
class has_negate
{
    template <typename U>
    using negate_t = decltype(math::negate(std::declval<U &>()));
    static const bool implementation_defined = is_detected<negate_t, T>::value;

public:
    static const bool value = implementation_defined;
};

template <typename T>
const bool has_negate<T>::value;

namespace detail
{

#if !defined(PIRANHA_DOXYGEN_INVOKED)

// Type definition and type checking for the output of Poisson brackets.
template <typename T>
using pbracket_type_tmp = decltype(math::partial(std::declval<const T &>(), std::string())
                                   * math::partial(std::declval<const T &>(), std::string()));

template <typename T, typename = void>
struct pbracket_type_ {
};

template <typename T>
struct pbracket_type_<
    T, typename std::enable_if<std::is_same<decltype(std::declval<const pbracket_type_tmp<T> &>()
                                                     + std::declval<const pbracket_type_tmp<T> &>()),
                                            pbracket_type_tmp<T>>::value
                               && std::is_same<decltype(std::declval<const pbracket_type_tmp<T> &>()
                                                        - std::declval<const pbracket_type_tmp<T> &>()),
                                               pbracket_type_tmp<T>>::value
                               && std::is_constructible<pbracket_type_tmp<T>, int>::value
                               && std::is_assignable<pbracket_type_tmp<T> &, pbracket_type_tmp<T>>::value>::type> {
    using type = pbracket_type_tmp<T>;
};

// The final typedef.
template <typename T>
using pbracket_type = typename pbracket_type_<T>::type;

#endif
}

namespace math
{


template <typename T>
inline detail::pbracket_type<T> pbracket(const T &f, const T &g, const std::vector<std::string> &p_list,
                                         const std::vector<std::string> &q_list)
{
    using return_type = detail::pbracket_type<T>;
    if (p_list.size() != q_list.size()) {
        piranha_throw(std::invalid_argument, "the number of coordinates is different from the number of momenta");
    }
    if (std::unordered_set<std::string>(p_list.begin(), p_list.end()).size() != p_list.size()) {
        piranha_throw(std::invalid_argument, "the list of momenta contains duplicate entries");
    }
    if (std::unordered_set<std::string>(q_list.begin(), q_list.end()).size() != q_list.size()) {
        piranha_throw(std::invalid_argument, "the list of coordinates contains duplicate entries");
    }
    return_type retval = return_type(0);
    for (decltype(p_list.size()) i = 0u; i < p_list.size(); ++i) {
        // NOTE: could use multadd/sub here, if we implement it for series.
        retval = retval + partial(f, q_list[i]) * partial(g, p_list[i]);
        retval = retval - partial(f, p_list[i]) * partial(g, q_list[i]);
    }
    return retval;
}
}


template <typename T>
class has_pbracket : detail::sfinae_types
{
    using v_string = std::vector<std::string>;
    template <typename T1>
    static auto test(const T1 &x)
        -> decltype(math::pbracket(x, x, std::declval<v_string const &>(), std::declval<v_string const &>()), void(),
                    yes());
    static no test(...);

public:
    static const bool value = std::is_same<decltype(test(std::declval<T>())), yes>::value;
};

template <typename T>
const bool has_pbracket<T>::value;

namespace detail
{

template <typename T>
using is_canonical_enabler = typename std::enable_if<has_pbracket<T>::value && has_is_zero<pbracket_type<T>>::value
                                                         && std::is_constructible<pbracket_type<T>, int>::value
                                                         && is_equality_comparable<pbracket_type<T>>::value,
                                                     int>::type;

template <typename T>
inline bool is_canonical_impl(const std::vector<T const *> &new_p, const std::vector<T const *> &new_q,
                              const std::vector<std::string> &p_list, const std::vector<std::string> &q_list)
{
    using p_type = decltype(math::pbracket(*new_q[0], *new_p[0], p_list, q_list));
    if (p_list.size() != q_list.size()) {
        piranha_throw(std::invalid_argument, "the number of coordinates is different from the number of momenta");
    }
    if (new_p.size() != new_q.size()) {
        piranha_throw(std::invalid_argument,
                      "the number of new coordinates is different from the number of new momenta");
    }
    if (p_list.size() != new_p.size()) {
        piranha_throw(std::invalid_argument, "the number of new momenta is different from the number of momenta");
    }
    if (std::unordered_set<std::string>(p_list.begin(), p_list.end()).size() != p_list.size()) {
        piranha_throw(std::invalid_argument, "the list of momenta contains duplicate entries");
    }
    if (std::unordered_set<std::string>(q_list.begin(), q_list.end()).size() != q_list.size()) {
        piranha_throw(std::invalid_argument, "the list of coordinates contains duplicate entries");
    }
    const auto size = new_p.size();
    for (decltype(new_p.size()) i = 0u; i < size; ++i) {
        for (decltype(new_p.size()) j = 0u; j < size; ++j) {
            // NOTE: no need for actually doing computations when i == j.
            if (i != j && !math::is_zero(math::pbracket(*new_p[i], *new_p[j], p_list, q_list))) {
                return false;
            }
            if (i != j && !math::is_zero(math::pbracket(*new_q[i], *new_q[j], p_list, q_list))) {
                return false;
            }
            // Poisson bracket needs to be zero for i != j, one for i == j.
            // NOTE: cast from bool to int is always 0 or 1.
            if (math::pbracket(*new_q[i], *new_p[j], p_list, q_list) != p_type(static_cast<int>(i == j))) {
                return false;
            }
        }
    }
    return true;
}
}

namespace math
{


template <typename T, detail::is_canonical_enabler<T> = 0>
inline bool transformation_is_canonical(const std::vector<T> &new_p, const std::vector<T> &new_q,
                                        const std::vector<std::string> &p_list, const std::vector<std::string> &q_list)
{
    std::vector<T const *> pv, qv;
    std::transform(new_p.begin(), new_p.end(), std::back_inserter(pv), [](const T &p) { return &p; });
    std::transform(new_q.begin(), new_q.end(), std::back_inserter(qv), [](const T &q) { return &q; });
    return detail::is_canonical_impl(pv, qv, p_list, q_list);
}

// clang-format off

// clang-format on
template <typename T, detail::is_canonical_enabler<T> = 0>
inline bool transformation_is_canonical(std::initializer_list<T> new_p, std::initializer_list<T> new_q,
                                        const std::vector<std::string> &p_list, const std::vector<std::string> &q_list)
{
    std::vector<T const *> pv, qv;
    std::transform(new_p.begin(), new_p.end(), std::back_inserter(pv), [](const T &p) { return &p; });
    std::transform(new_q.begin(), new_q.end(), std::back_inserter(qv), [](const T &q) { return &q; });
    return detail::is_canonical_impl(pv, qv, p_list, q_list);
}


template <typename T, typename Enable = void>
struct degree_impl {
};


template <typename T>
inline auto degree(const T &x) -> decltype(degree_impl<T>()(x))
{
    return degree_impl<T>()(x);
}


template <typename T>
inline auto degree(const T &x, const symbol_fset &s) -> decltype(degree_impl<T>()(x, s))
{
    return degree_impl<T>()(x, s);
}


template <typename T, typename Enable = void>
struct ldegree_impl {
};


template <typename T>
inline auto ldegree(const T &x) -> decltype(ldegree_impl<T>()(x))
{
    return ldegree_impl<T>()(x);
}


template <typename T>
inline auto ldegree(const T &x, const symbol_fset &s) -> decltype(ldegree_impl<T>()(x, s))
{
    return ldegree_impl<T>()(x, s);
}


template <typename T, typename Enable = void>
struct t_degree_impl {
};


template <typename T>
inline auto t_degree(const T &x) -> decltype(t_degree_impl<T>{}(x))
{
    return t_degree_impl<T>{}(x);
}


template <typename T>
inline auto t_degree(const T &x, const symbol_fset &names) -> decltype(t_degree_impl<T>{}(x, names))
{
    return t_degree_impl<T>{}(x, names);
}


template <typename T, typename Enable = void>
struct t_ldegree_impl {
};


template <typename T>
inline auto t_ldegree(const T &x) -> decltype(t_ldegree_impl<T>{}(x))
{
    return t_ldegree_impl<T>{}(x);
}


template <typename T>
inline auto t_ldegree(const T &x, const symbol_fset &names) -> decltype(t_ldegree_impl<T>{}(x, names))
{
    return t_ldegree_impl<T>{}(x, names);
}


template <typename T, typename Enable = void>
struct t_order_impl {
};


template <typename T>
inline auto t_order(const T &x) -> decltype(t_order_impl<T>{}(x))
{
    return t_order_impl<T>{}(x);
}


template <typename T>
inline auto t_order(const T &x, const symbol_fset &names) -> decltype(t_order_impl<T>{}(x, names))
{
    return t_order_impl<T>{}(x, names);
}


template <typename T, typename Enable = void>
struct t_lorder_impl {
};


template <typename T>
inline auto t_lorder(const T &x) -> decltype(t_lorder_impl<T>{}(x))
{
    return t_lorder_impl<T>{}(x);
}


template <typename T>
inline auto t_lorder(const T &x, const symbol_fset &names) -> decltype(t_lorder_impl<T>{}(x, names))
{
    return t_lorder_impl<T>{}(x, names);
}


template <typename T, typename U, typename = void>
struct truncate_degree_impl {
};
}

namespace detail
{

// Enablers for the degree truncation methods.
template <typename T, typename U>
using truncate_degree_enabler = typename std::enable_if<
    std::is_same<decltype(math::truncate_degree_impl<T, U>()(std::declval<const T &>(), std::declval<const U &>())),
                 T>::value,
    int>::type;

template <typename T, typename U>
using truncate_pdegree_enabler = typename std::enable_if<
    std::is_same<decltype(math::truncate_degree_impl<T, U>()(std::declval<const T &>(), std::declval<const U &>(),
                                                             std::declval<const symbol_fset &>())),
                 T>::value,
    int>::type;
}

namespace math
{


template <typename T, typename U, detail::truncate_degree_enabler<T, U> = 0>
inline T truncate_degree(const T &x, const U &max_degree)
{
    return truncate_degree_impl<T, U>()(x, max_degree);
}


template <typename T, typename U, detail::truncate_pdegree_enabler<T, U> = 0>
inline T truncate_degree(const T &x, const U &max_degree, const symbol_fset &names)
{
    return truncate_degree_impl<T, U>()(x, max_degree, names);
}
}


template <typename T, typename U>
class has_truncate_degree : detail::sfinae_types
{
    template <typename T1, typename U1>
    static auto test1(const T1 &t, const U1 &u) -> decltype(math::truncate_degree(t, u), void(), yes());
    static no test1(...);
    template <typename T1, typename U1>
    static auto test2(const T1 &t, const U1 &u)
        -> decltype(math::truncate_degree(t, u, std::declval<const symbol_fset &>()), void(), yes());
    static no test2(...);

public:
    static const bool value = std::is_same<decltype(test1(std::declval<T>(), std::declval<U>())), yes>::value
                              && std::is_same<decltype(test2(std::declval<T>(), std::declval<U>())), yes>::value;
};

template <typename T, typename U>
const bool has_truncate_degree<T, U>::value;


template <typename T>
class is_differentiable : detail::sfinae_types
{
    template <typename U>
    static auto test(const U &u) -> decltype(math::partial(u, ""), void(), yes());
    static no test(...);
    static const bool implementation_defined = std::is_same<decltype(test(std::declval<T>())), yes>::value;

public:
    static const bool value = implementation_defined;
};

// Static init.
template <typename T>
const bool is_differentiable<T>::value;

inline namespace impl
{

// A key is differentiable/integrable if the type returned by the partial/integrate method
// is a pair (sometype,key).
template <typename, typename>
struct is_differential_key_pair : std::false_type {
};

template <typename Key, typename T>
struct is_differential_key_pair<Key, std::pair<T, Key>> : std::true_type {
};
}


template <typename Key>
class key_is_differentiable
{
    PIRANHA_TT_CHECK(is_key, uncvref_t<Key>);
    template <typename U>
    using key_partial_t = decltype(
        std::declval<const U &>().partial(std::declval<const symbol_idx &>(), std::declval<const symbol_fset &>()));
    static const bool implementation_defined
        = is_differential_key_pair<uncvref_t<Key>, detected_t<key_partial_t, Key>>::value;

public:
    static const bool value = implementation_defined;
};

template <typename Key>
const bool key_is_differentiable<Key>::value;


template <typename T>
class is_integrable : detail::sfinae_types
{
    template <typename U>
    static auto test(const U &u) -> decltype(math::integrate(u, ""), void(), yes());
    static no test(...);
    static const bool implementation_defined = std::is_same<decltype(test(std::declval<T>())), yes>::value;

public:
    static const bool value = implementation_defined;
};

template <typename T>
const bool is_integrable<T>::value;


template <typename Key>
class key_is_integrable
{
    PIRANHA_TT_CHECK(is_key, uncvref_t<Key>);
    template <typename U>
    using key_integrate_t = decltype(
        std::declval<const U &>().integrate(std::declval<const std::string &>(), std::declval<const symbol_fset &>()));
    static const bool implementation_defined
        = is_differential_key_pair<uncvref_t<Key>, detected_t<key_integrate_t, Key>>::value;

public:
    static const bool value = implementation_defined;
};

template <typename T>
const bool key_is_integrable<T>::value;


template <typename T>
class has_degree : detail::sfinae_types
{
    template <typename U>
    static auto test1(const U &u) -> decltype(math::degree(u), void(), yes());
    static no test1(...);
    template <typename U>
    static auto test2(const U &u) -> decltype(math::degree(u, std::declval<const symbol_fset &>()), void(), yes());
    static no test2(...);

public:
    static const bool value = std::is_same<decltype(test1(std::declval<T>())), yes>::value
                              && std::is_same<decltype(test2(std::declval<T>())), yes>::value;
};

// Static init.
template <typename T>
const bool has_degree<T>::value;


template <typename T>
class has_ldegree : detail::sfinae_types
{
    template <typename U>
    static auto test1(const U &u) -> decltype(math::ldegree(u), void(), yes());
    static no test1(...);
    template <typename U>
    static auto test2(const U &u) -> decltype(math::ldegree(u, std::declval<const symbol_fset &>()), void(), yes());
    static no test2(...);

public:
    static const bool value = std::is_same<decltype(test1(std::declval<T>())), yes>::value
                              && std::is_same<decltype(test2(std::declval<T>())), yes>::value;
};

// Static init.
template <typename T>
const bool has_ldegree<T>::value;


template <typename T>
class has_t_degree : detail::sfinae_types
{
    template <typename U>
    static auto test1(const U &u) -> decltype(math::t_degree(u), void(), yes());
    static no test1(...);
    template <typename U>
    static auto test2(const U &u) -> decltype(math::t_degree(u, std::declval<const symbol_fset &>()), void(), yes());
    static no test2(...);

public:
    static const bool value = std::is_same<decltype(test1(std::declval<T>())), yes>::value
                              && std::is_same<decltype(test2(std::declval<T>())), yes>::value;
};

// Static init.
template <typename T>
const bool has_t_degree<T>::value;


template <typename T>
class has_t_ldegree : detail::sfinae_types
{
    template <typename U>
    static auto test1(const U &u) -> decltype(math::t_ldegree(u), void(), yes());
    static no test1(...);
    template <typename U>
    static auto test2(const U &u) -> decltype(math::t_ldegree(u, std::declval<const symbol_fset &>()), void(), yes());
    static no test2(...);

public:
    static const bool value = std::is_same<decltype(test1(std::declval<T>())), yes>::value
                              && std::is_same<decltype(test2(std::declval<T>())), yes>::value;
};

// Static init.
template <typename T>
const bool has_t_ldegree<T>::value;


template <typename T>
class has_t_order : detail::sfinae_types
{
    template <typename U>
    static auto test1(const U &u) -> decltype(math::t_order(u), void(), yes());
    static no test1(...);
    template <typename U>
    static auto test2(const U &u) -> decltype(math::t_order(u, std::declval<const symbol_fset &>()), void(), yes());
    static no test2(...);

public:
    static const bool value = std::is_same<decltype(test1(std::declval<T>())), yes>::value
                              && std::is_same<decltype(test2(std::declval<T>())), yes>::value;
};

// Static init.
template <typename T>
const bool has_t_order<T>::value;


template <typename T>
class has_t_lorder : detail::sfinae_types
{
    template <typename U>
    static auto test1(const U &u) -> decltype(math::t_lorder(u), void(), yes());
    static no test1(...);
    template <typename U>
    static auto test2(const U &u) -> decltype(math::t_lorder(u, std::declval<const symbol_fset &>()), void(), yes());
    static no test2(...);

public:
    static const bool value = std::is_same<decltype(test1(std::declval<T>())), yes>::value
                              && std::is_same<decltype(test2(std::declval<T>())), yes>::value;
};

// Static init.
template <typename T>
const bool has_t_lorder<T>::value;


template <typename Key>
class key_has_degree
{
    PIRANHA_TT_CHECK(is_key, uncvref_t<Key>);
    template <typename U>
    using total_degree_t = decltype(std::declval<const U &>().degree(std::declval<const symbol_fset &>()));
    template <typename U>
    using partial_degree_t = decltype(
        std::declval<const U &>().degree(std::declval<const symbol_idx_fset &>(), std::declval<const symbol_fset &>()));
    static const bool implementation_defined
        = conjunction<is_detected<total_degree_t, Key>, is_detected<partial_degree_t, Key>>::value;

public:
    static const bool value = implementation_defined;
};

template <typename Key>
const bool key_has_degree<Key>::value;


template <typename Key>
class key_has_ldegree
{
    PIRANHA_TT_CHECK(is_key, uncvref_t<Key>);
    template <typename U>
    using total_ldegree_t = decltype(std::declval<const U &>().ldegree(std::declval<const symbol_fset &>()));
    template <typename U>
    using partial_ldegree_t = decltype(std::declval<const U &>().ldegree(std::declval<const symbol_idx_fset &>(),
                                                                         std::declval<const symbol_fset &>()));
    static const bool implementation_defined
        = conjunction<is_detected<total_ldegree_t, Key>, is_detected<partial_ldegree_t, Key>>::value;

public:
    static const bool value = implementation_defined;
};

template <typename Key>
const bool key_has_ldegree<Key>::value;


template <typename Key>
class key_has_t_degree
{
    PIRANHA_TT_CHECK(is_key, uncvref_t<Key>);
    template <typename U>
    using total_t_degree_t = decltype(std::declval<const U &>().t_degree(std::declval<const symbol_fset &>()));
    template <typename U>
    using partial_t_degree_t = decltype(std::declval<const U &>().t_degree(std::declval<const symbol_idx_fset &>(),
                                                                           std::declval<const symbol_fset &>()));
    static const bool implementation_defined
        = conjunction<is_detected<total_t_degree_t, Key>, is_detected<partial_t_degree_t, Key>>::value;

public:
    static const bool value = implementation_defined;
};

// Static init.
template <typename T>
const bool key_has_t_degree<T>::value;


template <typename Key>
class key_has_t_ldegree
{
    PIRANHA_TT_CHECK(is_key, uncvref_t<Key>);
    template <typename U>
    using total_t_ldegree_t = decltype(std::declval<const U &>().t_ldegree(std::declval<const symbol_fset &>()));
    template <typename U>
    using partial_t_ldegree_t = decltype(std::declval<const U &>().t_ldegree(std::declval<const symbol_idx_fset &>(),
                                                                             std::declval<const symbol_fset &>()));
    static const bool implementation_defined
        = conjunction<is_detected<total_t_ldegree_t, Key>, is_detected<partial_t_ldegree_t, Key>>::value;

public:
    static const bool value = implementation_defined;
};

// Static init.
template <typename T>
const bool key_has_t_ldegree<T>::value;


template <typename Key>
class key_has_t_order
{
    PIRANHA_TT_CHECK(is_key, uncvref_t<Key>);
    template <typename U>
    using total_t_order_t = decltype(std::declval<const U &>().t_order(std::declval<const symbol_fset &>()));
    template <typename U>
    using partial_t_order_t = decltype(std::declval<const U &>().t_order(std::declval<const symbol_idx_fset &>(),
                                                                         std::declval<const symbol_fset &>()));
    static const bool implementation_defined
        = conjunction<is_detected<total_t_order_t, Key>, is_detected<partial_t_order_t, Key>>::value;

public:
    static const bool value = implementation_defined;
};

// Static init.
template <typename T>
const bool key_has_t_order<T>::value;


template <typename Key>
class key_has_t_lorder
{
    PIRANHA_TT_CHECK(is_key, uncvref_t<Key>);
    template <typename U>
    using total_t_lorder_t = decltype(std::declval<const U &>().t_lorder(std::declval<const symbol_fset &>()));
    template <typename U>
    using partial_t_lorder_t = decltype(std::declval<const U &>().t_lorder(std::declval<const symbol_idx_fset &>(),
                                                                           std::declval<const symbol_fset &>()));
    static const bool implementation_defined
        = conjunction<is_detected<total_t_lorder_t, Key>, is_detected<partial_t_lorder_t, Key>>::value;

public:
    static const bool value = implementation_defined;
};

// Static init.
template <typename T>
const bool key_has_t_lorder<T>::value;


template <typename T>
class has_is_unitary : detail::sfinae_types
{
    typedef typename std::decay<T>::type Td;
    template <typename T1>
    static auto test(const T1 &t) -> decltype(math::is_unitary(t), void(), yes());
    static no test(...);
    static const bool implementation_defined = std::is_same<decltype(test(std::declval<Td>())), yes>::value;

public:
    static const bool value = implementation_defined;
};

template <typename T>
const bool has_is_unitary<T>::value;


template <typename T, typename U>
class has_subs : detail::sfinae_types
{
    typedef typename std::decay<T>::type Td;
    typedef typename std::decay<U>::type Ud;
    template <typename T1, typename U1>
    static auto test(const T1 &t, const U1 &u)
        -> decltype(math::subs(t, std::declval<const symbol_fmap<U1> &>()), void(), yes());
    static no test(...);
    static const bool implementation_defined
        = std::is_same<decltype(test(std::declval<Td>(), std::declval<Ud>())), yes>::value;

public:
    static const bool value = implementation_defined;
};

// Static init.
template <typename T, typename U>
const bool has_subs<T, U>::value;


template <typename T, typename U, typename V>
class has_t_subs : detail::sfinae_types
{
    typedef typename std::decay<T>::type Td;
    typedef typename std::decay<U>::type Ud;
    typedef typename std::decay<V>::type Vd;
    template <typename T1, typename U1, typename V1>
    static auto test(const T1 &t, const U1 &u, const V1 &v)
        -> decltype(math::t_subs(t, std::declval<std::string const &>(), u, v), void(), yes());
    static no test(...);
    static const bool implementation_defined
        = std::is_same<decltype(test(std::declval<Td>(), std::declval<Ud>(), std::declval<Vd>())), yes>::value;

public:
    static const bool value = implementation_defined;
};

// Static init.
template <typename T, typename U, typename V>
const bool has_t_subs<T, U, V>::value;

inline namespace impl
{

// Check the result of key subs methods.
template <typename, typename>
struct key_subs_check_type : std::false_type {
};

template <typename Key, typename T>
struct key_subs_check_type<Key, std::vector<std::pair<T, Key>>> : std::true_type {
};
}


template <typename Key, typename T, typename U>
class key_has_t_subs
{
    PIRANHA_TT_CHECK(is_key, uncvref_t<Key>);
    template <typename Key1, typename T1, typename U1>
    using t_subs_t = decltype(
        std::declval<const Key1 &>().t_subs(std::declval<const symbol_idx &>(), std::declval<const T1 &>(),
                                            std::declval<const U1 &>(), std::declval<const symbol_fset &>()));
    static const bool implementation_defined
        = key_subs_check_type<uncvref_t<Key>, detected_t<t_subs_t, Key, T, U>>::value;

public:
    static const bool value = implementation_defined;
};

// Static init.
template <typename Key, typename T, typename U>
const bool key_has_t_subs<Key, T, U>::value;


template <typename Key, typename T>
class key_has_subs
{
    PIRANHA_TT_CHECK(is_key, uncvref_t<Key>);
    template <typename Key1, typename T1>
    using subs_t = decltype(std::declval<const Key1 &>().subs(std::declval<const symbol_idx_fmap<T1> &>(),
                                                              std::declval<const symbol_fset &>()));
    static const bool implementation_defined = key_subs_check_type<uncvref_t<Key>, detected_t<subs_t, Key, T>>::value;

public:
    static const bool value = implementation_defined;
};

// Static init.
template <typename Key, typename T>
const bool key_has_subs<Key, T>::value;


template <typename T>
class has_multiply_accumulate
{
    template <typename U>
    using multiply_accumulate_t = decltype(
        math::multiply_accumulate(std::declval<U &>(), std::declval<const U &>(), std::declval<const U &>()));
    static const bool implementation_defined = is_detected<multiply_accumulate_t, T>::value;

public:
    static const bool value = implementation_defined;
};

// Static init.
template <typename T>
const bool has_multiply_accumulate<T>::value;


template <typename T, typename U>
class is_evaluable
{
    template <typename T1, typename U1>
    using eval_t = decltype(math::evaluate(std::declval<const T1 &>(), std::declval<const symbol_fmap<U1> &>()));
    static const bool implementation_defined = is_detected<eval_t, T, U>::value;

public:
    static const bool value = implementation_defined;
};

template <typename T, typename U>
const bool is_evaluable<T, U>::value;


template <typename Key, typename T>
class key_is_evaluable
{
    PIRANHA_TT_CHECK(is_key, uncvref_t<Key>);
    template <typename Key1, typename T1>
    using evaluate_t = decltype(std::declval<const Key1 &>().evaluate(std::declval<const std::vector<T1> &>(),
                                                                      std::declval<const symbol_fset &>()));
    static const bool implementation_defined = is_detected<evaluate_t, Key, T>::value;

public:
    static const bool value = implementation_defined;
};

template <typename Key, typename T>
const bool key_is_evaluable<Key, T>::value;


template <typename T>
class has_sine : detail::sfinae_types
{
    template <typename T1>
    static auto test(const T1 &x) -> decltype(math::sin(x), void(), yes());
    static no test(...);

public:
    static const bool value = std::is_same<decltype(test(std::declval<T>())), yes>::value;
};

template <typename T>
const bool has_sine<T>::value;


template <typename T>
class has_cosine : detail::sfinae_types
{
    template <typename T1>
    static auto test(const T1 &x) -> decltype(math::cos(x), void(), yes());
    static no test(...);

public:
    static const bool value = std::is_same<decltype(test(std::declval<T>())), yes>::value;
};

template <typename T>
const bool has_cosine<T>::value;


template <typename T>
class has_transformation_is_canonical : detail::sfinae_types
{
    using v_string = std::vector<std::string>;
    template <typename T1>
    static auto test(const T1 &) -> decltype(math::transformation_is_canonical(std::declval<std::vector<T1> const &>(),
                                                                               std::declval<std::vector<T1> const &>(),
                                                                               std::declval<v_string const &>(),
                                                                               std::declval<v_string const &>()),
                                             void(), yes());
    static no test(...);

public:
    static const bool value = std::is_same<decltype(test(std::declval<T>())), yes>::value;
};

template <typename T>
const bool has_transformation_is_canonical<T>::value;

namespace math
{


template <typename T, typename Enable = void>
struct add3_impl {

    template <typename U>
    auto operator()(U &a, const U &b, const U &c) const -> decltype(a = b + c)
    {
        return a = b + c;
    }
};


template <typename T>
struct add3_impl<T, typename std::enable_if<std::is_integral<T>::value>::type> {

    T &operator()(T &a, const T &b, const T &c) const
    {
        return a = static_cast<T>(b + c);
    }
};


template <typename T>
inline auto add3(T &a, const T &b, const T &c) -> decltype(add3_impl<T>()(a, b, c))
{
    return add3_impl<T>()(a, b, c);
}


template <typename T, typename Enable = void>
struct sub3_impl {

    template <typename U>
    auto operator()(U &a, const U &b, const U &c) const -> decltype(a = b - c)
    {
        return a = b - c;
    }
};


template <typename T>
struct sub3_impl<T, typename std::enable_if<std::is_integral<T>::value>::type> {

    T &operator()(T &a, const T &b, const T &c) const
    {
        return a = static_cast<T>(b - c);
    }
};


template <typename T>
inline auto sub3(T &a, const T &b, const T &c) -> decltype(sub3_impl<T>()(a, b, c))
{
    return sub3_impl<T>()(a, b, c);
}


template <typename T, typename Enable = void>
struct mul3_impl {

    template <typename U>
    auto operator()(U &a, const U &b, const U &c) const -> decltype(a = b * c)
    {
        return a = b * c;
    }
};


template <typename T>
struct mul3_impl<T, typename std::enable_if<std::is_integral<T>::value>::type> {

    T &operator()(T &a, const T &b, const T &c) const
    {
        return a = static_cast<T>(b * c);
    }
};


template <typename T>
inline auto mul3(T &a, const T &b, const T &c) -> decltype(mul3_impl<T>()(a, b, c))
{
    return mul3_impl<T>()(a, b, c);
}


template <typename T, typename Enable = void>
struct div3_impl {

    template <typename U>
    auto operator()(U &a, const U &b, const U &c) const -> decltype(a = b / c)
    {
        return a = b / c;
    }
};


template <typename T>
struct div3_impl<T, typename std::enable_if<std::is_integral<T>::value>::type> {

    T &operator()(T &a, const T &b, const T &c) const
    {
        return a = static_cast<T>(b / c);
    }
};


template <typename T>
inline auto div3(T &a, const T &b, const T &c) -> decltype(div3_impl<T>()(a, b, c))
{
    return div3_impl<T>()(a, b, c);
}
}

namespace detail
{

// Greatest common divisor using the euclidean algorithm.
// NOTE: this can yield negative values, depending on the signs
// of a and b. Supports C++ integrals and mp_integer.
// NOTE: using this with C++ integrals unchecked on ranges can result in undefined
// behaviour.
template <typename T>
inline T gcd_euclidean(T a, T b)
{
    while (true) {
        if (math::is_zero(a)) {
            return b;
        }
        b %= a;
        if (math::is_zero(b)) {
            return a;
        }
        a %= b;
    }
}
}

namespace math
{


template <typename T, typename U, typename = void>
struct gcd_impl {
};


template <typename T, typename U>
struct gcd_impl<T, U, typename std::enable_if<std::is_integral<T>::value && std::is_integral<U>::value>::type> {
    using p_type = decltype(std::declval<const T &>() + std::declval<const U &>());

    p_type operator()(const T &a, const U &b) const
    {
        return detail::gcd_euclidean(static_cast<p_type>(a), static_cast<p_type>(b));
    }
};


template <typename T, typename U>
inline auto gcd(const T &a, const U &b) -> decltype(gcd_impl<T, U>()(a, b))
{
    return gcd_impl<T, U>()(a, b);
}


template <typename T, typename = void>
struct gcd3_impl {

    template <typename T1>
    auto operator()(T1 &out, const T1 &a, const T1 &b) const -> decltype(out = math::gcd(a, b))
    {
        return out = math::gcd(a, b);
    }
};


template <typename T>
struct gcd3_impl<T, typename std::enable_if<std::is_integral<T>::value>::type> {

    T &operator()(T &out, const T &a, const T &b) const
    {
        return out = static_cast<T>(math::gcd(a, b));
    }
};


template <typename T>
inline auto gcd3(T &out, const T &a, const T &b) -> decltype(gcd3_impl<T>()(out, a, b))
{
    return gcd3_impl<T>()(out, a, b);
}
}


template <typename T>
class has_add3
{
    template <typename U>
    using add3_t = decltype(math::add3(std::declval<U &>(), std::declval<const U &>(), std::declval<const U &>()));
    static const bool implementation_defined = is_detected<add3_t, T>::value;

public:
    static const bool value = implementation_defined;
};

template <typename T>
const bool has_add3<T>::value;


template <typename T>
class has_sub3
{
    template <typename U>
    using sub3_t = decltype(math::sub3(std::declval<U &>(), std::declval<const U &>(), std::declval<const U &>()));
    static const bool implementation_defined = is_detected<sub3_t, T>::value;

public:
    static const bool value = implementation_defined;
};

template <typename T>
const bool has_sub3<T>::value;


template <typename T>
class has_mul3
{
    template <typename U>
    using mul3_t = decltype(math::mul3(std::declval<U &>(), std::declval<const U &>(), std::declval<const U &>()));
    static const bool implementation_defined = is_detected<mul3_t, T>::value;

public:
    static const bool value = implementation_defined;
};

template <typename T>
const bool has_mul3<T>::value;


template <typename T>
class has_div3
{
    template <typename U>
    using div3_t = decltype(math::div3(std::declval<U &>(), std::declval<const U &>(), std::declval<const U &>()));
    static const bool implementation_defined = is_detected<div3_t, T>::value;

public:
    static const bool value = implementation_defined;
};

template <typename T>
const bool has_div3<T>::value;


template <typename T, typename U = T>
class has_gcd
{
    template <typename T1, typename U1>
    using gcd_t = decltype(math::gcd(std::declval<const T1 &>(), std::declval<const U1 &>()));
    static const bool implementation_defined = is_detected<gcd_t, T, U>::value;

public:
    static const bool value = implementation_defined;
};

template <typename T, typename U>
const bool has_gcd<T, U>::value;


template <typename T>
class has_gcd3
{
    template <typename U>
    using gcd3_t = decltype(math::gcd3(std::declval<U &>(), std::declval<const U &>(), std::declval<const U &>()));
    static const bool implementation_defined = is_detected<gcd3_t, T>::value;

public:
    static const bool value = implementation_defined;
};

template <typename T>
const bool has_gcd3<T>::value;
}

#endif