functions.cpp

#include <symengine/visitor.h>
#include <symengine/symengine_exception.h>
 
namespace SymEngine
{
 
extern RCP<const Basic> &i2;
extern RCP<const Basic> &i3;
extern RCP<const Basic> &i5;
extern RCP<const Basic> &im2;
extern RCP<const Basic> &im3;
extern RCP<const Basic> &im5;
 
RCP<const Basic> sqrt(RCP<const Basic> &arg)
{
    return pow(arg, div(one, i2));
}
RCP<const Basic> cbrt(RCP<const Basic> &arg)
{
    return pow(arg, div(one, i3));
}
 
extern RCP<const Basic> &sq3;
extern RCP<const Basic> &sq2;
extern RCP<const Basic> &sq5;
 
extern RCP<const Basic> &C0;
extern RCP<const Basic> &C1;
extern RCP<const Basic> &C2;
extern RCP<const Basic> &C3;
extern RCP<const Basic> &C4;
extern RCP<const Basic> &C5;
extern RCP<const Basic> &C6;
 
extern RCP<const Basic> &mC0;
extern RCP<const Basic> &mC1;
extern RCP<const Basic> &mC2;
extern RCP<const Basic> &mC3;
extern RCP<const Basic> &mC4;
extern RCP<const Basic> &mC5;
extern RCP<const Basic> &mC6;
 
// sin_table()[n] represents the value of sin(pi*n/12) for n = 0..23
static const RCP<const Basic> *sin_table()
{
    static const RCP<const Basic> table[]
        = {zero, C0,  C1,  C2,  C3,  C4,  one,       C4,  C3,  C2,  C1,  C0,
           zero, mC0, mC1, mC2, mC3, mC4, minus_one, mC4, mC3, mC2, mC1, mC0};
    return table;
}
 
static const umap_basic_basic &inverse_cst()
{
    static const umap_basic_basic inverse_cst_ = {
        {C3, i3},
        {mC3, im3},
        {C2, mul(i2, i2)},
        {mC2, mul(im2, i2)},
        {C4, integer(12)},
        {mC4, integer(-12)},
        {C5, i5},
        {mC5, im5},
        {C6, integer(10)},
        {mC6, integer(-10)},
        {div(one, i2), integer(6)},
        {div(minus_one, i2), integer(-6)},
    };
    return inverse_cst_;
}
 
static const umap_basic_basic &inverse_tct()
{
    static const umap_basic_basic inverse_tct_ = {
        {div(one, sq3), mul(i2, i3)},
        {div(minus_one, sq3), mul(im2, i3)},
        {sq3, i3},
        {mul(minus_one, sq3), im3},
        {add(one, sq2), div(pow(i2, i3), i3)},
        {mul(minus_one, add(one, sq2)), div(pow(i2, i3), im3)},
        {sub(sq2, one), pow(i2, i3)},
        {sub(one, sq2), pow(im2, i3)},
        {sub(i2, sq3), mul(mul(i2, i2), i3)},
        {sub(sq3, i2), mul(mul(im2, i2), i3)},
        {sqrt(add(i5, mul(i2, sqrt(i5)))), div(i5, i2)},
        {mul(minus_one, sqrt(add(i5, mul(i2, sqrt(i5))))), div(im5, i2)},
        {one, pow(i2, i2)},
        {minus_one, mul(minus_one, pow(i2, i2))},
    };
    return inverse_tct_;
}
 
Conjugate::Conjugate(const RCP<const Basic> &arg) : OneArgFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool Conjugate::is_canonical(const RCP<const Basic> &arg) const
{
    if (is_a_Number(*arg)) {
        if (eq(*arg, *ComplexInf)) {
            return true;
        }
        return false;
    }
    if (is_a<Constant>(*arg)) {
        return false;
    }
    if (is_a<Mul>(*arg)) {
        return false;
    }
    if (is_a<Pow>(*arg)) {
        if (is_a<Integer>(*down_cast<const Pow &>(*arg).get_exp())) {
            return false;
        }
    }
    // OneArgFunction classes
    if (is_a<Sign>(*arg) or is_a<Conjugate>(*arg) or is_a<Erf>(*arg)
        or is_a<Erfc>(*arg) or is_a<Gamma>(*arg) or is_a<LogGamma>(*arg)
        or is_a<Abs>(*arg)) {
        return false;
    }
    if (is_a<Sin>(*arg) or is_a<Cos>(*arg) or is_a<Tan>(*arg) or is_a<Cot>(*arg)
        or is_a<Sec>(*arg) or is_a<Csc>(*arg)) {
        return false;
    }
    if (is_a<Sinh>(*arg) or is_a<Cosh>(*arg) or is_a<Tanh>(*arg)
        or is_a<Coth>(*arg) or is_a<Sech>(*arg) or is_a<Csch>(*arg)) {
        return false;
    }
    // TwoArgFunction classes
    if (is_a<KroneckerDelta>(*arg) or is_a<ATan2>(*arg)
        or is_a<LowerGamma>(*arg) or is_a<UpperGamma>(*arg)
        or is_a<Beta>(*arg)) {
        return false;
    }
    // MultiArgFunction class
    if (is_a<LeviCivita>(*arg)) {
        return false;
    }
    return true;
}
 
RCP<const Basic> Conjugate::create(const RCP<const Basic> &arg) const
{
    return conjugate(arg);
}
 
RCP<const Basic> conjugate(const RCP<const Basic> &arg)
{
    if (is_a_Number(*arg)) {
        return down_cast<const Number &>(*arg).conjugate();
    }
    if (is_a<Constant>(*arg) or is_a<Abs>(*arg) or is_a<KroneckerDelta>(*arg)
        or is_a<LeviCivita>(*arg)) {
        return arg;
    }
    if (is_a<Mul>(*arg)) {
        const map_basic_basic &dict = down_cast<const Mul &>(*arg).get_dict();
        map_basic_basic new_dict;
        RCP<const Number> coef = rcp_static_cast<const Number>(
            conjugate(down_cast<const Mul &>(*arg).get_coef()));
        for (const auto &p : dict) {
            if (is_a<Integer>(*p.second)) {
                Mul::dict_add_term_new(outArg(coef), new_dict, p.second,
                                       conjugate(p.first));
            } else {
                Mul::dict_add_term_new(
                    outArg(coef), new_dict, one,
                    conjugate(Mul::from_dict(one, {{p.first, p.second}})));
            }
        }
        return Mul::from_dict(coef, std::move(new_dict));
    }
    if (is_a<Pow>(*arg)) {
        RCP<const Basic> base = down_cast<const Pow &>(*arg).get_base();
        RCP<const Basic> exp = down_cast<const Pow &>(*arg).get_exp();
        if (is_a<Integer>(*exp)) {
            return pow(conjugate(base), exp);
        }
    }
    if (is_a<Conjugate>(*arg)) {
        return down_cast<const Conjugate &>(*arg).get_arg();
    }
    if (is_a<Sign>(*arg) or is_a<Erf>(*arg) or is_a<Erfc>(*arg)
        or is_a<Gamma>(*arg) or is_a<LogGamma>(*arg) or is_a<Sin>(*arg)
        or is_a<Cos>(*arg) or is_a<Tan>(*arg) or is_a<Cot>(*arg)
        or is_a<Sec>(*arg) or is_a<Csc>(*arg) or is_a<Sinh>(*arg)
        or is_a<Cosh>(*arg) or is_a<Tanh>(*arg) or is_a<Coth>(*arg)
        or is_a<Sech>(*arg) or is_a<Csch>(*arg)) {
        const OneArgFunction &func = down_cast<const OneArgFunction &>(*arg);
        return func.create(conjugate(func.get_arg()));
    }
    if (is_a<ATan2>(*arg) or is_a<LowerGamma>(*arg) or is_a<UpperGamma>(*arg)
        or is_a<Beta>(*arg)) {
        const TwoArgFunction &func = down_cast<const TwoArgFunction &>(*arg);
        return func.create(conjugate(func.get_arg1()),
                           conjugate(func.get_arg2()));
    }
    return make_rcp<const Conjugate>(arg);
}
 
bool get_pi_shift(const RCP<const Basic> &arg, const Ptr<RCP<const Number>> &n,
                  const Ptr<RCP<const Basic>> &x)
{
    if (is_a<Add>(*arg)) {
        const Add &s = down_cast<const Add &>(*arg);
        RCP<const Basic> coef = s.get_coef();
        auto size = s.get_dict().size();
        if (size > 1) {
            // arg should be of form `x + n*pi`
            // `n` is an integer
            // `x` is an `Expression`
            bool check_pi = false;
            RCP<const Basic> temp;
            *x = coef;
            for (const auto &p : s.get_dict()) {
                if (eq(*p.first, *pi)
                    and (is_a<Integer>(*p.second)
                         or is_a<Rational>(*p.second))) {
                    check_pi = true;
                    *n = p.second;
                } else {
                    *x = add(mul(p.first, p.second), *x);
                }
            }
            if (check_pi)
                return true;
            else // No term with `pi` found
                return false;
        } else if (size == 1) {
            // arg should be of form `a + n*pi`
            // where `a` is a `Number`.
            auto p = s.get_dict().begin();
            if (eq(*p->first, *pi)
                and (is_a<Integer>(*p->second) or is_a<Rational>(*p->second))) {
                *n = p->second;
                *x = coef;
                return true;
            } else {
                return false;
            }
        } else { // Should never reach here though!
            // Dict of size < 1
            return false;
        }
    } else if (is_a<Mul>(*arg)) {
        // `arg` is of the form `k*pi/12`
        const Mul &s = down_cast<const Mul &>(*arg);
        auto p = s.get_dict().begin();
        // dict should contain symbol `pi` only
        if (s.get_dict().size() == 1 and eq(*p->first, *pi)
            and eq(*p->second, *one)
            and (is_a<Integer>(*s.get_coef())
                 or is_a<Rational>(*s.get_coef()))) {
            *n = s.get_coef();
            *x = zero;
            return true;
        } else {
            return false;
        }
    } else if (eq(*arg, *pi)) {
        *n = one;
        *x = zero;
        return true;
    } else if (eq(*arg, *zero)) {
        *n = zero;
        *x = zero;
        return true;
    } else {
        return false;
    }
}
 
// Return true if arg is of form a+b*pi, with b integer or rational
// with denominator 2. The a may be zero or any expression.
bool trig_has_basic_shift(const RCP<const Basic> &arg)
{
    if (is_a<Add>(*arg)) {
        const Add &s = down_cast<const Add &>(*arg);
        for (const auto &p : s.get_dict()) {
            const auto &temp = mul(p.second, integer(2));
            if (eq(*p.first, *pi)) {
                if (is_a<Integer>(*temp)) {
                    return true;
                }
                if (is_a<Rational>(*temp)) {
                    auto m = down_cast<const Rational &>(*temp)
                                 .as_rational_class();
                    return (m < 0) or (m > 1);
                }
                return false;
            }
        }
        return false;
    } else if (is_a<Mul>(*arg)) {
        // is `arg` of the form `k*pi/2`?
        // dict should contain symbol `pi` only
        // and `k` should be a rational s.t. 0 < k < 1
        const Mul &s = down_cast<const Mul &>(*arg);
        RCP<const Basic> coef = mul(s.get_coef(), integer(2));
        auto p = s.get_dict().begin();
        if (s.get_dict().size() == 1 and eq(*p->first, *pi)
            and eq(*p->second, *one)) {
            if (is_a<Integer>(*coef)) {
                return true;
            }
            if (is_a<Rational>(*coef)) {
                auto m = down_cast<const Rational &>(*coef).as_rational_class();
                return (m < 0) or (m > 1);
            }
            return false;
        } else {
            return false;
        }
    } else if (eq(*arg, *pi)) {
        return true;
    } else if (eq(*arg, *zero)) {
        return true;
    } else {
        return false;
    }
}
 
bool could_extract_minus(const Basic &arg)
{
    if (is_a_Number(arg)) {
        if (down_cast<const Number &>(arg).is_negative()) {
            return true;
        } else if (is_a_Complex(arg)) {
            const ComplexBase &c = down_cast<const ComplexBase &>(arg);
            RCP<const Number> real_part = c.real_part();
            return (real_part->is_negative())
                   or (eq(*real_part, *zero)
                       and c.imaginary_part()->is_negative());
        } else {
            return false;
        }
    } else if (is_a<Mul>(arg)) {
        const Mul &s = down_cast<const Mul &>(arg);
        return could_extract_minus(*s.get_coef());
    } else if (is_a<Add>(arg)) {
        const Add &s = down_cast<const Add &>(arg);
        if (s.get_coef()->is_zero()) {
            map_basic_num d(s.get_dict().begin(), s.get_dict().end());
            return could_extract_minus(*d.begin()->second);
        } else {
            return could_extract_minus(*s.get_coef());
        }
    } else {
        return false;
    }
}
 
bool handle_minus(const RCP<const Basic> &arg,
                  const Ptr<RCP<const Basic>> &rarg)
{
    if (is_a<Mul>(*arg)) {
        const Mul &s = down_cast<const Mul &>(*arg);
        // Check for -Add instances to transform -(-x + 2*y) to (x - 2*y)
        if (s.get_coef()->is_minus_one() && s.get_dict().size() == 1
            && eq(*s.get_dict().begin()->second, *one)) {
            return not handle_minus(mul(minus_one, arg), rarg);
        } else if (could_extract_minus(*s.get_coef())) {
            *rarg = mul(minus_one, arg);
            return true;
        }
    } else if (is_a<Add>(*arg)) {
        if (could_extract_minus(*arg)) {
            const Add &s = down_cast<const Add &>(*arg);
            umap_basic_num d = s.get_dict();
            for (auto &p : d) {
                p.second = p.second->mul(*minus_one);
            }
            *rarg = Add::from_dict(s.get_coef()->mul(*minus_one), std::move(d));
            return true;
        }
    } else if (could_extract_minus(*arg)) {
        *rarg = mul(minus_one, arg);
        return true;
    }
    *rarg = arg;
    return false;
}
 
// \return true if conjugate has to be returned finally else false
bool trig_simplify(const RCP<const Basic> &arg, unsigned period, bool odd,
                   bool conj_odd, // input
                   const Ptr<RCP<const Basic>> &rarg, int &index,
                   int &sign) // output
{
    bool check;
    RCP<const Number> n;
    RCP<const Basic> r;
    RCP<const Basic> ret_arg;
    check = get_pi_shift(arg, outArg(n), outArg(r));
    if (check) {
        RCP<const Number> t = mulnum(n, integer(12));
        sign = 1;
        if (is_a<Integer>(*t)) {
            int m = numeric_cast<int>(
                mod_f(down_cast<const Integer &>(*t), *integer(12 * period))
                    ->as_int());
            if (eq(*r, *zero)) {
                index = m;
                *rarg = zero;
                return false;
            } else if (m == 0) {
                index = 0;
                bool b = handle_minus(r, outArg(ret_arg));
                *rarg = ret_arg;
                if (odd and b)
                    sign = -1;
                return false;
            }
        }
 
        rational_class m;
        if (is_a<Integer>(*n)) {
            // 2*pi periodic => f(r + pi * n) = f(r - pi * n)
            m = mp_abs(down_cast<const Integer &>(*n).as_integer_class());
            m /= period;
        } else {
            SYMENGINE_ASSERT(is_a<Rational>(*n));
            m = down_cast<const Rational &>(*n).as_rational_class() / period;
            integer_class t;
#if SYMENGINE_INTEGER_CLASS != SYMENGINE_BOOSTMP
            mp_fdiv_r(t, get_num(m), get_den(m));
            get_num(m) = t;
#else
            integer_class quo;
            mp_fdiv_qr(quo, t, get_num(m), get_den(m));
            m -= rational_class(quo);
#endif
            // m = a / b => m = (a % b / b)
        }
        // Now, arg = r + 2 * pi * m  where 0 <= m < 1
        m *= 2 * period;
        // Now, arg = r + pi * m / 2  where 0 <= m < 4
        if (m >= 2 and m < 3) {
            sign = -1;
            r = add(r, mul(pi, Rational::from_mpq((m - 2) / 2)));
            bool b = handle_minus(r, outArg(ret_arg));
            *rarg = ret_arg;
            if (odd and b)
                sign = -1 * sign;
            return false;
        } else if (m >= 1) {
            if (m < 2) {
                // 1 <= m < 2
                sign = 1;
                r = add(r, mul(pi, Rational::from_mpq((m - 1) / 2)));
            } else {
                // 3 <= m < 4
                sign = -1;
                r = add(r, mul(pi, Rational::from_mpq((m - 3) / 2)));
            }
            bool b = handle_minus(r, outArg(ret_arg));
            *rarg = ret_arg;
            if (not b and conj_odd)
                sign = -sign;
            return true;
        } else {
            *rarg = add(r, mul(pi, Rational::from_mpq(m / 2)));
            index = -1;
            return false;
        }
    } else {
        bool b = handle_minus(arg, outArg(ret_arg));
        *rarg = ret_arg;
        index = -1;
        if (odd and b)
            sign = -1;
        else
            sign = 1;
        return false;
    }
}
 
bool inverse_lookup(const umap_basic_basic &d, const RCP<const Basic> &t,
                    const Ptr<RCP<const Basic>> &index)
{
    auto it = d.find(t);
    if (it == d.end()) {
        // Not found in lookup
        return false;
    } else {
        *index = (it->second);
        return true;
    }
}
 
Sign::Sign(const RCP<const Basic> &arg) : OneArgFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool Sign::is_canonical(const RCP<const Basic> &arg) const
{
    if (is_a_Number(*arg)) {
        if (eq(*arg, *ComplexInf)) {
            return true;
        }
        return false;
    }
    if (is_a<Constant>(*arg)) {
        return false;
    }
    if (is_a<Sign>(*arg)) {
        return false;
    }
    if (is_a<Mul>(*arg)) {
        if (neq(*down_cast<const Mul &>(*arg).get_coef(), *one)
            and neq(*down_cast<const Mul &>(*arg).get_coef(), *minus_one)) {
            return false;
        }
    }
    return true;
}
 
RCP<const Basic> Sign::create(const RCP<const Basic> &arg) const
{
    return sign(arg);
}
 
RCP<const Basic> sign(const RCP<const Basic> &arg)
{
    if (is_a_Number(*arg)) {
        if (is_a<NaN>(*arg)) {
            return Nan;
        }
        if (down_cast<const Number &>(*arg).is_zero()) {
            return zero;
        }
        if (down_cast<const Number &>(*arg).is_positive()) {
            return one;
        }
        if (down_cast<const Number &>(*arg).is_negative()) {
            return minus_one;
        }
        if (is_a_Complex(*arg)
            and down_cast<const ComplexBase &>(*arg).is_re_zero()) {
            RCP<const Number> r
                = down_cast<const ComplexBase &>(*arg).imaginary_part();
            if (down_cast<const Number &>(*r).is_positive()) {
                return I;
            }
            if (down_cast<const Number &>(*r).is_negative()) {
                return mul(minus_one, I);
            }
        }
    }
    if (is_a<Constant>(*arg)) {
        if (eq(*arg, *pi) or eq(*arg, *E) or eq(*arg, *EulerGamma)
            or eq(*arg, *Catalan) or eq(*arg, *GoldenRatio))
            return one;
    }
    if (is_a<Sign>(*arg)) {
        return arg;
    }
    if (is_a<Mul>(*arg)) {
        RCP<const Basic> s = sign(down_cast<const Mul &>(*arg).get_coef());
        map_basic_basic dict = down_cast<const Mul &>(*arg).get_dict();
        return mul(s,
                   make_rcp<const Sign>(Mul::from_dict(one, std::move(dict))));
    }
    if (is_a<Pow>(*arg)) {
        RCP<const Pow> pow_arg = rcp_static_cast<const Pow>(arg);
        RCP<const Basic> s = sign(pow_arg->get_base());
        if (not is_a<Sign>(*s) and not eq(*s, *pow_arg->get_base())) {
            return sign(pow(s, pow_arg->get_exp()));
        }
    }
    return make_rcp<const Sign>(arg);
}
 
Floor::Floor(const RCP<const Basic> &arg) : OneArgFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool Floor::is_canonical(const RCP<const Basic> &arg) const
{
    if (is_a_Number(*arg)) {
        return false;
    }
    if (is_a<Constant>(*arg)) {
        return false;
    }
    if (is_a<Floor>(*arg)) {
        return false;
    }
    if (is_a<Ceiling>(*arg)) {
        return false;
    }
    if (is_a<Truncate>(*arg)) {
        return false;
    }
    if (is_a<BooleanAtom>(*arg) or is_a_Relational(*arg)) {
        return false;
    }
    if (is_a<Add>(*arg)) {
        RCP<const Number> s = down_cast<const Add &>(*arg).get_coef();
        if (neq(*zero, *s) and is_a<Integer>(*s)) {
            return false;
        }
    }
    return true;
}
 
RCP<const Basic> Floor::create(const RCP<const Basic> &arg) const
{
    return floor(arg);
}
 
RCP<const Basic> floor(const RCP<const Basic> &arg)
{
    if (is_a_Number(*arg)) {
        if (down_cast<const Number &>(*arg).is_exact()) {
            if (is_a<Rational>(*arg)) {
                const Rational &s = down_cast<const Rational &>(*arg);
                integer_class quotient;
                mp_fdiv_q(quotient, SymEngine::get_num(s.as_rational_class()),
                          SymEngine::get_den(s.as_rational_class()));
                return integer(std::move(quotient));
            }
            return arg;
        }
        RCP<const Number> _arg = rcp_static_cast<const Number>(arg);
        return _arg->get_eval().floor(*_arg);
    }
    if (is_a<Constant>(*arg)) {
        if (eq(*arg, *pi)) {
            return integer(3);
        }
        if (eq(*arg, *E)) {
            return integer(2);
        }
        if (eq(*arg, *GoldenRatio)) {
            return integer(1);
        }
        if (eq(*arg, *Catalan) or eq(*arg, *EulerGamma)) {
            return integer(0);
        }
    }
    if (is_a<Floor>(*arg)) {
        return arg;
    }
    if (is_a<Ceiling>(*arg)) {
        return arg;
    }
    if (is_a<Truncate>(*arg)) {
        return arg;
    }
    if (is_a<BooleanAtom>(*arg) or is_a_Relational(*arg)) {
        throw SymEngineException(
            "Boolean objects not allowed in this context.");
    }
    if (is_a<Add>(*arg)) {
        RCP<const Number> s = down_cast<const Add &>(*arg).get_coef();
        umap_basic_num d = down_cast<const Add &>(*arg).get_dict();
        if (is_a<Integer>(*s)
            and not down_cast<const Integer &>(*s).is_zero()) {
            return add(s, floor(Add::from_dict(zero, std::move(d))));
        }
    }
    return make_rcp<const Floor>(arg);
}
 
Ceiling::Ceiling(const RCP<const Basic> &arg) : OneArgFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool Ceiling::is_canonical(const RCP<const Basic> &arg) const
{
    if (is_a_Number(*arg)) {
        return false;
    }
    if (is_a<Constant>(*arg)) {
        return false;
    }
    if (is_a<Floor>(*arg)) {
        return false;
    }
    if (is_a<Ceiling>(*arg)) {
        return false;
    }
    if (is_a<Truncate>(*arg)) {
        return false;
    }
    if (is_a<BooleanAtom>(*arg) or is_a_Relational(*arg)) {
        return false;
    }
    if (is_a<Add>(*arg)) {
        RCP<const Number> s = down_cast<const Add &>(*arg).get_coef();
        if (neq(*zero, *s) and is_a<Integer>(*s)) {
            return false;
        }
    }
    return true;
}
 
RCP<const Basic> Ceiling::create(const RCP<const Basic> &arg) const
{
    return ceiling(arg);
}
 
RCP<const Basic> ceiling(const RCP<const Basic> &arg)
{
    if (is_a_Number(*arg)) {
        if (down_cast<const Number &>(*arg).is_exact()) {
            if (is_a<Rational>(*arg)) {
                const Rational &s = down_cast<const Rational &>(*arg);
                integer_class quotient;
                mp_cdiv_q(quotient, SymEngine::get_num(s.as_rational_class()),
                          SymEngine::get_den(s.as_rational_class()));
                return integer(std::move(quotient));
            }
            return arg;
        }
        RCP<const Number> _arg = rcp_static_cast<const Number>(arg);
        return _arg->get_eval().ceiling(*_arg);
    }
    if (is_a<Constant>(*arg)) {
        if (eq(*arg, *pi)) {
            return integer(4);
        }
        if (eq(*arg, *E)) {
            return integer(3);
        }
        if (eq(*arg, *GoldenRatio)) {
            return integer(2);
        }
        if (eq(*arg, *Catalan) or eq(*arg, *EulerGamma)) {
            return integer(1);
        }
    }
    if (is_a<Floor>(*arg)) {
        return arg;
    }
    if (is_a<Ceiling>(*arg)) {
        return arg;
    }
    if (is_a<Truncate>(*arg)) {
        return arg;
    }
    if (is_a<BooleanAtom>(*arg) or is_a_Relational(*arg)) {
        throw SymEngineException(
            "Boolean objects not allowed in this context.");
    }
    if (is_a<Add>(*arg)) {
        RCP<const Number> s = down_cast<const Add &>(*arg).get_coef();
        umap_basic_num d = down_cast<const Add &>(*arg).get_dict();
        if (is_a<Integer>(*s)) {
            return add(
                s, make_rcp<const Ceiling>(Add::from_dict(zero, std::move(d))));
        }
    }
    return make_rcp<const Ceiling>(arg);
}
 
Truncate::Truncate(const RCP<const Basic> &arg) : OneArgFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool Truncate::is_canonical(const RCP<const Basic> &arg) const
{
    if (is_a_Number(*arg)) {
        return false;
    }
    if (is_a<Constant>(*arg)) {
        return false;
    }
    if (is_a<Floor>(*arg)) {
        return false;
    }
    if (is_a<Ceiling>(*arg)) {
        return false;
    }
    if (is_a<Truncate>(*arg)) {
        return false;
    }
    if (is_a<BooleanAtom>(*arg) or is_a_Relational(*arg)) {
        return false;
    }
    if (is_a<Add>(*arg)) {
        RCP<const Number> s = down_cast<const Add &>(*arg).get_coef();
        if (neq(*zero, *s) and is_a<Integer>(*s)) {
            return false;
        }
    }
    return true;
}
 
RCP<const Basic> Truncate::create(const RCP<const Basic> &arg) const
{
    return truncate(arg);
}
 
RCP<const Basic> truncate(const RCP<const Basic> &arg)
{
    if (is_a_Number(*arg)) {
        if (down_cast<const Number &>(*arg).is_exact()) {
            if (is_a<Rational>(*arg)) {
                const Rational &s = down_cast<const Rational &>(*arg);
                integer_class quotient;
                mp_tdiv_q(quotient, SymEngine::get_num(s.as_rational_class()),
                          SymEngine::get_den(s.as_rational_class()));
                return integer(std::move(quotient));
            }
            return arg;
        }
        RCP<const Number> _arg = rcp_static_cast<const Number>(arg);
        return _arg->get_eval().truncate(*_arg);
    }
    if (is_a<Constant>(*arg)) {
        if (eq(*arg, *pi)) {
            return integer(3);
        }
        if (eq(*arg, *E)) {
            return integer(2);
        }
        if (eq(*arg, *GoldenRatio)) {
            return integer(1);
        }
        if (eq(*arg, *Catalan) or eq(*arg, *EulerGamma)) {
            return integer(0);
        }
    }
    if (is_a<Floor>(*arg)) {
        return arg;
    }
    if (is_a<Ceiling>(*arg)) {
        return arg;
    }
    if (is_a<Truncate>(*arg)) {
        return arg;
    }
    if (is_a<BooleanAtom>(*arg) or is_a_Relational(*arg)) {
        throw SymEngineException(
            "Boolean objects not allowed in this context.");
    }
    if (is_a<Add>(*arg)) {
        RCP<const Number> s = down_cast<const Add &>(*arg).get_coef();
        umap_basic_num d = down_cast<const Add &>(*arg).get_dict();
        if (is_a<Integer>(*s)) {
            return add(s, make_rcp<const Truncate>(
                              Add::from_dict(zero, std::move(d))));
        }
    }
    return make_rcp<const Truncate>(arg);
}
 
Sin::Sin(const RCP<const Basic> &arg) : TrigFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool Sin::is_canonical(const RCP<const Basic> &arg) const
{
    // e.g. sin(0)
    if (is_a<Integer>(*arg) and down_cast<const Integer &>(*arg).is_zero())
        return false;
    // e.g sin(7*pi/2+y)
    if (trig_has_basic_shift(arg)) {
        return false;
    }
    if (is_a_Number(*arg) and not down_cast<const Number &>(*arg).is_exact()) {
        return false;
    }
    return true;
}
 
RCP<const Basic> sin(const RCP<const Basic> &arg)
{
    if (eq(*arg, *zero))
        return zero;
    if (is_a_Number(*arg) and not down_cast<const Number &>(*arg).is_exact()) {
        return down_cast<const Number &>(*arg).get_eval().sin(*arg);
    }
 
    if (is_a<ASin>(*arg)) {
        return down_cast<const ASin &>(*arg).get_arg();
    } else if (is_a<ACsc>(*arg)) {
        return div(one, down_cast<const ACsc &>(*arg).get_arg());
    }
 
    RCP<const Basic> ret_arg;
    int index, sign;
    bool conjugate = trig_simplify(arg, 2, true, false,           // input
                                   outArg(ret_arg), index, sign); // output
 
    if (conjugate) {
        // cos has to be returned
        if (sign == 1) {
            return cos(ret_arg);
        } else {
            return mul(minus_one, cos(ret_arg));
        }
    } else {
        if (eq(*ret_arg, *zero)) {
            return mul(integer(sign), sin_table()[index]);
        } else {
            // If ret_arg is the same as arg, a `Sin` instance is returned
            // Or else `sin` is called again.
            if (sign == 1) {
                if (neq(*ret_arg, *arg)) {
                    return sin(ret_arg);
                } else {
                    return make_rcp<const Sin>(arg);
                }
            } else {
                return mul(minus_one, sin(ret_arg));
            }
        }
    }
}
 
/* ---------------------------- */
 
Cos::Cos(const RCP<const Basic> &arg) : TrigFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool Cos::is_canonical(const RCP<const Basic> &arg) const
{
    // e.g. cos(0)
    if (is_a<Integer>(*arg) and down_cast<const Integer &>(*arg).is_zero())
        return false;
    // e.g cos(k*pi/2)
    if (trig_has_basic_shift(arg)) {
        return false;
    }
    if (is_a_Number(*arg) and not down_cast<const Number &>(*arg).is_exact()) {
        return false;
    }
    return true;
}
 
RCP<const Basic> cos(const RCP<const Basic> &arg)
{
    if (eq(*arg, *zero))
        return one;
    if (is_a_Number(*arg) and not down_cast<const Number &>(*arg).is_exact()) {
        return down_cast<const Number &>(*arg).get_eval().cos(*arg);
    }
 
    if (is_a<ACos>(*arg)) {
        return down_cast<const ACos &>(*arg).get_arg();
    } else if (is_a<ASec>(*arg)) {
        return div(one, down_cast<const ASec &>(*arg).get_arg());
    }
 
    RCP<const Basic> ret_arg;
    int index, sign;
    bool conjugate = trig_simplify(arg, 2, false, true,           // input
                                   outArg(ret_arg), index, sign); // output
 
    if (conjugate) {
        // sin has to be returned
        if (sign == 1) {
            return sin(ret_arg);
        } else {
            return mul(minus_one, sin(ret_arg));
        }
    } else {
        if (eq(*ret_arg, *zero)) {
            return mul(integer(sign), sin_table()[(index + 6) % 24]);
        } else {
            if (sign == 1) {
                if (neq(*ret_arg, *arg)) {
                    return cos(ret_arg);
                } else {
                    return make_rcp<const Cos>(ret_arg);
                }
            } else {
                return mul(minus_one, cos(ret_arg));
            }
        }
    }
}
 
/* ---------------------------- */
 
Tan::Tan(const RCP<const Basic> &arg) : TrigFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool Tan::is_canonical(const RCP<const Basic> &arg) const
{
    if (is_a<Integer>(*arg) and down_cast<const Integer &>(*arg).is_zero())
        return false;
    // e.g tan(k*pi/2)
    if (trig_has_basic_shift(arg)) {
        return false;
    }
    if (is_a_Number(*arg) and not down_cast<const Number &>(*arg).is_exact()) {
        return false;
    }
    return true;
}
 
RCP<const Basic> tan(const RCP<const Basic> &arg)
{
    if (eq(*arg, *zero))
        return zero;
    if (is_a_Number(*arg) and not down_cast<const Number &>(*arg).is_exact()) {
        return down_cast<const Number &>(*arg).get_eval().tan(*arg);
    }
 
    if (is_a<ATan>(*arg)) {
        return down_cast<const ATan &>(*arg).get_arg();
    } else if (is_a<ACot>(*arg)) {
        return div(one, down_cast<const ACot &>(*arg).get_arg());
    }
 
    RCP<const Basic> ret_arg;
    int index, sign;
    bool conjugate = trig_simplify(arg, 1, true, true,            // input
                                   outArg(ret_arg), index, sign); // output
 
    if (conjugate) {
        // cot has to be returned
        if (sign == 1) {
            return cot(ret_arg);
        } else {
            return mul(minus_one, cot(ret_arg));
        }
    } else {
        if (eq(*ret_arg, *zero)) {
            return mul(integer(sign),
                       div(sin_table()[index], sin_table()[(index + 6) % 24]));
        } else {
            if (sign == 1) {
                if (neq(*ret_arg, *arg)) {
                    return tan(ret_arg);
                } else {
                    return make_rcp<const Tan>(ret_arg);
                }
            } else {
                return mul(minus_one, tan(ret_arg));
            }
        }
    }
}
 
/* ---------------------------- */
 
Cot::Cot(const RCP<const Basic> &arg) : TrigFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool Cot::is_canonical(const RCP<const Basic> &arg) const
{
    if (is_a<Integer>(*arg) and down_cast<const Integer &>(*arg).is_zero())
        return false;
    // e.g cot(k*pi/2)
    if (trig_has_basic_shift(arg)) {
        return false;
    }
    if (is_a_Number(*arg) and not down_cast<const Number &>(*arg).is_exact()) {
        return false;
    }
    return true;
}
 
RCP<const Basic> cot(const RCP<const Basic> &arg)
{
    if (is_a_Number(*arg) and not down_cast<const Number &>(*arg).is_exact()) {
        return down_cast<const Number &>(*arg).get_eval().cot(*arg);
    }
 
    if (is_a<ACot>(*arg)) {
        return down_cast<const ACot &>(*arg).get_arg();
    } else if (is_a<ATan>(*arg)) {
        return div(one, down_cast<const ATan &>(*arg).get_arg());
    }
 
    RCP<const Basic> ret_arg;
    int index, sign;
    bool conjugate = trig_simplify(arg, 1, true, true,            // input
                                   outArg(ret_arg), index, sign); // output
 
    if (conjugate) {
        // tan has to be returned
        if (sign == 1) {
            return tan(ret_arg);
        } else {
            return mul(minus_one, tan(ret_arg));
        }
    } else {
        if (eq(*ret_arg, *zero)) {
            return mul(integer(sign),
                       div(sin_table()[(index + 6) % 24], sin_table()[index]));
        } else {
            if (sign == 1) {
                if (neq(*ret_arg, *arg)) {
                    return cot(ret_arg);
                } else {
                    return make_rcp<const Cot>(ret_arg);
                }
            } else {
                return mul(minus_one, cot(ret_arg));
            }
        }
    }
}
 
/* ---------------------------- */
 
Csc::Csc(const RCP<const Basic> &arg) : TrigFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool Csc::is_canonical(const RCP<const Basic> &arg) const
{
    // e.g. Csc(0)
    if (is_a<Integer>(*arg) and down_cast<const Integer &>(*arg).is_zero())
        return false;
    // e.g csc(k*pi/2)
    if (trig_has_basic_shift(arg)) {
        return false;
    }
    if (is_a_Number(*arg) and not down_cast<const Number &>(*arg).is_exact()) {
        return false;
    }
    return true;
}
 
RCP<const Basic> csc(const RCP<const Basic> &arg)
{
    if (is_a_Number(*arg) and not down_cast<const Number &>(*arg).is_exact()) {
        return down_cast<const Number &>(*arg).get_eval().csc(*arg);
    }
 
    if (is_a<ACsc>(*arg)) {
        return down_cast<const ACsc &>(*arg).get_arg();
    } else if (is_a<ASin>(*arg)) {
        return div(one, down_cast<const ASin &>(*arg).get_arg());
    }
 
    RCP<const Basic> ret_arg;
    int index, sign;
    bool conjugate = trig_simplify(arg, 2, true, false,           // input
                                   outArg(ret_arg), index, sign); // output
 
    if (conjugate) {
        // cos has to be returned
        if (sign == 1) {
            return sec(ret_arg);
        } else {
            return mul(minus_one, sec(ret_arg));
        }
    } else {
        if (eq(*ret_arg, *zero)) {
            return mul(integer(sign), div(one, sin_table()[index]));
        } else {
            if (sign == 1) {
                if (neq(*ret_arg, *arg)) {
                    return csc(ret_arg);
                } else {
                    return make_rcp<const Csc>(ret_arg);
                }
            } else {
                return mul(minus_one, csc(ret_arg));
            }
        }
    }
}
 
/* ---------------------------- */
 
Sec::Sec(const RCP<const Basic> &arg) : TrigFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool Sec::is_canonical(const RCP<const Basic> &arg) const
{
    // e.g. Sec(0)
    if (is_a<Integer>(*arg) and down_cast<const Integer &>(*arg).is_zero())
        return false;
    // e.g sec(k*pi/2)
    if (trig_has_basic_shift(arg)) {
        return false;
    }
    if (is_a_Number(*arg) and not down_cast<const Number &>(*arg).is_exact()) {
        return false;
    }
    return true;
}
 
RCP<const Basic> sec(const RCP<const Basic> &arg)
{
    if (is_a_Number(*arg) and not down_cast<const Number &>(*arg).is_exact()) {
        return down_cast<const Number &>(*arg).get_eval().sec(*arg);
    }
 
    if (is_a<ASec>(*arg)) {
        return down_cast<const ASec &>(*arg).get_arg();
    } else if (is_a<ACos>(*arg)) {
        return div(one, down_cast<const ACos &>(*arg).get_arg());
    }
 
    RCP<const Basic> ret_arg;
    int index, sign;
    bool conjugate = trig_simplify(arg, 2, false, true,           // input
                                   outArg(ret_arg), index, sign); // output
 
    if (conjugate) {
        // csc has to be returned
        if (sign == 1) {
            return csc(ret_arg);
        } else {
            return mul(minus_one, csc(ret_arg));
        }
    } else {
        if (eq(*ret_arg, *zero)) {
            return mul(integer(sign), div(one, sin_table()[(index + 6) % 24]));
        } else {
            if (sign == 1) {
                if (neq(*ret_arg, *arg)) {
                    return sec(ret_arg);
                } else {
                    return make_rcp<const Sec>(ret_arg);
                }
            } else {
                return mul(minus_one, sec(ret_arg));
            }
        }
    }
}
/* ---------------------------- */
 
// simplifies trigonometric functions wherever possible
// currently deals with simplifications of type sin(acos())
RCP<const Basic> trig_to_sqrt(const RCP<const Basic> &arg)
{
    RCP<const Basic> i_arg;
 
    if (is_a<Sin>(*arg)) {
        if (is_a<ACos>(*arg->get_args()[0])) {
            i_arg = down_cast<const ACos &>(*(arg->get_args()[0])).get_arg();
            return sqrt(sub(one, pow(i_arg, i2)));
        } else if (is_a<ATan>(*arg->get_args()[0])) {
            i_arg = down_cast<const ATan &>(*(arg->get_args()[0])).get_arg();
            return div(i_arg, sqrt(add(one, pow(i_arg, i2))));
        } else if (is_a<ASec>(*arg->get_args()[0])) {
            i_arg = down_cast<const ASec &>(*(arg->get_args()[0])).get_arg();
            return sqrt(sub(one, pow(i_arg, im2)));
        } else if (is_a<ACot>(*arg->get_args()[0])) {
            i_arg = down_cast<const ACot &>(*(arg->get_args()[0])).get_arg();
            return div(one, mul(i_arg, sqrt(add(one, pow(i_arg, im2)))));
        }
    } else if (is_a<Cos>(*arg)) {
        if (is_a<ASin>(*arg->get_args()[0])) {
            i_arg = down_cast<const ASin &>(*(arg->get_args()[0])).get_arg();
            return sqrt(sub(one, pow(i_arg, i2)));
        } else if (is_a<ATan>(*arg->get_args()[0])) {
            i_arg = down_cast<const ATan &>(*(arg->get_args()[0])).get_arg();
            return div(one, sqrt(add(one, pow(i_arg, i2))));
        } else if (is_a<ACsc>(*arg->get_args()[0])) {
            i_arg = down_cast<const ACsc &>(*(arg->get_args()[0])).get_arg();
            return sqrt(sub(one, pow(i_arg, im2)));
        } else if (is_a<ACot>(*arg->get_args()[0])) {
            i_arg = down_cast<const ACot &>(*(arg->get_args()[0])).get_arg();
            return div(one, sqrt(add(one, pow(i_arg, im2))));
        }
    } else if (is_a<Tan>(*arg)) {
        if (is_a<ASin>(*arg->get_args()[0])) {
            i_arg = down_cast<const ASin &>(*(arg->get_args()[0])).get_arg();
            return div(i_arg, sqrt(sub(one, pow(i_arg, i2))));
        } else if (is_a<ACos>(*arg->get_args()[0])) {
            i_arg = down_cast<const ACos &>(*(arg->get_args()[0])).get_arg();
            return div(sqrt(sub(one, pow(i_arg, i2))), i_arg);
        } else if (is_a<ACsc>(*arg->get_args()[0])) {
            i_arg = down_cast<const ACsc &>(*(arg->get_args()[0])).get_arg();
            return div(one, mul(i_arg, sqrt(sub(one, pow(i_arg, im2)))));
        } else if (is_a<ASec>(*arg->get_args()[0])) {
            i_arg = down_cast<const ASec &>(*(arg->get_args()[0])).get_arg();
            return mul(i_arg, sqrt(sub(one, pow(i_arg, im2))));
        }
    } else if (is_a<Csc>(*arg)) {
        if (is_a<ACos>(*arg->get_args()[0])) {
            i_arg = down_cast<const ACos &>(*(arg->get_args()[0])).get_arg();
            return div(one, sqrt(sub(one, pow(i_arg, i2))));
        } else if (is_a<ATan>(*arg->get_args()[0])) {
            i_arg = down_cast<const ATan &>(*(arg->get_args()[0])).get_arg();
            return div(sqrt(add(one, pow(i_arg, i2))), i_arg);
        } else if (is_a<ASec>(*arg->get_args()[0])) {
            i_arg = down_cast<const ASec &>(*(arg->get_args()[0])).get_arg();
            return div(one, sqrt(sub(one, pow(i_arg, im2))));
        } else if (is_a<ACot>(*arg->get_args()[0])) {
            i_arg = down_cast<const ACot &>(*(arg->get_args()[0])).get_arg();
            return mul(i_arg, sqrt(add(one, pow(i_arg, im2))));
        }
    } else if (is_a<Sec>(*arg)) {
        if (is_a<ASin>(*arg->get_args()[0])) {
            i_arg = down_cast<const ASin &>(*(arg->get_args()[0])).get_arg();
            return div(one, sqrt(sub(one, pow(i_arg, i2))));
        } else if (is_a<ATan>(*arg->get_args()[0])) {
            i_arg = down_cast<const ATan &>(*(arg->get_args()[0])).get_arg();
            return sqrt(add(one, pow(i_arg, i2)));
        } else if (is_a<ACsc>(*arg->get_args()[0])) {
            i_arg = down_cast<const ACsc &>(*(arg->get_args()[0])).get_arg();
            return div(one, sqrt(sub(one, pow(i_arg, im2))));
        } else if (is_a<ACot>(*arg->get_args()[0])) {
            i_arg = down_cast<const ACot &>(*(arg->get_args()[0])).get_arg();
            return sqrt(add(one, pow(i_arg, im2)));
        }
    } else if (is_a<Cot>(*arg)) {
        if (is_a<ASin>(*arg->get_args()[0])) {
            i_arg = down_cast<const ASin &>(*(arg->get_args()[0])).get_arg();
            return div(sqrt(sub(one, pow(i_arg, i2))), i_arg);
        } else if (is_a<ACos>(*arg->get_args()[0])) {
            i_arg = down_cast<const ACos &>(*(arg->get_args()[0])).get_arg();
            return div(i_arg, sqrt(sub(one, pow(i_arg, i2))));
        } else if (is_a<ACsc>(*arg->get_args()[0])) {
            i_arg = down_cast<const ACsc &>(*(arg->get_args()[0])).get_arg();
            return mul(i_arg, sqrt(sub(one, pow(i_arg, im2))));
        } else if (is_a<ASec>(*arg->get_args()[0])) {
            i_arg = down_cast<const ASec &>(*(arg->get_args()[0])).get_arg();
            return div(one, mul(i_arg, sqrt(sub(one, pow(i_arg, im2)))));
        }
    }
 
    return arg;
}
 
/* ---------------------------- */
ASin::ASin(const RCP<const Basic> &arg) : InverseTrigFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool ASin::is_canonical(const RCP<const Basic> &arg) const
{
    if (eq(*arg, *zero) or eq(*arg, *one) or eq(*arg, *minus_one))
        return false;
    RCP<const Basic> index;
    if (inverse_lookup(inverse_cst(), get_arg(), outArg(index))) {
        return false;
    }
    if (is_a_Number(*arg) and not down_cast<const Number &>(*arg).is_exact()) {
        return false;
    }
    return true;
}
 
RCP<const Basic> asin(const RCP<const Basic> &arg)
{
    if (eq(*arg, *zero))
        return zero;
    else if (eq(*arg, *one))
        return div(pi, i2);
    else if (eq(*arg, *minus_one))
        return mul(minus_one, div(pi, i2));
    else if (is_a_Number(*arg)
             and not down_cast<const Number &>(*arg).is_exact()) {
        return down_cast<const Number &>(*arg).get_eval().asin(*arg);
    }
 
    RCP<const Basic> index;
    bool b = inverse_lookup(inverse_cst(), arg, outArg(index));
    if (b) {
        return div(pi, index);
    } else {
        return make_rcp<const ASin>(arg);
    }
}
 
ACos::ACos(const RCP<const Basic> &arg) : InverseTrigFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool ACos::is_canonical(const RCP<const Basic> &arg) const
{
    if (eq(*arg, *zero) or eq(*arg, *one) or eq(*arg, *minus_one))
        return false;
    RCP<const Basic> index;
    if (inverse_lookup(inverse_cst(), get_arg(), outArg(index))) {
        return false;
    }
    if (is_a_Number(*arg) and not down_cast<const Number &>(*arg).is_exact()) {
        return false;
    }
    return true;
}
 
RCP<const Basic> acos(const RCP<const Basic> &arg)
{
    if (eq(*arg, *zero))
        return div(pi, i2);
    else if (eq(*arg, *one))
        return zero;
    else if (eq(*arg, *minus_one))
        return pi;
    else if (is_a_Number(*arg)
             and not down_cast<const Number &>(*arg).is_exact()) {
        return down_cast<const Number &>(*arg).get_eval().acos(*arg);
    }
 
    RCP<const Basic> index;
    bool b = inverse_lookup(inverse_cst(), arg, outArg(index));
    if (b) {
        return sub(div(pi, i2), div(pi, index));
    } else {
        return make_rcp<const ACos>(arg);
    }
}
 
ASec::ASec(const RCP<const Basic> &arg) : InverseTrigFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool ASec::is_canonical(const RCP<const Basic> &arg) const
{
    if (eq(*arg, *one) or eq(*arg, *minus_one))
        return false;
    RCP<const Basic> index;
    if (inverse_lookup(inverse_cst(), div(one, get_arg()), outArg(index))) {
        return false;
    } else if (is_a_Number(*arg)
               and not down_cast<const Number &>(*arg).is_exact()) {
        return false;
    }
    return true;
}
 
RCP<const Basic> asec(const RCP<const Basic> &arg)
{
    if (eq(*arg, *one))
        return zero;
    else if (eq(*arg, *minus_one))
        return pi;
    else if (is_a_Number(*arg)
             and not down_cast<const Number &>(*arg).is_exact()) {
        return down_cast<const Number &>(*arg).get_eval().asec(*arg);
    }
 
    RCP<const Basic> index;
    bool b = inverse_lookup(inverse_cst(), div(one, arg), outArg(index));
    if (b) {
        return sub(div(pi, i2), div(pi, index));
    } else {
        return make_rcp<const ASec>(arg);
    }
}
 
ACsc::ACsc(const RCP<const Basic> &arg) : InverseTrigFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool ACsc::is_canonical(const RCP<const Basic> &arg) const
{
    if (eq(*arg, *one) or eq(*arg, *minus_one))
        return false;
    RCP<const Basic> index;
    if (inverse_lookup(inverse_cst(), div(one, arg), outArg(index))) {
        return false;
    }
    if (is_a_Number(*arg) and not down_cast<const Number &>(*arg).is_exact()) {
        return false;
    }
    return true;
}
 
RCP<const Basic> acsc(const RCP<const Basic> &arg)
{
    if (eq(*arg, *one))
        return div(pi, i2);
    else if (eq(*arg, *minus_one))
        return div(pi, im2);
    else if (is_a_Number(*arg)
             and not down_cast<const Number &>(*arg).is_exact()) {
        return down_cast<const Number &>(*arg).get_eval().acsc(*arg);
    }
 
    RCP<const Basic> index;
    bool b = inverse_lookup(inverse_cst(), div(one, arg), outArg(index));
    if (b) {
        return div(pi, index);
    } else {
        return make_rcp<const ACsc>(arg);
    }
}
 
ATan::ATan(const RCP<const Basic> &arg) : InverseTrigFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool ATan::is_canonical(const RCP<const Basic> &arg) const
{
    if (eq(*arg, *zero) or eq(*arg, *one) or eq(*arg, *minus_one))
        return false;
    RCP<const Basic> index;
    if (inverse_lookup(inverse_tct(), get_arg(), outArg(index))) {
        return false;
    }
    if (is_a_Number(*arg) and not down_cast<const Number &>(*arg).is_exact()) {
        return false;
    }
    return true;
}
 
RCP<const Basic> atan(const RCP<const Basic> &arg)
{
    if (eq(*arg, *zero))
        return zero;
    else if (eq(*arg, *one))
        return div(pi, mul(i2, i2));
    else if (eq(*arg, *minus_one))
        return mul(minus_one, div(pi, mul(i2, i2)));
    else if (is_a_Number(*arg)
             and not down_cast<const Number &>(*arg).is_exact()) {
        return down_cast<const Number &>(*arg).get_eval().atan(*arg);
    }
 
    RCP<const Basic> index;
    bool b = inverse_lookup(inverse_tct(), arg, outArg(index));
    if (b) {
        return div(pi, index);
    } else {
        return make_rcp<const ATan>(arg);
    }
}
 
ACot::ACot(const RCP<const Basic> &arg) : InverseTrigFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool ACot::is_canonical(const RCP<const Basic> &arg) const
{
    if (eq(*arg, *zero) or eq(*arg, *one) or eq(*arg, *minus_one))
        return false;
    RCP<const Basic> index;
    if (inverse_lookup(inverse_tct(), arg, outArg(index))) {
        return false;
    }
    if (is_a_Number(*arg) and not down_cast<const Number &>(*arg).is_exact()) {
        return false;
    }
    return true;
}
 
RCP<const Basic> acot(const RCP<const Basic> &arg)
{
    if (eq(*arg, *zero))
        return div(pi, i2);
    else if (eq(*arg, *one))
        return div(pi, mul(i2, i2));
    else if (eq(*arg, *minus_one))
        return mul(i3, div(pi, mul(i2, i2)));
    else if (is_a_Number(*arg)
             and not down_cast<const Number &>(*arg).is_exact()) {
        return down_cast<const Number &>(*arg).get_eval().acot(*arg);
    }
 
    RCP<const Basic> index;
    bool b = inverse_lookup(inverse_tct(), arg, outArg(index));
    if (b) {
        return sub(div(pi, i2), div(pi, index));
    } else {
        return make_rcp<const ACot>(arg);
    }
}
 
ATan2::ATan2(const RCP<const Basic> &num, const RCP<const Basic> &den)
    : TwoArgFunction(num, den)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(num, den))
}
 
