//
// Created by trotfunky on 09/05/19.
//

#ifndef SNIPPETS_POLYNOMIAL_TPP
#define SNIPPETS_POLYNOMIAL_TPP

#include <vector>
#include <map>
#include <iostream>
#include <math.h>

template <typename T> class Polynomial;

/**
 * Class allowing creation, manipulation and calculation of polynomials
 * Factors are stored from the lowest power to the highest
 * The polynomial has at least one factor, which defaults to zero
 * @tparam T Has to be an arithmetic type (std::is_arithmetic<T>::value is true)
 */
template <typename T>
class Polynomial {
public:
    Polynomial();
    explicit Polynomial(const std::vector<T>&);
    explicit Polynomial(const std::map<int,T>&);

    template <typename T1>
    explicit Polynomial(const Polynomial<T1>&);

    int getDegree() const;
    Polynomial<T> getNthDerivative(int n) const;

    template <typename T1>
    bool equals(const Polynomial<T1>&) const;

    template <typename T1>
    auto add(const Polynomial<T1>& operand) const -> Polynomial<decltype(T{} + T1{})>;

    template <typename T1>
    friend std::ostream& operator<<(std::ostream&, const Polynomial<T1>&);

    template<typename T1>
    T1 operator()(const T1&) const;

    T operator[](int) const;

private:
    std::vector<T> factors;
    double factorial(unsigned int n) const;

    template <typename T2>
    friend class Polynomial;
};

////////////////////
///              ///
/// Definitions  ///
///              ///
////////////////////

template<typename T>
Polynomial<T>::Polynomial()
{
    // Not bad, but could be better with C++20. Not that useful either
    static_assert(std::is_arithmetic<T>::value,"Polynomial must be of an arithmetic type!");
    factors = {0};
}

template<typename T>
Polynomial<T>::Polynomial(const std::vector<T>& polynomialFactors) : Polynomial<T>()
{
    if(polynomialFactors.size()>0)
    {
        factors = polynomialFactors;
    }
}

template<typename T>
Polynomial<T>::Polynomial(const std::map<int, T>& polynomialFactors) : Polynomial<T>()
{
    int degree = polynomialFactors.rbegin()->first;
    if(degree >= 0)
    {
        factors.pop_back();
    }
    factors.reserve(degree);

    for(int i = 0;i<=degree;i++)
    {
        auto nextFactor = polynomialFactors.find(i);
        if(nextFactor != polynomialFactors.end())
        {
            factors.push_back(nextFactor->second);
        }
        else
        {
            factors.push_back(0);
        }
    }
}

template <typename T>
template <typename T1>
Polynomial<T>::Polynomial(const Polynomial<T1>& copied)
{
    static_assert(std::is_arithmetic<T>::value,"Polynomial must be of an arithmetic type!");

    for(const T1& factor : copied.factors)
    {
        factors.push_back(static_cast<T>(factor));
    }
}

/**
 * Returns the degree of the polynomial.
 * @tparam T Type of the polynomial factors
 * @return Degree of the polynomial, -1 if it is the null-polynomial
 */
template<typename T>
int Polynomial<T>::getDegree() const
{
    for(int i = factors.size()-1;i>=0;i--)
    {
        if(factors[i] != 0)
        {
            return(i);
        }
    }
    return(-1);
}

/**
 * Computes the nth derivative of the polynomial
 * @tparam T Polyniaml factor type
 * @param  n Order of the derivative
 * @return Derived polynomial, null-polynomial if the order is greater than the degree
 */
template<typename T>
Polynomial<T> Polynomial<T>::getNthDerivative(int n) const
{
    if(n<=0)
    {
        return(Polynomial<T>(*this));
    }

    if(n>getDegree())
    {
        return(Polynomial<T>());
    }

    std::vector<T> newFactors;

    for(int i = n;i<factors.size();i++)
    {
        newFactors.push_back(factors[i]*factorial(i)/factorial(i-n));
    }

    return(Polynomial<T>(newFactors));
}

