#ifndef SCALAR_HPP
#define SCALAR_HPP

#include "symbol.hpp"
#define MPFR_REAL_DATA_PUBLIC
#include "mpfr/real.hpp"
#include <sstream>
#include <iostream>
#include <iomanip>
#include <complex>
#include <cstdlib>

template<typename>
class scalar_t;
template<typename>
class complex_t;

typedef double          real_type;
typedef mpfr::real<128> large_real;

typedef scalar_t<real_type>   scalar;
typedef scalar_t<large_real>  scalar_hp;
typedef complex_t<real_type>  complex;
typedef complex_t<large_real> complex_hp;

template<typename T>
struct numeric_traits {
    static std::string to_string(const T& v, const interpreter_context& context) {
        std::stringstream ss;
        ss << std::setprecision(context.digits) << v;
        return ss.str();
    }
};

template<>
struct numeric_traits<large_real> {
    static std::string to_string(const large_real& v, const interpreter_context& context);
};

template<typename ntype>
ntype factorial(ntype d) {
    if (d > 1.0)
        return d*factorial(d - 1.0);
    else
        return 1.0;
}

template<typename T>
class scalar_t : public symbol_impl {
public :
    typedef T ntype;

    ntype value;

    explicit scalar_t(ntype d) : value(d) {}

    scalar_t* clone() const override {
        return new scalar_t(value);
    }

    scalar_t* evaluate(const interpreter_context&) const override {
        return clone();
    }

    std::string serialize(const interpreter_context& context) const override {
        return numeric_traits<ntype>::to_string(value, context);
    }

    const std::string& type() const override {
        return type_name;
    }

    bool is_litteral() const override {
        return true;
    }

    static const std::string type_name;

    static void reg(interpreter_context& context);
};

template<typename ntype>
const std::string scalar_t<ntype>::type_name = "scalar";

template<typename T>
class complex_t : public symbol_impl {
public :
    typedef T ntype;
    typedef std::complex<ntype> ctype;
    typedef ::scalar_t<ntype> scalar_t;

    ctype  value;
    ntype& rvalue;
    ntype& ivalue;

    explicit complex_t(ctype c) : value(c),
        rvalue(reinterpret_cast<ntype*>(&value)[0]),
        ivalue(reinterpret_cast<ntype*>(&value)[1]) {}

    explicit complex_t(ntype r, ntype i) : value(r, i),
        rvalue(reinterpret_cast<ntype*>(&value)[0]),
        ivalue(reinterpret_cast<ntype*>(&value)[1]) {}

    complex_t* clone() const override {
        return new complex_t(value);
    }

    symbol_impl* evaluate(const interpreter_context&) const override {
        return clone();
    }

    std::string serialize(const interpreter_context& context) const override {
        std::string s;
        if (ivalue == 0.0) {
            s = numeric_traits<ntype>::to_string(rvalue, context);
        } else if (ivalue == 1.0 || ivalue == -1.0) {
            if (rvalue != 0.0) {
                if (ivalue > 0.0)
                    s = numeric_traits<ntype>::to_string(rvalue, context) + "+i";
                else
                    s = numeric_traits<ntype>::to_string(rvalue, context) + "-i";
            } else {
                if (ivalue > 0.0)
                    s = "i";
                else
                    s = "-i";
            }
        } else {
            if (rvalue != 0.0) {
                if (ivalue > 0.0) {
                    s = numeric_traits<ntype>::to_string(rvalue, context) + "+"
                        + numeric_traits<ntype>::to_string(ivalue, context) + "i";
                } else {
                    s = numeric_traits<ntype>::to_string(rvalue, context) +
                        numeric_traits<ntype>::to_string(ivalue, context) + "i";
                }
            } else {
                s = numeric_traits<ntype>::to_string(ivalue, context) + "i";
            }
        }
        return s;
    }

    const std::string& type() const override {
        return type_name;
    }

    bool is_litteral() const override {
        return true;
    }

    static const std::string type_name;

    static void reg(interpreter_context& context);
};

template<typename ntype>
const std::string complex_t<ntype>::type_name = "complex";