bool ATan2::is_canonical(const RCP<const Basic> &num,
                         const RCP<const Basic> &den) const
{
    if (eq(*num, *den) or eq(*num, *mul(minus_one, den)))
        return false;
    if (eq(*num, *zero)) {
        return not(is_a_Number(*den) or eq(*den, *zero));
    }
    if (eq(*den, *zero)) {
        return not is_a_Number(*num);
    }
    RCP<const Basic> index;
    bool b = inverse_lookup(inverse_tct(), div(num, den), outArg(index));
    if (b)
        return false;
    else
        return true;
}
 
RCP<const Basic> ATan2::create(const RCP<const Basic> &a,
                               const RCP<const Basic> &b) const
{
    return atan2(a, b);
}
 
RCP<const Basic> atan2(const RCP<const Basic> &num, const RCP<const Basic> &den)
{
    if (eq(*num, *zero)) {
        if (eq(*den, *zero)) {
            return Nan;
        }
        if (is_a_Number(*den)) {
            RCP<const Number> den_new = rcp_static_cast<const Number>(den);
            if (den_new->is_negative())
                return pi;
            else if (den_new->is_positive())
                return zero;
            else {
                return Nan;
            }
        }
        return make_rcp<const ATan2>(num, den);
    } else if (eq(*den, *zero)) {
        if (is_a_Number(*num)) {
            RCP<const Number> num_new = rcp_static_cast<const Number>(num);
            if (num_new->is_negative())
                return div(pi, im2);
            else
                return div(pi, i2);
        }
        return make_rcp<const ATan2>(num, den);
    }
    RCP<const Basic> index;
    bool b = inverse_lookup(inverse_tct(), div(num, den), outArg(index));
    if (b) {
        SYMENGINE_ASSERT(is_a_Number(*index));
        RCP<const Number> index_num = rcp_static_cast<const Number>(index);
        if (index_num->is_positive()) {
            return add(div(pi, index), mul(div(pi, i2), sub(sign(den), one)));
        } else {
            return sub(div(pi, index), mul(div(pi, i2), sub(sign(den), one)));
        }
    } else {
        return make_rcp<const ATan2>(num, den);
    }
}
 