template<typename T>
double Polynomial<T>::factorial(unsigned int n) const
{
    double value = 1;
    for(unsigned int i = 2;i<=n;i++)
    {
        value *= i;
    }

    return(value);
}

template<typename T>
template<typename T1>
bool Polynomial<T>::equals(const Polynomial<T1>& operand) const
{
    if(getDegree() != operand.getDegree())
    {
        return(false);
    }
    if(getDegree() == -1)
    {
        return(true);
    }
    const std::vector<T1>& p2Factors = operand.factors;

    for(int i = factors.size();i>=0;i--)
    {
        if(factors[i] != p2Factors[i])
        {
            return(false);
        }
    }

    return(true);
}

template<typename T>
template<typename T1>
auto Polynomial<T>::add(const Polynomial<T1>& operand) const -> Polynomial<decltype(T{} + T1{})>
{
    bool isLargest = true;

    int largestSize = factors.size();
    int smallestSize = operand.factors.size();

    if(getDegree() < operand.getDegree())
    {
        isLargest = false;

        smallestSize = factors.size();
        largestSize = operand.factors.size();
    }

    std::vector<decltype(T{} + T1{})> resultPolynomial = {};
    resultPolynomial.reserve(largestSize);

    for(int i = 0;i<smallestSize;i++)
    {
        resultPolynomial.push_back(factors[i]+operand[i]);
    }

    for(int i = smallestSize;i<largestSize;i++)
    {
        resultPolynomial.push_back((isLargest ? factors[i] : operand[i]));
    }

    return Polynomial<decltype(T{} + T1{})>(resultPolynomial);
}

////////////////////
///              ///
///   Operators  ///
///              ///
////////////////////

template<typename T1,typename T2>
bool operator==(const Polynomial<T1>& p1, const Polynomial<T2>& p2)
{
    return(p1.equals(p2));
}


/**
 * Computes the sum of the two polynomials by adding factor by factor at first and pushing back
 * the leftover factors from the largest polynomial if their are different in degree.
 * @tparam T  Type of the first polynomial
 * @tparam R  Type of the second polynomial
 * @param p1  First polynomial
 * @param p2  Second polynomial
 * @return    Sum of the two polynomials, null polynomial if both of them are null. The type is deduced from the addition
 *              of the first factors of the two polynomials.
 */
template<typename T1, typename T2>
auto operator+(const Polynomial<T1>& p1, const Polynomial<T2>& p2) -> Polynomial<decltype(p1[0] + p2[0])>
{
    return(p1.add(p2));
}

template <typename T>
std::ostream& operator<<(std::ostream& ostream, const Polynomial<T>& polynomial)
{
    if(polynomial.getDegree() < 0)
    {
        ostream << "0 ";
        return ostream;
    }

    const std::vector<T>& factors = polynomial.factors;

    if(factors[0] > 0)
    {
        ostream << factors[0] << " ";
    }
    else if(factors[0] < 0)
    {
        ostream << std::showpos << factors[0] << " ";
    }

    for(int i = 1;i<factors.size();i++)
    {
        if(factors[i] != 0)
        {
            ostream << std::showpos << factors[i] << "*X^" << std::noshowpos << i << " ";
        }
    }

    return(ostream);
}

template<typename T>
template<typename T1>
T1 Polynomial<T>::operator()(const T1& input) const
{
    static_assert(std::is_arithmetic<T1>::value,"Cannot evaluate polynomial at a non arithmetic value !");

    if(getDegree() == -1)
    {
        return(0);
    }

    T1 returnValue = 0;
    for(int i = 0;i<factors.size();i++)
    {
        returnValue += factors[i] * pow(input,i);
    }

    return(returnValue);
}

template<typename T>
T Polynomial<T>::operator[](const int index) const
{
    return(factors.at(index));
}


#endif //SNIPPETS_POLYNOMIAL_TPP