template<typename ntype>
void scalar_t<ntype>::reg(interpreter_context& context) {
        context.register_operator<scalar_t>("+", [](const symbol& s) {
            return s.as<scalar_t>().clone();
        });
        context.register_operator<scalar_t>("-", [](const symbol& s) {
            return new scalar_t(-s.as<scalar_t>().value);
        });
        context.register_operator<scalar_t, scalar_t>("+", [](const symbol& s1, const symbol& s2) {
            return new scalar_t(s1.as<scalar_t>().value + s2.as<scalar_t>().value);
        });
        context.register_operator<scalar_t, scalar_t>("-", [](const symbol& s1, const symbol& s2) {
            return new scalar_t(s1.as<scalar_t>().value - s2.as<scalar_t>().value);
        });
        context.register_operator<scalar_t, scalar_t>("*", [](const symbol& s1, const symbol& s2) {
            return new scalar_t(s1.as<scalar_t>().value * s2.as<scalar_t>().value);
        });
        context.register_operator<scalar_t, scalar_t>("/", [](const symbol& s1, const symbol& s2) {
            return new scalar_t(s1.as<scalar_t>().value / s2.as<scalar_t>().value);
        });
        context.register_operator<scalar_t, scalar_t>("^", [](const symbol& s1, const symbol& s2) {
            return new scalar_t(pow(s1.as<scalar_t>().value, s2.as<scalar_t>().value));
        });
        context.register_operator<scalar_t, scalar_t>("%", [](const symbol& s1, const symbol& s2) {
            return new scalar_t(fmod(s1.as<scalar_t>().value, s2.as<scalar_t>().value));
        });

        context.register_operator<scalar_t>("!", [](const symbol& s) {
            return new scalar_t(factorial(s.as<scalar_t>().value));
        });

        context.register_function("rand", [](const std::vector<symbol>& args) {
            return new scalar_t(2.0*rand()/RAND_MAX - 1.0);
        });
        context.register_function<scalar_t, scalar_t>("rand", [](const std::vector<symbol>& args) {
            typename scalar_t::ntype rmin = args[0].as<scalar_t>().value, rmax = args[1].as<scalar_t>().value;
            return new scalar_t((rmax - rmin)*rand()/RAND_MAX + rmin);
        });
        context.register_function<scalar_t, scalar_t>("min", [](const std::vector<symbol>& args) {
            return new scalar_t(std::min(args[0].as<scalar_t>().value, args[1].as<scalar_t>().value));
        });
        context.register_function<scalar_t, scalar_t>("max", [](const std::vector<symbol>& args) {
            return new scalar_t(std::max(args[0].as<scalar_t>().value, args[1].as<scalar_t>().value));
        });
        context.register_function<scalar_t>("round", [](const std::vector<symbol>& args) {
            return new scalar_t(round(args[0].as<scalar_t>().value));
        });
        context.register_function<scalar_t>({"round_up", "ceil"}, [](const std::vector<symbol>& args) {
            return new scalar_t(ceil(args[0].as<scalar_t>().value));
        });
        context.register_function<scalar_t>({"round_down", "floor"}, [](const std::vector<symbol>& args) {
            return new scalar_t(floor(args[0].as<scalar_t>().value));
        });
        context.register_function<scalar_t>({"square_root", "sqrt"}, [](const std::vector<symbol>& args) {
            return new scalar_t(sqrt(args[0].as<scalar_t>().value));
        });
        context.register_function<scalar_t>("exp", [](const std::vector<symbol>& args) {
            return new scalar_t(exp(args[0].as<scalar_t>().value));
        });
        context.register_function<scalar_t>("log", [](const std::vector<symbol>& args) {
            return new scalar_t(log(args[0].as<scalar_t>().value));
        });
        context.register_function<scalar_t>("cos", [](const std::vector<symbol>& args) {
            return new scalar_t(cos(args[0].as<scalar_t>().value));
        });
        context.register_function<scalar_t>({"cosh", "ch"}, [](const std::vector<symbol>& args) {
            return new scalar_t(cosh(args[0].as<scalar_t>().value));
        });
        context.register_function<scalar_t>("acos", [](const std::vector<symbol>& args) {
            return new scalar_t(acos(args[0].as<scalar_t>().value));
        });
        context.register_function<scalar_t>("sin", [](const std::vector<symbol>& args) {
            return new scalar_t(sin(args[0].as<scalar_t>().value));
        });
        context.register_function<scalar_t>({"sinh", "sh"}, [](const std::vector<symbol>& args) {
            return new scalar_t(sinh(args[0].as<scalar_t>().value));
        });
        context.register_function<scalar_t>("asin", [](const std::vector<symbol>& args) {
            return new scalar_t(asin(args[0].as<scalar_t>().value));
        });
        context.register_function<scalar_t>("tan", [](const std::vector<symbol>& args) {
            return new scalar_t(tan(args[0].as<scalar_t>().value));
        });
        context.register_function<scalar_t>({"tanh", "th"}, [](const std::vector<symbol>& args) {
            return new scalar_t(tanh(args[0].as<scalar_t>().value));
        });
        context.register_function<scalar_t>("atan", [](const std::vector<symbol>& args) {
            return new scalar_t(atan(args[0].as<scalar_t>().value));
        });
        context.register_function<scalar_t,scalar_t>("atan", [](const std::vector<symbol>& args) {
            return new scalar_t(atan2(args[0].as<scalar_t>().value, args[1].as<scalar_t>().value));
        });
        context.register_function<scalar_t>("abs", [](const std::vector<symbol>& args) {
            return new scalar_t(fabs(args[0].as<scalar_t>().value));
        });

        context.register_function("pi", [](const std::vector<symbol>& args) {
            return new scalar_t(acos(-1.0));
        });
    }