/* ---------------------------- */
 
RCP<const Basic> Sin::create(const RCP<const Basic> &arg) const
{
    return sin(arg);
}
 
RCP<const Basic> Cos::create(const RCP<const Basic> &arg) const
{
    return cos(arg);
}
 
RCP<const Basic> Tan::create(const RCP<const Basic> &arg) const
{
    return tan(arg);
}
 
RCP<const Basic> Cot::create(const RCP<const Basic> &arg) const
{
    return cot(arg);
}
 
RCP<const Basic> Sec::create(const RCP<const Basic> &arg) const
{
    return sec(arg);
}
 
RCP<const Basic> Csc::create(const RCP<const Basic> &arg) const
{
    return csc(arg);
}
 
RCP<const Basic> ASin::create(const RCP<const Basic> &arg) const
{
    return asin(arg);
}
 
RCP<const Basic> ACos::create(const RCP<const Basic> &arg) const
{
    return acos(arg);
}
 
RCP<const Basic> ATan::create(const RCP<const Basic> &arg) const
{
    return atan(arg);
}
 
RCP<const Basic> ACot::create(const RCP<const Basic> &arg) const
{
    return acot(arg);
}
 
RCP<const Basic> ASec::create(const RCP<const Basic> &arg) const
{
    return asec(arg);
}
 
RCP<const Basic> ACsc::create(const RCP<const Basic> &arg) const
{
    return acsc(arg);
}
 
/* ---------------------------- */
 
Log::Log(const RCP<const Basic> &arg) : OneArgFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool Log::is_canonical(const RCP<const Basic> &arg) const
{
    //  log(0)
    if (is_a<Integer>(*arg) and down_cast<const Integer &>(*arg).is_zero())
        return false;
    //  log(1)
    if (is_a<Integer>(*arg) and down_cast<const Integer &>(*arg).is_one())
        return false;
    // log(E)
    if (eq(*arg, *E))
        return false;
 
    if (is_a_Number(*arg) and down_cast<const Number &>(*arg).is_negative())
        return false;
 
    // log(Inf) is also handled here.
    if (is_a_Number(*arg) and not down_cast<const Number &>(*arg).is_exact())
        return false;
 
    // log(3I) should be expanded to log(3) + I*pi/2
    if (is_a<Complex>(*arg) and down_cast<const Complex &>(*arg).is_re_zero())
        return false;
    // log(num/den) = log(num) - log(den)
    if (is_a<Rational>(*arg))
        return false;
    return true;
}
 
RCP<const Basic> Log::create(const RCP<const Basic> &a) const
{
    return log(a);
}
 
RCP<const Basic> log(const RCP<const Basic> &arg)
{
    if (eq(*arg, *zero))
        return ComplexInf;
    if (eq(*arg, *one))
        return zero;
    if (eq(*arg, *E))
        return one;
 
    if (is_a_Number(*arg)) {
        RCP<const Number> _arg = rcp_static_cast<const Number>(arg);
        if (not _arg->is_exact()) {
            return _arg->get_eval().log(*_arg);
        } else if (_arg->is_negative()) {
            return add(log(mul(minus_one, _arg)), mul(pi, I));
        }
    }
 
    if (is_a<Rational>(*arg)) {
        RCP<const Integer> num, den;
        get_num_den(down_cast<const Rational &>(*arg), outArg(num),
                    outArg(den));
        return sub(log(num), log(den));
    }
 
    if (is_a<Complex>(*arg)) {
        RCP<const Complex> _arg = rcp_static_cast<const Complex>(arg);
        if (_arg->is_re_zero()) {
            RCP<const Number> arg_img = _arg->imaginary_part();
            if (arg_img->is_negative()) {
                return sub(log(mul(minus_one, arg_img)),
                           mul(I, div(pi, integer(2))));
            } else if (arg_img->is_zero()) {
                return ComplexInf;
            } else if (arg_img->is_positive()) {
                return add(log(arg_img), mul(I, div(pi, integer(2))));
            }
        }
    }
 
    return make_rcp<const Log>(arg);
}
 