template<typename ntype>
void complex_t<ntype>::reg(interpreter_context& context) {
    context.register_operator<complex_t>("+", [](const symbol& s) {
        return s.as<complex_t>().clone();
    });
    context.register_operator<complex_t>("-", [](const symbol& s) {
        const complex_t& c = s.as<complex_t>();
        return new complex_t(-c.rvalue, -c.ivalue);
    });

    context.register_operator<complex_t, complex_t>("+", [](const symbol& s1, const symbol& s2) {
        const complex_t& c1 = s1.as<complex_t>();
        const complex_t& c2 = s2.as<complex_t>();
        return new complex_t(c1.rvalue + c2.rvalue, c1.ivalue + c2.ivalue);
    });
    context.register_operator<complex_t, complex_t>("-", [](const symbol& s1, const symbol& s2) {
        const complex_t& c1 = s1.as<complex_t>();
        const complex_t& c2 = s2.as<complex_t>();
        return new complex_t(c1.rvalue - c2.rvalue, c1.ivalue - c2.ivalue);
    });
    context.register_operator<complex_t, complex_t>("*", [](const symbol& s1, const symbol& s2) {
        const complex_t& c1 = s1.as<complex_t>();
        const complex_t& c2 = s2.as<complex_t>();
        return new complex_t(c1.rvalue*c2.rvalue - c1.ivalue*c2.ivalue, c1.ivalue*c2.rvalue + c1.rvalue*c2.ivalue);
    });
    context.register_operator<complex_t, complex_t>("/", [](const symbol& s1, const symbol& s2) {
        const complex_t& c1 = s1.as<complex_t>();
        const complex_t& c2 = s2.as<complex_t>();
        typename scalar_t::ntype denom = c2.rvalue*c2.rvalue + c2.ivalue*c2.ivalue;
        return new complex_t((c1.rvalue*c2.rvalue + c1.ivalue*c2.ivalue)/denom, (c1.ivalue*c2.rvalue - c1.rvalue*c2.ivalue)/denom);
    });

    context.register_commutative_operator<complex_t, scalar_t>("+", [](const symbol& s1, const symbol& s2) {
        const scalar_t& s = s2.as<scalar_t>();
        const complex_t& c = s1.as<complex_t>();
        return new complex_t(c.rvalue + s.value, c.ivalue);
    });
    context.register_operator<complex_t, scalar_t>("-", [](const symbol& s1, const symbol& s2) {
        const scalar_t& s = s2.as<scalar_t>();
        const complex_t& c = s1.as<complex_t>();
        return new complex_t(c.rvalue - s.value, c.ivalue);
    });
    context.register_operator<scalar_t, complex_t>("-", [](const symbol& s1, const symbol& s2) {
        const scalar_t& s = s1.as<scalar_t>();
        const complex_t& c = s2.as<complex_t>();
        return new complex_t(s.value - c.rvalue, -c.ivalue);
    });
    context.register_commutative_operator<complex_t, scalar_t>("*", [](const symbol& s1, const symbol& s2) {
        const scalar_t& s = s2.as<scalar_t>();
        const complex_t& c = s1.as<complex_t>();
        return new complex_t(c.rvalue*s.value, c.ivalue*s.value);
    });
    context.register_operator<complex_t, scalar_t>("/", [](const symbol& s1, const symbol& s2) {
        const scalar_t& s = s2.as<scalar_t>();
        const complex_t& c = s1.as<complex_t>();
        return new complex_t(c.rvalue/s.value, c.ivalue/s.value);
    });
    context.register_operator<scalar_t, complex_t>("/", [](const symbol& s1, const symbol& s2) {
        const scalar_t& s = s1.as<scalar_t>();
        const complex_t& c = s2.as<complex_t>();
        typename scalar_t::ntype denom = c.rvalue*c.rvalue + c.ivalue*c.ivalue;
        return new complex_t((c.rvalue*s.value)/denom, -c.ivalue/denom);
    });

    context.register_function<complex_t, complex_t>("min", [](const std::vector<symbol>& args) {
        const complex_t& c1 = args[0].as<complex_t>();
        const complex_t& c2 = args[1].as<complex_t>();
        return new complex_t(
            std::min(c1.rvalue, c2.rvalue),
            std::min(c1.ivalue, c2.ivalue)
        );
    });
    context.register_function<complex_t, scalar_t>("min", [](const std::vector<symbol>& args) {
        const complex_t& c1 = args[0].as<complex_t>();
        const scalar_t&  s2 = args[1].as<scalar_t>();
        return new complex_t(
            std::min(c1.rvalue, s2.value),
            std::min(c1.ivalue, ntype(0.0))
        );
    });
    context.register_function<scalar_t, complex_t>("min", [](const std::vector<symbol>& args) {
        const complex_t& c1 = args[1].as<complex_t>();
        const scalar_t&  s2 = args[0].as<scalar_t>();
        return new complex_t(
            std::min(c1.rvalue, s2.value),
            std::min(c1.ivalue, ntype(0.0))
        );
    });
    context.register_function<complex_t, complex_t>("max", [](const std::vector<symbol>& args) {
        const complex_t& c1 = args[0].as<complex_t>();
        const complex_t& c2 = args[1].as<complex_t>();
        return new complex_t(
            std::max(c1.rvalue, c2.rvalue),
            std::max(c1.ivalue, c2.ivalue)
        );
    });
    context.register_function<complex_t, scalar_t>("max", [](const std::vector<symbol>& args) {
        const complex_t& c1 = args[0].as<complex_t>();
        const scalar_t&  s2 = args[1].as<scalar_t>();
        return new complex_t(
            std::max(c1.rvalue, s2.value),
            std::max(c1.ivalue, ntype(0.0))
        );
    });
    context.register_function<scalar_t, complex_t>("max", [](const std::vector<symbol>& args) {
        const complex_t& c1 = args[1].as<complex_t>();
        const scalar_t&  s2 = args[0].as<scalar_t>();
        return new complex_t(
            std::max(c1.rvalue, s2.value),
            std::max(c1.ivalue, ntype(0.0))
        );
    });
    context.register_function<complex_t>("round", [](const std::vector<symbol>& args) {
        const complex_t& c = args[0].as<complex_t>();
        return new complex_t(round(c.rvalue), round(c.ivalue));
    });
    context.register_function<complex_t>({"round_up", "ceil"}, [](const std::vector<symbol>& args) {
        const complex_t& c = args[0].as<complex_t>();
        return new complex_t(ceil(c.rvalue), ceil(c.ivalue));
    });
    context.register_function<complex_t>({"round_down", "floor"}, [](const std::vector<symbol>& args) {
        const complex_t& c = args[0].as<complex_t>();
        return new complex_t(floor(c.rvalue), floor(c.ivalue));
    });
    context.register_function<complex_t>({"square_root", "sqrt"}, [](const std::vector<symbol>& args) {
        return new complex_t(sqrt(args[0].as<complex_t>().value));
    });
    context.register_function<complex_t>("exp", [](const std::vector<symbol>& args) {
        return new complex_t(exp(args[0].as<complex_t>().value));
    });
    context.register_function<complex_t>("log", [](const std::vector<symbol>& args) {
        return new complex_t(log(args[0].as<complex_t>().value));
    });
    context.register_function<complex_t>("cos", [](const std::vector<symbol>& args) {
        return new complex_t(cos(args[0].as<complex_t>().value));
    });
    context.register_function<complex_t>({"cosh", "ch"}, [](const std::vector<symbol>& args) {
        return new complex_t(cosh(args[0].as<complex_t>().value));
    });
    context.register_function<complex_t>("acos", [](const std::vector<symbol>& args) {
        return new complex_t(acos(args[0].as<complex_t>().value));
    });
    context.register_function<complex_t>("sin", [](const std::vector<symbol>& args) {
        return new complex_t(sin(args[0].as<complex_t>().value));
    });
    context.register_function<complex_t>({"sinh", "sh"}, [](const std::vector<symbol>& args) {
        return new complex_t(sinh(args[0].as<complex_t>().value));
    });
    context.register_function<complex_t>("asin", [](const std::vector<symbol>& args) {
        return new complex_t(asin(args[0].as<complex_t>().value));
    });
    context.register_function<complex_t>("tan", [](const std::vector<symbol>& args) {
        return new complex_t(tan(args[0].as<complex_t>().value));
    });
    context.register_function<complex_t>({"tanh", "th"}, [](const std::vector<symbol>& args) {
        return new complex_t(tanh(args[0].as<complex_t>().value));
    });
    context.register_function<complex_t>("atan", [](const std::vector<symbol>& args) {
        return new complex_t(atan(args[0].as<complex_t>().value));
    });
    context.register_function<complex_t>({"abs", "norm"}, [](const std::vector<symbol>& args) {
        return new complex_t(norm(args[0].as<complex_t>().value));
    });

    context.register_function("i", [](const std::vector<symbol>& args) {
        return new complex_t(ntype(0.0), ntype(1.0));
    });
}
#endif