RCP<const Basic> log(const RCP<const Basic> &arg, const RCP<const Basic> &base)
{
    return div(log(arg), log(base));
}
 
LambertW::LambertW(const RCP<const Basic> &arg) : OneArgFunction{arg}
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool LambertW::is_canonical(const RCP<const Basic> &arg) const
{
    if (eq(*arg, *zero))
        return false;
    if (eq(*arg, *E))
        return false;
    if (eq(*arg, *div(neg(one), E)))
        return false;
    if (eq(*arg, *div(log(i2), im2)))
        return false;
    return true;
}
 
RCP<const Basic> LambertW::create(const RCP<const Basic> &arg) const
{
    return lambertw(arg);
}
 
RCP<const Basic> lambertw(const RCP<const Basic> &arg)
{
    if (eq(*arg, *zero))
        return zero;
    if (eq(*arg, *E))
        return one;
    if (eq(*arg, *div(neg(one), E)))
        return minus_one;
    if (eq(*arg, *div(log(i2), im2)))
        return mul(minus_one, log(i2));
    return make_rcp<const LambertW>(arg);
}
 
FunctionSymbol::FunctionSymbol(std::string name, const RCP<const Basic> &arg)
    : MultiArgFunction({arg}), name_{name} {SYMENGINE_ASSIGN_TYPEID()
                                                SYMENGINE_ASSERT(
                                                    is_canonical(get_vec()))}
 
      FunctionSymbol::FunctionSymbol(std::string name, const vec_basic &arg)
    : MultiArgFunction(arg), name_{name}
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(get_vec()))
}
 
bool FunctionSymbol::is_canonical(const vec_basic &arg) const
{
    return true;
}
 
hash_t FunctionSymbol::__hash__() const
{
    hash_t seed = SYMENGINE_FUNCTIONSYMBOL;
    for (const auto &a : get_vec())
        hash_combine<Basic>(seed, *a);
    hash_combine<std::string>(seed, name_);
    return seed;
}
 
bool FunctionSymbol::__eq__(const Basic &o) const
{
    if (is_a<FunctionSymbol>(o)
        and name_ == down_cast<const FunctionSymbol &>(o).name_
        and unified_eq(get_vec(),
                       down_cast<const FunctionSymbol &>(o).get_vec()))
        return true;
    return false;
}
 
int FunctionSymbol::compare(const Basic &o) const
{
    SYMENGINE_ASSERT(is_a<FunctionSymbol>(o))
    const FunctionSymbol &s = down_cast<const FunctionSymbol &>(o);
    if (name_ == s.name_)
        return unified_compare(get_vec(), s.get_vec());
    else
        return name_ < s.name_ ? -1 : 1;
}
 
RCP<const Basic> FunctionSymbol::create(const vec_basic &x) const
{
    return make_rcp<const FunctionSymbol>(name_, x);
}
 
RCP<const Basic> function_symbol(std::string name, const vec_basic &arg)
{
    return make_rcp<const FunctionSymbol>(name, arg);
}
 
RCP<const Basic> function_symbol(std::string name, const RCP<const Basic> &arg)
{
    return make_rcp<const FunctionSymbol>(name, arg);
}
 
FunctionWrapper::FunctionWrapper(std::string name, const RCP<const Basic> &arg)
    : FunctionSymbol(name, arg){SYMENGINE_ASSIGN_TYPEID()}
 
      FunctionWrapper::FunctionWrapper(std::string name, const vec_basic &vec)
    : FunctionSymbol(name, vec){SYMENGINE_ASSIGN_TYPEID()}
 
      /* ---------------------------- */
 
      Derivative::Derivative(const RCP<const Basic> &arg,
                             const multiset_basic &x)
    : arg_{arg}, x_{x}
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg, x))
}
 
bool Derivative::is_canonical(const RCP<const Basic> &arg,
                              const multiset_basic &x) const
{
    // Check that 'x' are Symbols:
    for (const auto &a : x)
        if (not is_a<Symbol>(*a))
            return false;
    if (is_a<FunctionSymbol>(*arg) or is_a<LeviCivita>(*arg)) {
        for (auto &p : x) {
            RCP<const Symbol> s = rcp_static_cast<const Symbol>(p);
            RCP<const MultiArgFunction> f
                = rcp_static_cast<const MultiArgFunction>(arg);
            bool found_s = false;
            // 's' should be one of the args of the function
            // and should not appear anywhere else.
            for (const auto &a : f->get_args()) {
                if (eq(*a, *s)) {
                    if (found_s) {
                        return false;
                    } else {
                        found_s = true;
                    }
                } else if (neq(*a->diff(s), *zero)) {
                    return false;
                }
            }
            if (!found_s) {
                return false;
            }
        }
        return true;
    } else if (is_a<Abs>(*arg)) {
        return true;
    } else if (is_a<FunctionWrapper>(*arg)) {
        return true;
    } else if (is_a<PolyGamma>(*arg) or is_a<Zeta>(*arg)
               or is_a<UpperGamma>(*arg) or is_a<LowerGamma>(*arg)
               or is_a<Dirichlet_eta>(*arg)) {
        bool found = false;
        auto v = arg->get_args();
        for (auto &p : x) {
            if (has_symbol(*v[0], *rcp_static_cast<const Symbol>(p))) {
                found = true;
                break;
            }
        }
        return found;
    } else if (is_a<KroneckerDelta>(*arg)) {
        bool found = false;
        auto v = arg->get_args();
        for (auto &p : x) {
            if (has_symbol(*v[0], *rcp_static_cast<const Symbol>(p))
                or has_symbol(*v[1], *rcp_static_cast<const Symbol>(p))) {
                found = true;
                break;
            }
        }
        return found;
    }
    return false;
}
 
hash_t Derivative::__hash__() const
{
    hash_t seed = SYMENGINE_DERIVATIVE;
    hash_combine<Basic>(seed, *arg_);
    for (auto &p : x_) {
        hash_combine<Basic>(seed, *p);
    }
    return seed;
}
 
bool Derivative::__eq__(const Basic &o) const
{
    if (is_a<Derivative>(o)
        and eq(*arg_, *(down_cast<const Derivative &>(o).arg_))
        and unified_eq(x_, down_cast<const Derivative &>(o).x_))
        return true;
    return false;
}
 
int Derivative::compare(const Basic &o) const
{
    SYMENGINE_ASSERT(is_a<Derivative>(o))
    const Derivative &s = down_cast<const Derivative &>(o);
    int cmp = arg_->__cmp__(*(s.arg_));
    if (cmp != 0)
        return cmp;
    cmp = unified_compare(x_, s.x_);
    return cmp;
}
 
// Subs class
Subs::Subs(const RCP<const Basic> &arg, const map_basic_basic &dict)
    : arg_{arg}, dict_{dict}
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg, dict))
}
 
bool Subs::is_canonical(const RCP<const Basic> &arg,
                        const map_basic_basic &dict) const
{
    if (is_a<Derivative>(*arg)) {
        return true;
    }
    return false;
}
 
hash_t Subs::__hash__() const
{
    hash_t seed = SYMENGINE_SUBS;
    hash_combine<Basic>(seed, *arg_);
    for (const auto &p : dict_) {
        hash_combine<Basic>(seed, *p.first);
        hash_combine<Basic>(seed, *p.second);
    }
    return seed;
}
 
bool Subs::__eq__(const Basic &o) const
{
    if (is_a<Subs>(o) and eq(*arg_, *(down_cast<const Subs &>(o).arg_))
        and unified_eq(dict_, down_cast<const Subs &>(o).dict_))
        return true;
    return false;
}
 
int Subs::compare(const Basic &o) const
{
    SYMENGINE_ASSERT(is_a<Subs>(o))
    const Subs &s = down_cast<const Subs &>(o);
    int cmp = arg_->__cmp__(*(s.arg_));
    if (cmp != 0)
        return cmp;
    cmp = unified_compare(dict_, s.dict_);
    return cmp;
}
 
vec_basic Subs::get_variables() const
{
    vec_basic v;
    for (const auto &p : dict_) {
        v.push_back(p.first);
    }
    return v;
}
 
vec_basic Subs::get_point() const
{
    vec_basic v;
    for (const auto &p : dict_) {
        v.push_back(p.second);
    }
    return v;
}
 
vec_basic Subs::get_args() const
{
    vec_basic v = {arg_};
    for (const auto &p : dict_) {
        v.push_back(p.first);
    }
    for (const auto &p : dict_) {
        v.push_back(p.second);
    }
    return v;
}
 
Sinh::Sinh(const RCP<const Basic> &arg) : HyperbolicFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool Sinh::is_canonical(const RCP<const Basic> &arg) const
{
    if (eq(*arg, *zero))
        return false;
    if (is_a_Number(*arg)) {
        if (down_cast<const Number &>(*arg).is_negative()) {
            return false;
        } else if (not down_cast<const Number &>(*arg).is_exact()) {
            return false;
        }
    }
    if (could_extract_minus(*arg))
        return false;
    return true;
}
 
RCP<const Basic> sinh(const RCP<const Basic> &arg)
{
    if (eq(*arg, *zero))
        return zero;
    if (is_a_Number(*arg)) {
        RCP<const Number> _arg = rcp_static_cast<const Number>(arg);
        if (not _arg->is_exact()) {
            return _arg->get_eval().sinh(*_arg);
        } else if (_arg->is_negative()) {
            return neg(sinh(zero->sub(*_arg)));
        }
    }
    RCP<const Basic> d;
    bool b = handle_minus(arg, outArg(d));
    if (b) {
        return neg(sinh(d));
    }
    return make_rcp<const Sinh>(d);
}
 
Csch::Csch(const RCP<const Basic> &arg) : HyperbolicFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool Csch::is_canonical(const RCP<const Basic> &arg) const
{
    if (eq(*arg, *zero))
        return false;
    if (is_a_Number(*arg)) {
        if (down_cast<const Number &>(*arg).is_negative()) {
            return false;
        } else if (not down_cast<const Number &>(*arg).is_exact()) {
            return false;
        }
    }
    if (could_extract_minus(*arg))
        return false;
    return true;
}
 
RCP<const Basic> csch(const RCP<const Basic> &arg)
{
    if (eq(*arg, *zero)) {
        return ComplexInf;
    }
    if (is_a_Number(*arg)) {
        RCP<const Number> _arg = rcp_static_cast<const Number>(arg);
        if (not _arg->is_exact()) {
            return _arg->get_eval().csch(*_arg);
        } else if (_arg->is_negative()) {
            return neg(csch(zero->sub(*_arg)));
        }
    }
    RCP<const Basic> d;
    bool b = handle_minus(arg, outArg(d));
    if (b) {
        return neg(csch(d));
    }
    return make_rcp<const Csch>(d);
}
 
Cosh::Cosh(const RCP<const Basic> &arg) : HyperbolicFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool Cosh::is_canonical(const RCP<const Basic> &arg) const
{
    if (eq(*arg, *zero))
        return false;
    if (is_a_Number(*arg)) {
        if (down_cast<const Number &>(*arg).is_negative()) {
            return false;
        } else if (not down_cast<const Number &>(*arg).is_exact()) {
            return false;
        }
    }
    if (could_extract_minus(*arg))
        return false;
    return true;
}
 
RCP<const Basic> cosh(const RCP<const Basic> &arg)
{
    if (eq(*arg, *zero))
        return one;
    if (is_a_Number(*arg)) {
        RCP<const Number> _arg = rcp_static_cast<const Number>(arg);
        if (not _arg->is_exact()) {
            return _arg->get_eval().cosh(*_arg);
        } else if (_arg->is_negative()) {
            return cosh(zero->sub(*_arg));
        }
    }
    RCP<const Basic> d;
    handle_minus(arg, outArg(d));
    return make_rcp<const Cosh>(d);
}
 
Sech::Sech(const RCP<const Basic> &arg) : HyperbolicFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool Sech::is_canonical(const RCP<const Basic> &arg) const
{
    if (eq(*arg, *zero))
        return false;
    if (is_a_Number(*arg)) {
        if (down_cast<const Number &>(*arg).is_negative()) {
            return false;
        } else if (not down_cast<const Number &>(*arg).is_exact()) {
            return false;
        }
    }
    if (could_extract_minus(*arg))
        return false;
    return true;
}
 
RCP<const Basic> sech(const RCP<const Basic> &arg)
{
    if (eq(*arg, *zero))
        return one;
    if (is_a_Number(*arg)) {
        RCP<const Number> _arg = rcp_static_cast<const Number>(arg);
        if (not _arg->is_exact()) {
            return _arg->get_eval().sech(*_arg);
        } else if (_arg->is_negative()) {
            return sech(zero->sub(*_arg));
        }
    }
    RCP<const Basic> d;
    handle_minus(arg, outArg(d));
    return make_rcp<const Sech>(d);
}
 
Tanh::Tanh(const RCP<const Basic> &arg) : HyperbolicFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool Tanh::is_canonical(const RCP<const Basic> &arg) const
{
    if (eq(*arg, *zero))
        return false;
    if (is_a_Number(*arg)) {
        if (down_cast<const Number &>(*arg).is_negative()) {
            return false;
        } else if (not down_cast<const Number &>(*arg).is_exact()) {
            return false;
        }
    }
    if (could_extract_minus(*arg))
        return false;
    return true;
}
 
RCP<const Basic> tanh(const RCP<const Basic> &arg)
{
    if (eq(*arg, *zero))
        return zero;
    if (is_a_Number(*arg)) {
        RCP<const Number> _arg = rcp_static_cast<const Number>(arg);
        if (not _arg->is_exact()) {
            return _arg->get_eval().tanh(*_arg);
        } else if (_arg->is_negative()) {
            return neg(tanh(zero->sub(*_arg)));
        }
    }
 
    RCP<const Basic> d;
    bool b = handle_minus(arg, outArg(d));
    if (b) {
        return neg(tanh(d));
    }
    return make_rcp<const Tanh>(d);
}
 
Coth::Coth(const RCP<const Basic> &arg) : HyperbolicFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool Coth::is_canonical(const RCP<const Basic> &arg) const
{
    if (eq(*arg, *zero))
        return false;
    if (is_a_Number(*arg)) {
        if (down_cast<const Number &>(*arg).is_negative()) {
            return false;
        } else if (not down_cast<const Number &>(*arg).is_exact()) {
            return false;
        }
    }
    if (could_extract_minus(*arg))
        return false;
    return true;
}
 
RCP<const Basic> coth(const RCP<const Basic> &arg)
{
    if (eq(*arg, *zero)) {
        return ComplexInf;
    }
    if (is_a_Number(*arg)) {
        RCP<const Number> _arg = rcp_static_cast<const Number>(arg);
        if (not _arg->is_exact()) {
            return _arg->get_eval().coth(*_arg);
        } else if (_arg->is_negative()) {
            return neg(coth(zero->sub(*_arg)));
        }
    }
    RCP<const Basic> d;
    bool b = handle_minus(arg, outArg(d));
    if (b) {
        return neg(coth(d));
    }
    return make_rcp<const Coth>(d);
}
 
ASinh::ASinh(const RCP<const Basic> &arg) : InverseHyperbolicFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool ASinh::is_canonical(const RCP<const Basic> &arg) const
{
    if (eq(*arg, *zero) or eq(*arg, *one) or eq(*arg, *minus_one))
        return false;
    if (is_a_Number(*arg)) {
        if (down_cast<const Number &>(*arg).is_negative()) {
            return false;
        } else if (not down_cast<const Number &>(*arg).is_exact()) {
            return false;
        }
    }
    if (could_extract_minus(*arg))
        return false;
    return true;
}
 
RCP<const Basic> asinh(const RCP<const Basic> &arg)
{
    if (eq(*arg, *zero))
        return zero;
    if (eq(*arg, *one))
        return log(add(one, sq2));
    if (eq(*arg, *minus_one))
        return log(sub(sq2, one));
    if (is_a_Number(*arg)) {
        RCP<const Number> _arg = rcp_static_cast<const Number>(arg);
        if (not _arg->is_exact()) {
            return _arg->get_eval().asinh(*_arg);
        } else if (_arg->is_negative()) {
            return neg(asinh(zero->sub(*_arg)));
        }
    }
    RCP<const Basic> d;
    bool b = handle_minus(arg, outArg(d));
    if (b) {
        return neg(asinh(d));
    }
    return make_rcp<const ASinh>(d);
}
 
ACsch::ACsch(const RCP<const Basic> &arg) : InverseHyperbolicFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool ACsch::is_canonical(const RCP<const Basic> &arg) const
{
    if (eq(*arg, *one) or eq(*arg, *minus_one))
        return false;
    if (is_a_Number(*arg)) {
        if (down_cast<const Number &>(*arg).is_negative()) {
            return false;
        } else if (not down_cast<const Number &>(*arg).is_exact()) {
            return false;
        }
    }
    if (could_extract_minus(*arg))
        return false;
    return true;
}
 
RCP<const Basic> acsch(const RCP<const Basic> &arg)
{
    if (eq(*arg, *one))
        return log(add(one, sq2));
    if (eq(*arg, *minus_one))
        return log(sub(sq2, one));
 
    if (is_a_Number(*arg)) {
        RCP<const Number> _arg = rcp_static_cast<const Number>(arg);
        if (not _arg->is_exact()) {
            return _arg->get_eval().acsch(*_arg);
        }
    }
 
    RCP<const Basic> d;
    bool b = handle_minus(arg, outArg(d));
    if (b) {
        return neg(acsch(d));
    }
    return make_rcp<const ACsch>(d);
}
 
ACosh::ACosh(const RCP<const Basic> &arg) : InverseHyperbolicFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool ACosh::is_canonical(const RCP<const Basic> &arg) const
{
    // TODO: Lookup into a cst table once complex is implemented
    if (eq(*arg, *one))
        return false;
    if (is_a_Number(*arg) and not down_cast<const Number &>(*arg).is_exact()) {
        return false;
    }
    return true;
}
 
RCP<const Basic> acosh(const RCP<const Basic> &arg)
{
    // TODO: Lookup into a cst table once complex is implemented
    if (eq(*arg, *one))
        return zero;
    if (is_a_Number(*arg) and not down_cast<const Number &>(*arg).is_exact()) {
        return down_cast<const Number &>(*arg).get_eval().acosh(*arg);
    }
    return make_rcp<const ACosh>(arg);
}
 
ATanh::ATanh(const RCP<const Basic> &arg) : InverseHyperbolicFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool ATanh::is_canonical(const RCP<const Basic> &arg) const
{
    if (eq(*arg, *zero))
        return false;
    if (is_a_Number(*arg)) {
        if (down_cast<const Number &>(*arg).is_negative()) {
            return false;
        } else if (not down_cast<const Number &>(*arg).is_exact()) {
            return false;
        }
    }
    if (could_extract_minus(*arg))
        return false;
    return true;
}
 
RCP<const Basic> atanh(const RCP<const Basic> &arg)
{
    if (eq(*arg, *zero))
        return zero;
    if (is_a_Number(*arg)) {
        RCP<const Number> _arg = rcp_static_cast<const Number>(arg);
        if (not _arg->is_exact()) {
            return _arg->get_eval().atanh(*_arg);
        } else if (_arg->is_negative()) {
            return neg(atanh(zero->sub(*_arg)));
        }
    }
    RCP<const Basic> d;
    bool b = handle_minus(arg, outArg(d));
    if (b) {
        return neg(atanh(d));
    }
    return make_rcp<const ATanh>(d);
}
 
ACoth::ACoth(const RCP<const Basic> &arg) : InverseHyperbolicFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool ACoth::is_canonical(const RCP<const Basic> &arg) const
{
    if (is_a_Number(*arg)) {
        if (down_cast<const Number &>(*arg).is_negative()) {
            return false;
        } else if (not down_cast<const Number &>(*arg).is_exact()) {
            return false;
        }
    }
    if (could_extract_minus(*arg))
        return false;
    return true;
}
 
RCP<const Basic> acoth(const RCP<const Basic> &arg)
{
    if (is_a_Number(*arg)) {
        RCP<const Number> _arg = rcp_static_cast<const Number>(arg);
        if (not _arg->is_exact()) {
            return _arg->get_eval().acoth(*_arg);
        } else if (_arg->is_negative()) {
            return neg(acoth(zero->sub(*_arg)));
        }
    }
    RCP<const Basic> d;
    bool b = handle_minus(arg, outArg(d));
    if (b) {
        return neg(acoth(d));
    }
    return make_rcp<const ACoth>(d);
}
 
ASech::ASech(const RCP<const Basic> &arg) : InverseHyperbolicFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool ASech::is_canonical(const RCP<const Basic> &arg) const
{
    // TODO: Lookup into a cst table once complex is implemented
    if (eq(*arg, *one))
        return false;
    if (eq(*arg, *zero))
        return false;
    if (is_a_Number(*arg) and not down_cast<const Number &>(*arg).is_exact()) {
        return false;
    }
    return true;
}
 
RCP<const Basic> asech(const RCP<const Basic> &arg)
{
    // TODO: Lookup into a cst table once complex is implemented
    if (eq(*arg, *one))
        return zero;
    if (eq(*arg, *zero))
        return Inf;
    if (is_a_Number(*arg)) {
        RCP<const Number> _arg = rcp_static_cast<const Number>(arg);
        if (not _arg->is_exact()) {
            return _arg->get_eval().asech(*_arg);
        }
    }
    return make_rcp<const ASech>(arg);
}
 
RCP<const Basic> Sinh::create(const RCP<const Basic> &arg) const
{
    return sinh(arg);
}
 
RCP<const Basic> Csch::create(const RCP<const Basic> &arg) const
{
    return csch(arg);
}
 
RCP<const Basic> Cosh::create(const RCP<const Basic> &arg) const
{
    return cosh(arg);
}
 
RCP<const Basic> Sech::create(const RCP<const Basic> &arg) const
{
    return sech(arg);
}
 
RCP<const Basic> Tanh::create(const RCP<const Basic> &arg) const
{
    return tanh(arg);
}
 
RCP<const Basic> Coth::create(const RCP<const Basic> &arg) const
{
    return coth(arg);
}
 
RCP<const Basic> ASinh::create(const RCP<const Basic> &arg) const
{
    return asinh(arg);
}
 
RCP<const Basic> ACsch::create(const RCP<const Basic> &arg) const
{
    return acsch(arg);
}
 
RCP<const Basic> ACosh::create(const RCP<const Basic> &arg) const
{
    return acosh(arg);
}
 
RCP<const Basic> ATanh::create(const RCP<const Basic> &arg) const
{
    return atanh(arg);
}
 
RCP<const Basic> ACoth::create(const RCP<const Basic> &arg) const
{
    return acoth(arg);
}
 
RCP<const Basic> ASech::create(const RCP<const Basic> &arg) const
{
    return asech(arg);
}
 
KroneckerDelta::KroneckerDelta(const RCP<const Basic> &i,
                               const RCP<const Basic> &j)
    : TwoArgFunction(i, j)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(i, j))
}
 
bool KroneckerDelta::is_canonical(const RCP<const Basic> &i,
                                  const RCP<const Basic> &j) const
{
    RCP<const Basic> diff = expand(sub(i, j));
    if (eq(*diff, *zero)) {
        return false;
    } else if (is_a_Number(*diff)) {
        return false;
    } else {
        // TODO: SymPy uses default key sorting to return in order
        return true;
    }
}
 
RCP<const Basic> KroneckerDelta::create(const RCP<const Basic> &a,
                                        const RCP<const Basic> &b) const
{
    return kronecker_delta(a, b);
}
 
RCP<const Basic> kronecker_delta(const RCP<const Basic> &i,
                                 const RCP<const Basic> &j)
{
    // Expand is needed to simplify things like `i-(i+1)` to `-1`
    RCP<const Basic> diff = expand(sub(i, j));
    if (eq(*diff, *zero)) {
        return one;
    } else if (is_a_Number(*diff)) {
        return zero;
    } else {
        // SymPy uses default key sorting to return in order
        return make_rcp<const KroneckerDelta>(i, j);
    }
}
 
bool has_dup(const vec_basic &arg)
{
    map_basic_basic d;
    auto it = d.end();
    for (const auto &p : arg) {
        it = d.find(p);
        if (it == d.end()) {
            insert(d, p, one);
        } else {
            return true;
        }
    }
    return false;
}
 
LeviCivita::LeviCivita(const vec_basic &&arg) : MultiArgFunction(std::move(arg))
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(get_vec()))
}
 
bool LeviCivita::is_canonical(const vec_basic &arg) const
{
    bool are_int = true;
    for (const auto &p : arg) {
        if (not(is_a_Number(*p))) {
            are_int = false;
            break;
        }
    }
    if (are_int) {
        return false;
    } else if (has_dup(arg)) {
        return false;
    } else {
        return true;
    }
}
 
RCP<const Basic> LeviCivita::create(const vec_basic &a) const
{
    return levi_civita(a);
}
 
RCP<const Basic> eval_levicivita(const vec_basic &arg, int len)
{
    int i, j;
    RCP<const Basic> res = one;
    for (i = 0; i < len; i++) {
        for (j = i + 1; j < len; j++) {
            res = mul(sub(arg[j], arg[i]), res);
        }
        res = div(res, factorial(i));
    }
    return res;
}
 
RCP<const Basic> levi_civita(const vec_basic &arg)
{
    bool are_int = true;
    int len = 0;
    for (const auto &p : arg) {
        if (not(is_a_Number(*p))) {
            are_int = false;
            break;
        } else {
            len++;
        }
    }
    if (are_int) {
        return eval_levicivita(arg, len);
    } else if (has_dup(arg)) {
        return zero;
    } else {
        return make_rcp<const LeviCivita>(std::move(arg));
    }
}
 
Zeta::Zeta(const RCP<const Basic> &s, const RCP<const Basic> &a)
    : TwoArgFunction(s, a){SYMENGINE_ASSIGN_TYPEID()
                               SYMENGINE_ASSERT(is_canonical(s, a))}
 
      Zeta::Zeta(const RCP<const Basic> &s)
    : TwoArgFunction(s, one)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(s, one))
}
 
bool Zeta::is_canonical(const RCP<const Basic> &s,
                        const RCP<const Basic> &a) const
{
    if (eq(*s, *zero))
        return false;
    if (eq(*s, *one))
        return false;
    if (is_a<Integer>(*s) and is_a<Integer>(*a)) {
        auto s_ = down_cast<const Integer &>(*s).as_int();
        if (s_ < 0 || s_ % 2 == 0)
            return false;
    }
    return true;
}
 
RCP<const Basic> Zeta::create(const RCP<const Basic> &a,
                              const RCP<const Basic> &b) const
{
    return zeta(a, b);
}
 
RCP<const Basic> zeta(const RCP<const Basic> &s, const RCP<const Basic> &a)
{
    if (is_a_Number(*s)) {
        if (down_cast<const Number &>(*s).is_zero()) {
            return sub(div(one, i2), a);
        } else if (down_cast<const Number &>(*s).is_one()) {
            return infty(0);
        } else if (is_a<Integer>(*s) and is_a<Integer>(*a)) {
            auto s_ = down_cast<const Integer &>(*s).as_int();
            auto a_ = down_cast<const Integer &>(*a).as_int();
            RCP<const Basic> zeta;
            if (s_ < 0) {
                RCP<const Number> res = (s_ % 2 == 0) ? one : minus_one;
                zeta
                    = mulnum(res, divnum(bernoulli(-s_ + 1), integer(-s_ + 1)));
            } else if (s_ % 2 == 0) {
                RCP<const Number> b = bernoulli(s_);
                RCP<const Number> f = factorial(s_);
                zeta = divnum(pownum(integer(2), integer(s_ - 1)), f);
                zeta = mul(zeta, mul(pow(pi, s), abs(b)));
            } else {
                return make_rcp<const Zeta>(s, a);
            }
            if (a_ < 0)
                return add(zeta, harmonic(-a_, s_));
            return sub(zeta, harmonic(a_ - 1, s_));
        }
    }
    return make_rcp<const Zeta>(s, a);
}
 
RCP<const Basic> zeta(const RCP<const Basic> &s)
{
    return zeta(s, one);
}
 
Dirichlet_eta::Dirichlet_eta(const RCP<const Basic> &s) : OneArgFunction(s)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(s))
}
 
bool Dirichlet_eta::is_canonical(const RCP<const Basic> &s) const
{
    if (eq(*s, *one))
        return false;
    if (not(is_a<Zeta>(*zeta(s))))
        return false;
    return true;
}
 
RCP<const Basic> Dirichlet_eta::rewrite_as_zeta() const
{
    return mul(sub(one, pow(i2, sub(one, get_arg()))), zeta(get_arg()));
}
 
RCP<const Basic> Dirichlet_eta::create(const RCP<const Basic> &arg) const
{
    return dirichlet_eta(arg);
}
 
RCP<const Basic> dirichlet_eta(const RCP<const Basic> &s)
{
    if (is_a_Number(*s) and down_cast<const Number &>(*s).is_one()) {
        return log(i2);
    }
    RCP<const Basic> z = zeta(s);
    if (is_a<Zeta>(*z)) {
        return make_rcp<const Dirichlet_eta>(s);
    } else {
        return mul(sub(one, pow(i2, sub(one, s))), z);
    }
}
 
bool Erf::is_canonical(const RCP<const Basic> &arg) const
{
    if (is_a<Integer>(*arg) and down_cast<const Integer &>(*arg).is_zero())
        return false;
    if (could_extract_minus(*arg))
        return false;
    if (is_a_Number(*arg) and not down_cast<const Number &>(*arg).is_exact()) {
        return false;
    }
    return true;
}
 
RCP<const Basic> Erf::create(const RCP<const Basic> &arg) const
{
    return erf(arg);
}
 
RCP<const Basic> erf(const RCP<const Basic> &arg)
{
    if (is_a<Integer>(*arg) and down_cast<const Integer &>(*arg).is_zero()) {
        return zero;
    }
    if (is_a_Number(*arg)) {
        RCP<const Number> _arg = rcp_static_cast<const Number>(arg);
        if (not _arg->is_exact()) {
            return _arg->get_eval().erf(*_arg);
        }
    }
    RCP<const Basic> d;
    bool b = handle_minus(arg, outArg(d));
    if (b) {
        return neg(erf(d));
    }
    return make_rcp<const Erf>(d);
}
 
bool Erfc::is_canonical(const RCP<const Basic> &arg) const
{
    if (is_a<Integer>(*arg) and down_cast<const Integer &>(*arg).is_zero())
        return false;
    if (could_extract_minus(*arg))
        return false;
    if (is_a_Number(*arg) and not down_cast<const Number &>(*arg).is_exact()) {
        return false;
    }
    return true;
}
 
RCP<const Basic> Erfc::create(const RCP<const Basic> &arg) const
{
    return erfc(arg);
}
 
RCP<const Basic> erfc(const RCP<const Basic> &arg)
{
    if (is_a<Integer>(*arg) and down_cast<const Integer &>(*arg).is_zero()) {
        return one;
    }
    if (is_a_Number(*arg)) {
        RCP<const Number> _arg = rcp_static_cast<const Number>(arg);
        if (not _arg->is_exact()) {
            return _arg->get_eval().erfc(*_arg);
        }
    }
 
    RCP<const Basic> d;
    bool b = handle_minus(arg, outArg(d));
    if (b) {
        return add(integer(2), neg(erfc(d)));
    }
    return make_rcp<const Erfc>(d);
}
 
Gamma::Gamma(const RCP<const Basic> &arg) : OneArgFunction{arg}
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool Gamma::is_canonical(const RCP<const Basic> &arg) const
{
    if (is_a<Integer>(*arg))
        return false;
    if (is_a<Rational>(*arg)
        and (get_den(down_cast<const Rational &>(*arg).as_rational_class()))
                == 2) {
        return false;
    }
    if (is_a_Number(*arg) and not down_cast<const Number &>(*arg).is_exact()) {
        return false;
    }
    return true;
}
 
RCP<const Basic> Gamma::create(const RCP<const Basic> &arg) const
{
    return gamma(arg);
}
 
RCP<const Basic> gamma_positive_int(const RCP<const Basic> &arg)
{
    SYMENGINE_ASSERT(is_a<Integer>(*arg))
    RCP<const Integer> arg_ = rcp_static_cast<const Integer>(arg);
    SYMENGINE_ASSERT(arg_->is_positive())
    return factorial((arg_->subint(*one))->as_int());
}
 
RCP<const Basic> gamma_multiple_2(const RCP<const Basic> &arg)
{
    SYMENGINE_ASSERT(is_a<Rational>(*arg))
    RCP<const Rational> arg_ = rcp_static_cast<const Rational>(arg);
    SYMENGINE_ASSERT(get_den(arg_->as_rational_class()) == 2)
    RCP<const Integer> n, k;
    RCP<const Number> coeff;
    n = quotient_f(*(integer(mp_abs(get_num(arg_->as_rational_class())))),
                   *(integer(get_den(arg_->as_rational_class()))));
    if (arg_->is_positive()) {
        k = n;
        coeff = one;
    } else {
        n = n->addint(*one);
        k = n;
        if ((n->as_int() & 1) == 0) {
            coeff = one;
        } else {
            coeff = minus_one;
        }
    }
    int j = 1;
    for (int i = 3; i < 2 * k->as_int(); i = i + 2) {
        j = j * i;
    }
    coeff = mulnum(coeff, integer(j));
    if (arg_->is_positive()) {
        return div(mul(coeff, sqrt(pi)), pow(i2, n));
    } else {
        return div(mul(pow(i2, n), sqrt(pi)), coeff);
    }
}
 
RCP<const Basic> gamma(const RCP<const Basic> &arg)
{
    if (is_a<Integer>(*arg)) {
        RCP<const Integer> arg_ = rcp_static_cast<const Integer>(arg);
        if (arg_->is_positive()) {
            return gamma_positive_int(arg);
        } else {
            return ComplexInf;
        }
    } else if (is_a<Rational>(*arg)) {
        RCP<const Rational> arg_ = rcp_static_cast<const Rational>(arg);
        if ((get_den(arg_->as_rational_class())) == 2) {
            return gamma_multiple_2(arg);
        } else {
            return make_rcp<const Gamma>(arg);
        }
    } else if (is_a_Number(*arg)
               and not down_cast<const Number &>(*arg).is_exact()) {
        return down_cast<const Number &>(*arg).get_eval().gamma(*arg);
    }
    return make_rcp<const Gamma>(arg);
}
 
LowerGamma::LowerGamma(const RCP<const Basic> &s, const RCP<const Basic> &x)
    : TwoArgFunction(s, x)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(s, x))
}
 
bool LowerGamma::is_canonical(const RCP<const Basic> &s,
                              const RCP<const Basic> &x) const
{
    // Only special values are evaluated
    if (eq(*s, *one))
        return false;
    if (is_a<Integer>(*s)
        and down_cast<const Integer &>(*s).as_integer_class() > 1)
        return false;
    if (is_a<Integer>(*mul(i2, s)))
        return false;
#ifdef HAVE_SYMENGINE_MPFR
#if MPFR_VERSION_MAJOR > 3
    if (is_a<RealMPFR>(*s) && is_a<RealMPFR>(*x))
        return false;
#endif
#endif
    return true;
}
 
RCP<const Basic> LowerGamma::create(const RCP<const Basic> &a,
                                    const RCP<const Basic> &b) const
{
    return lowergamma(a, b);
}
 
RCP<const Basic> lowergamma(const RCP<const Basic> &s,
                            const RCP<const Basic> &x)
{
    // Only special values are being evaluated
    if (is_a<Integer>(*s)) {
        RCP<const Integer> s_int = rcp_static_cast<const Integer>(s);
        if (s_int->is_one()) {
            return sub(one, exp(mul(minus_one, x)));
        } else if (s_int->as_integer_class() > 1) {
            s_int = s_int->subint(*one);
            return sub(mul(s_int, lowergamma(s_int, x)),
                       mul(pow(x, s_int), exp(mul(minus_one, x))));
        } else {
            return make_rcp<const LowerGamma>(s, x);
        }
    } else if (is_a<Integer>(*(mul(i2, s)))) {
        RCP<const Number> s_num = rcp_static_cast<const Number>(s);
        s_num = subnum(s_num, one);
        if (eq(*s, *div(one, integer(2)))) {
            return mul(sqrt(pi),
                       erf(sqrt(x))); // base case for s of the form n/2
        } else if (s_num->is_positive()) {
            return sub(mul(s_num, lowergamma(s_num, x)),
                       mul(pow(x, s_num), exp(mul(minus_one, x))));
        } else {
            return div(add(lowergamma(add(s, one), x),
                           mul(pow(x, s), exp(mul(minus_one, x)))),
                       s);
        }
#ifdef HAVE_SYMENGINE_MPFR
#if MPFR_VERSION_MAJOR > 3
    } else if (is_a<RealMPFR>(*s) && is_a<RealMPFR>(*x)) {
        const auto &s_ = down_cast<const RealMPFR &>(*s).i.get_mpfr_t();
        const auto &x_ = down_cast<const RealMPFR &>(*x).i.get_mpfr_t();
        if (mpfr_cmp_si(x_, 0) >= 0) {
            mpfr_class t(std::max(mpfr_get_prec(s_), mpfr_get_prec(x_)));
            mpfr_class u(std::max(mpfr_get_prec(s_), mpfr_get_prec(x_)));
            mpfr_gamma_inc(t.get_mpfr_t(), s_, x_, MPFR_RNDN);
            mpfr_gamma(u.get_mpfr_t(), s_, MPFR_RNDN);
            mpfr_sub(t.get_mpfr_t(), u.get_mpfr_t(), t.get_mpfr_t(), MPFR_RNDN);
            return real_mpfr(std::move(t));
        } else {
            throw NotImplementedError("Not implemented.");
        }
#endif
#endif
    }
    return make_rcp<const LowerGamma>(s, x);
}
 
UpperGamma::UpperGamma(const RCP<const Basic> &s, const RCP<const Basic> &x)
    : TwoArgFunction(s, x)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(s, x))
}
 
bool UpperGamma::is_canonical(const RCP<const Basic> &s,
                              const RCP<const Basic> &x) const
{
    // Only special values are evaluated
    if (eq(*s, *one))
        return false;
    if (is_a<Integer>(*s)
        and down_cast<const Integer &>(*s).as_integer_class() > 1)
        return false;
    if (is_a<Integer>(*mul(i2, s)))
        return false;
#ifdef HAVE_SYMENGINE_MPFR
#if MPFR_VERSION_MAJOR > 3
    if (is_a<RealMPFR>(*s) && is_a<RealMPFR>(*x))
        return false;
#endif
#endif
    return true;
}
 
RCP<const Basic> UpperGamma::create(const RCP<const Basic> &a,
                                    const RCP<const Basic> &b) const
{
    return uppergamma(a, b);
}
 
RCP<const Basic> uppergamma(const RCP<const Basic> &s,
                            const RCP<const Basic> &x)
{
    // Only special values are being evaluated
    if (is_a<Integer>(*s)) {
        RCP<const Integer> s_int = rcp_static_cast<const Integer>(s);
        if (s_int->is_one()) {
            return exp(mul(minus_one, x));
        } else if (s_int->as_integer_class() > 1) {
            s_int = s_int->subint(*one);
            return add(mul(s_int, uppergamma(s_int, x)),
                       mul(pow(x, s_int), exp(mul(minus_one, x))));
        } else {
            // TODO: implement unpolarfy to handle this case
            return make_rcp<const LowerGamma>(s, x);
        }
    } else if (is_a<Integer>(*(mul(i2, s)))) {
        RCP<const Number> s_num = rcp_static_cast<const Number>(s);
        s_num = subnum(s_num, one);
        if (eq(*s, *div(one, integer(2)))) {
            return mul(sqrt(pi),
                       erfc(sqrt(x))); // base case for s of the form n/2
        } else if (s_num->is_positive()) {
            return add(mul(s_num, uppergamma(s_num, x)),
                       mul(pow(x, s_num), exp(mul(minus_one, x))));
        } else {
            return div(sub(uppergamma(add(s, one), x),
                           mul(pow(x, s), exp(mul(minus_one, x)))),
                       s);
        }
#ifdef HAVE_SYMENGINE_MPFR
#if MPFR_VERSION_MAJOR > 3
    } else if (is_a<RealMPFR>(*s) && is_a<RealMPFR>(*x)) {
        const auto &s_ = down_cast<const RealMPFR &>(*s).i.get_mpfr_t();
        const auto &x_ = down_cast<const RealMPFR &>(*x).i.get_mpfr_t();
        if (mpfr_cmp_si(x_, 0) >= 0) {
            mpfr_class t(std::max(mpfr_get_prec(s_), mpfr_get_prec(x_)));
            mpfr_gamma_inc(t.get_mpfr_t(), s_, x_, MPFR_RNDN);
            return real_mpfr(std::move(t));
        } else {
            throw NotImplementedError("Not implemented.");
        }
#endif
#endif
    }
    return make_rcp<const UpperGamma>(s, x);
}
 
bool LogGamma::is_canonical(const RCP<const Basic> &arg) const
{
    if (is_a<Integer>(*arg)) {
        RCP<const Integer> arg_int = rcp_static_cast<const Integer>(arg);
        if (not arg_int->is_positive()) {
            return false;
        }
        if (eq(*integer(1), *arg_int) or eq(*integer(2), *arg_int)
            or eq(*integer(3), *arg_int)) {
            return false;
        }
    }
    return true;
}
 
RCP<const Basic> LogGamma::rewrite_as_gamma() const
{
    return log(gamma(get_arg()));
}
 
RCP<const Basic> LogGamma::create(const RCP<const Basic> &arg) const
{
    return loggamma(arg);
}
 
RCP<const Basic> loggamma(const RCP<const Basic> &arg)
{
    if (is_a<Integer>(*arg)) {
        RCP<const Integer> arg_int = rcp_static_cast<const Integer>(arg);
        if (not arg_int->is_positive()) {
            return Inf;
        }
        if (eq(*integer(1), *arg_int) or eq(*integer(2), *arg_int)) {
            return zero;
        } else if (eq(*integer(3), *arg_int)) {
            return log(integer(2));
        }
    }
    return make_rcp<const LogGamma>(arg);
}
 
RCP<const Beta> Beta::from_two_basic(const RCP<const Basic> &x,
                                     const RCP<const Basic> &y)
{
    if (x->__cmp__(*y) == -1) {
        return make_rcp<const Beta>(y, x);
    }
    return make_rcp<const Beta>(x, y);
}
 
bool Beta::is_canonical(const RCP<const Basic> &x, const RCP<const Basic> &y)
{
    if (x->__cmp__(*y) == -1) {
        return false;
    }
    if (is_a<Integer>(*x)
        or (is_a<Rational>(*x)
            and (get_den(down_cast<const Rational &>(*x).as_rational_class()))
                    == 2)) {
        if (is_a<Integer>(*y)
            or (is_a<Rational>(*y)
                and (get_den(
                        down_cast<const Rational &>(*y).as_rational_class()))
                        == 2)) {
            return false;
        }
    }
    return true;
}
 
RCP<const Basic> Beta::rewrite_as_gamma() const
{
    return div(mul(gamma(get_arg1()), gamma(get_arg2())),
               gamma(add(get_arg1(), get_arg2())));
}
 
RCP<const Basic> Beta::create(const RCP<const Basic> &a,
                              const RCP<const Basic> &b) const
{
    return beta(a, b);
}
 
RCP<const Basic> beta(const RCP<const Basic> &x, const RCP<const Basic> &y)
{
    // Only special values are being evaluated
    if (eq(*add(x, y), *one)) {
        return ComplexInf;
    }
 
    if (is_a<Integer>(*x)) {
        RCP<const Integer> x_int = rcp_static_cast<const Integer>(x);
        if (x_int->is_positive()) {
            if (is_a<Integer>(*y)) {
                RCP<const Integer> y_int = rcp_static_cast<const Integer>(y);
                if (y_int->is_positive()) {
                    return div(
                        mul(gamma_positive_int(x), gamma_positive_int(y)),
                        gamma_positive_int(add(x, y)));
                } else {
                    return ComplexInf;
                }
            } else if (is_a<Rational>(*y)) {
                RCP<const Rational> y_ = rcp_static_cast<const Rational>(y);
                if (get_den(y_->as_rational_class()) == 2) {
                    return div(mul(gamma_positive_int(x), gamma_multiple_2(y)),
                               gamma_multiple_2(add(x, y)));
                } else {
                    return Beta::from_two_basic(x, y);
                }
            }
        } else {
            return ComplexInf;
        }
    }
 
    if (is_a<Integer>(*y)) {
        RCP<const Integer> y_int = rcp_static_cast<const Integer>(y);
        if (y_int->is_positive()) {
            if (is_a<Rational>(*x)) {
                RCP<const Rational> x_ = rcp_static_cast<const Rational>(x);
                if (get_den(x_->as_rational_class()) == 2) {
                    return div(mul(gamma_positive_int(y), gamma_multiple_2(x)),
                               gamma_multiple_2(add(x, y)));
                } else {
                    return Beta::from_two_basic(x, y);
                }
            }
        } else {
            return ComplexInf;
        }
    }
 
    if (is_a<const Rational>(*x)
        and get_den(down_cast<const Rational &>(*x).as_rational_class()) == 2) {
        if (is_a<Integer>(*y)) {
            RCP<const Integer> y_int = rcp_static_cast<const Integer>(y);
            if (y_int->is_positive()) {
                return div(mul(gamma_multiple_2(x), gamma_positive_int(y)),
                           gamma_multiple_2(add(x, y)));
            } else {
                return ComplexInf;
            }
        }
        if (is_a<const Rational>(*y)
            and get_den((down_cast<const Rational &>(*y)).as_rational_class())
                    == 2) {
            return div(mul(gamma_multiple_2(x), gamma_multiple_2(y)),
                       gamma_positive_int(add(x, y)));
        }
    }
    return Beta::from_two_basic(x, y);
}
 
bool PolyGamma::is_canonical(const RCP<const Basic> &n,
                             const RCP<const Basic> &x)
{
    if (is_a_Number(*x) and not(down_cast<const Number &>(*x)).is_positive()) {
        return false;
    }
    if (eq(*n, *zero)) {
        if (eq(*x, *one)) {
            return false;
        }
        if (is_a<Rational>(*x)) {
            auto x_ = rcp_static_cast<const Rational>(x);
            auto den = get_den(x_->as_rational_class());
            if (den == 2 or den == 3 or den == 4) {
                return false;
            }
        }
    }
    return true;
}
 
RCP<const Basic> PolyGamma::rewrite_as_zeta() const
{
    if (not is_a<Integer>(*get_arg1())) {
        return rcp_from_this();
    }
    RCP<const Integer> n = rcp_static_cast<const Integer>(get_arg1());
    if (not(n->is_positive())) {
        return rcp_from_this();
    }
    if ((n->as_int() & 1) == 0) {
        return neg(mul(factorial(n->as_int()), zeta(add(n, one), get_arg2())));
    } else {
        return mul(factorial(n->as_int()), zeta(add(n, one), get_arg2()));
    }
}
 
RCP<const Basic> PolyGamma::create(const RCP<const Basic> &a,
                                   const RCP<const Basic> &b) const
{
    return polygamma(a, b);
}
 
RCP<const Basic> polygamma(const RCP<const Basic> &n_,
                           const RCP<const Basic> &x_)
{
    // Only special values are being evaluated
    if (is_a_Number(*x_)
        and not(down_cast<const Number &>(*x_)).is_positive()) {
        return ComplexInf;
    }
    if (is_a<Integer>(*n_) and is_a<Integer>(*x_)) {
        auto n = down_cast<const Integer &>(*n_).as_int();
        auto x = down_cast<const Integer &>(*x_).as_int();
        if (n == 0) {
            return sub(harmonic(x - 1, 1), EulerGamma);
        } else if (n % 2 == 1) {
            return mul(factorial(n), zeta(add(n_, one), x_));
        }
    }
    if (eq(*n_, *zero)) {
        if (eq(*x_, *one)) {
            return neg(EulerGamma);
        }
        if (is_a<Rational>(*x_)) {
            RCP<const Rational> x = rcp_static_cast<const Rational>(x_);
            const auto den = get_den(x->as_rational_class());
            const auto num = get_num(x->as_rational_class());
            const integer_class r = num % den;
            RCP<const Basic> res;
            if (den == 2) {
                res = sub(mul(im2, log(i2)), EulerGamma);
            } else if (den == 3) {
                if (num == 1) {
                    res = add(neg(div(div(pi, i2), sqrt(i3))),
                              sub(div(mul(im3, log(i3)), i2), EulerGamma));
                } else {
                    res = add(div(div(pi, i2), sqrt(i3)),
                              sub(div(mul(im3, log(i3)), i2), EulerGamma));
                }
            } else if (den == 4) {
                if (num == 1) {
                    res = add(div(pi, im2), sub(mul(im3, log(i2)), EulerGamma));
                } else {
                    res = add(div(pi, i2), sub(mul(im3, log(i2)), EulerGamma));
                }
            } else {
                return make_rcp<const PolyGamma>(n_, x_);
            }
            rational_class a(0), f(r, den);
            for (unsigned long i = 0; i < (num - r) / den; ++i) {
                a += 1 / (f + i);
            }
            return add(Rational::from_mpq(a), res);
        }
    }
    return make_rcp<const PolyGamma>(n_, x_);
}
 
RCP<const Basic> digamma(const RCP<const Basic> &x)
{
    return polygamma(zero, x);
}
 
RCP<const Basic> trigamma(const RCP<const Basic> &x)
{
    return polygamma(one, x);
}
 
Abs::Abs(const RCP<const Basic> &arg) : OneArgFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool Abs::is_canonical(const RCP<const Basic> &arg) const
{
    if (is_a<Integer>(*arg) or is_a<Rational>(*arg) or is_a<Complex>(*arg))
        return false;
    if (is_a_Number(*arg) and not down_cast<const Number &>(*arg).is_exact()) {
        return false;
    }
    if (is_a<Abs>(*arg)) {
        return false;
    }
 
    if (could_extract_minus(*arg)) {
        return false;
    }
 
    return true;
}
 
RCP<const Basic> Abs::create(const RCP<const Basic> &arg) const
{
    return abs(arg);
}
 
RCP<const Basic> abs(const RCP<const Basic> &arg)
{
    if (is_a<Integer>(*arg)) {
        RCP<const Integer> arg_ = rcp_static_cast<const Integer>(arg);
        if (arg_->is_negative()) {
            return arg_->neg();
        } else {
            return arg_;
        }
    } else if (is_a<Rational>(*arg)) {
        RCP<const Rational> arg_ = rcp_static_cast<const Rational>(arg);
        if (arg_->is_negative()) {
            return arg_->neg();
        } else {
            return arg_;
        }
    } else if (is_a<Complex>(*arg)) {
        RCP<const Complex> arg_ = rcp_static_cast<const Complex>(arg);
        return sqrt(Rational::from_mpq(arg_->real_ * arg_->real_
                                       + arg_->imaginary_ * arg_->imaginary_));
    } else if (is_a_Number(*arg)
               and not down_cast<const Number &>(*arg).is_exact()) {
        return down_cast<const Number &>(*arg).get_eval().abs(*arg);
    }
    if (is_a<Abs>(*arg)) {
        return arg;
    }
 
    RCP<const Basic> d;
    handle_minus(arg, outArg(d));
    return make_rcp<const Abs>(d);
}
 
Max::Max(const vec_basic &&arg) : MultiArgFunction(std::move(arg))
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(get_vec()))
}
 
bool Max::is_canonical(const vec_basic &arg) const
{
    if (arg.size() < 2)
        return false;
 
    bool non_number_exists = false;
 
    for (const auto &p : arg) {
        if (is_a<Complex>(*p) or is_a<Max>(*p))
            return false;
        if (not is_a_Number(*p))
            non_number_exists = true;
    }
    if (not std::is_sorted(arg.begin(), arg.end(), RCPBasicKeyLess()))
        return false;
 
    return non_number_exists; // all arguments cant be numbers
}
 
RCP<const Basic> Max::create(const vec_basic &a) const
{
    return max(a);
}
 
RCP<const Basic> max(const vec_basic &arg)
{
    bool number_set = false;
    RCP<const Number> max_number, difference;
    set_basic new_args;
 
    for (const auto &p : arg) {
        if (is_a<Complex>(*p))
            throw SymEngineException("Complex can't be passed to max!");
 
        if (is_a_Number(*p)) {
            if (not number_set) {
                max_number = rcp_static_cast<const Number>(p);
 
            } else {
                if (eq(*p, *Inf)) {
                    return Inf;
                } else if (eq(*p, *NegInf)) {
                    continue;
                }
                difference = down_cast<const Number &>(*p).sub(*max_number);
 
                if (difference->is_zero() and not difference->is_exact()) {
                    if (max_number->is_exact())
                        max_number = rcp_static_cast<const Number>(p);
                } else if (difference->is_positive()) {
                    max_number = rcp_static_cast<const Number>(p);
                }
            }
            number_set = true;
 
        } else if (is_a<Max>(*p)) {
            for (const auto &l : down_cast<const Max &>(*p).get_args()) {
                if (is_a_Number(*l)) {
                    if (not number_set) {
                        max_number = rcp_static_cast<const Number>(l);
 
                    } else {
                        difference = rcp_static_cast<const Number>(l)->sub(
                            *max_number);
 
                        if (difference->is_zero()
                            and not difference->is_exact()) {
                            if (max_number->is_exact())
                                max_number = rcp_static_cast<const Number>(l);
                        } else if (difference->is_positive()) {
                            max_number = rcp_static_cast<const Number>(l);
                        }
                    }
                    number_set = true;
                } else {
                    new_args.insert(l);
                }
            }
        } else {
            new_args.insert(p);
        }
    }
 
    if (number_set)
        new_args.insert(max_number);
 
    vec_basic final_args(new_args.size());
    std::copy(new_args.begin(), new_args.end(), final_args.begin());
 
    if (final_args.size() > 1) {
        return make_rcp<const Max>(std::move(final_args));
    } else if (final_args.size() == 1) {
        return final_args[0];
    } else {
        throw SymEngineException("Empty vec_basic passed to max!");
    }
}
 
Min::Min(const vec_basic &&arg) : MultiArgFunction(std::move(arg))
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(get_vec()))
}
 
bool Min::is_canonical(const vec_basic &arg) const
{
    if (arg.size() < 2)
        return false;
 
    bool non_number_exists = false;
 
    for (const auto &p : arg) {
        if (is_a<Complex>(*p) or is_a<Min>(*p))
            return false;
        if (not is_a_Number(*p))
            non_number_exists = true;
    }
    if (not std::is_sorted(arg.begin(), arg.end(), RCPBasicKeyLess()))
        return false;
 
    return non_number_exists; // all arguments cant be numbers
}
 
RCP<const Basic> Min::create(const vec_basic &a) const
{
    return min(a);
}
 
RCP<const Basic> min(const vec_basic &arg)
{
    bool number_set = false;
    RCP<const Number> min_number, difference;
    set_basic new_args;
 
    for (const auto &p : arg) {
        if (is_a<Complex>(*p))
            throw SymEngineException("Complex can't be passed to min!");
 
        if (is_a_Number(*p)) {
            if (not number_set) {
                min_number = rcp_static_cast<const Number>(p);
 
            } else {
                if (eq(*p, *Inf)) {
                    continue;
                } else if (eq(*p, *NegInf)) {
                    return NegInf;
                }
                difference = min_number->sub(*rcp_static_cast<const Number>(p));
 
                if (difference->is_zero() and not difference->is_exact()) {
                    if (min_number->is_exact())
                        min_number = rcp_static_cast<const Number>(p);
                } else if (difference->is_positive()) {
                    min_number = rcp_static_cast<const Number>(p);
                }
            }
            number_set = true;
 
        } else if (is_a<Min>(*p)) {
            for (const auto &l : down_cast<const Min &>(*p).get_args()) {
                if (is_a_Number(*l)) {
                    if (not number_set) {
                        min_number = rcp_static_cast<const Number>(l);
 
                    } else {
                        difference = min_number->sub(
                            *rcp_static_cast<const Number>(l));
 
                        if (difference->is_zero()
                            and not difference->is_exact()) {
                            if (min_number->is_exact())
                                min_number = rcp_static_cast<const Number>(l);
                        } else if (difference->is_positive()) {
                            min_number = rcp_static_cast<const Number>(l);
                        }
                    }
                    number_set = true;
                } else {
                    new_args.insert(l);
                }
            }
        } else {
            new_args.insert(p);
        }
    }
 
    if (number_set)
        new_args.insert(min_number);
 
    vec_basic final_args(new_args.size());
    std::copy(new_args.begin(), new_args.end(), final_args.begin());
 
    if (final_args.size() > 1) {
        return make_rcp<const Min>(std::move(final_args));
    } else if (final_args.size() == 1) {
        return final_args[0];
    } else {
        throw SymEngineException("Empty vec_basic passed to min!");
    }
}
 
UnevaluatedExpr::UnevaluatedExpr(const RCP<const Basic> &arg)
    : OneArgFunction(arg)
{
    SYMENGINE_ASSIGN_TYPEID()
    SYMENGINE_ASSERT(is_canonical(arg))
}
 
bool UnevaluatedExpr::is_canonical(const RCP<const Basic> &arg) const
{
    return true;
}
 
RCP<const Basic> UnevaluatedExpr::create(const RCP<const Basic> &arg) const
{
    return make_rcp<const UnevaluatedExpr>(arg);
}
 
RCP<const Basic> unevaluated_expr(const RCP<const Basic> &arg)
{
    return make_rcp<const UnevaluatedExpr>(arg);
}
 
} // namespace SymEngine