#include <poincare/addition.h>
  #include <poincare/multiplication.h>
  #include <poincare/subtraction.h>
  #include <poincare/power.h>
  #include <poincare/opposite.h>
  #include <poincare/undefined.h>
  #include <poincare/matrix.h>
  extern "C" {
  #include <assert.h>
  #include <stdlib.h>
  }
  
  namespace Poincare {
  
  Expression::Type Addition::type() const {
    return Type::Addition;
  }
  
  Expression * Addition::clone() const {
    if (numberOfOperands() == 0) {
      return new Addition();
    }
    return new Addition(operands(), numberOfOperands(), true);
  }
  
  int Addition::polynomialDegree(char symbolName) const {
    int degree = 0;
    for (int i = 0; i < numberOfOperands(); i++) {
      int d = operand(i)->polynomialDegree(symbolName);
      if (d < 0) {
        return -1;
      }
      degree = d > degree ? d : degree;
    }
    return degree;
  }
  
  int Addition::privateGetPolynomialCoefficients(char symbolName, Expression * coefficients[]) const {
    int deg = polynomialDegree(symbolName);
    if (deg < 0 || deg > k_maxPolynomialDegree) {
      return -1;
    }
    for (int k = 0; k < deg+1; k++) {
      coefficients[k] = new Addition();
    }
    Expression * intermediateCoefficients[k_maxNumberOfPolynomialCoefficients];
    for (int i = 0; i < numberOfOperands(); i++) {
      int d = operand(i)->privateGetPolynomialCoefficients(symbolName, intermediateCoefficients);
      assert(d < k_maxNumberOfPolynomialCoefficients);
      for (int j = 0; j < d+1; j++) {
        static_cast<Addition *>(coefficients[j])->addOperand(intermediateCoefficients[j]);
      }
    }
    return deg;
  }
  
  /* Layout */
  
  bool Addition::needParenthesisWithParent(const Expression * e) const {
    Type types[] = {Type::Subtraction, Type::Opposite, Type::Multiplication, Type::Division, Type::Power, Type::Factorial};
    return e->isOfType(types, 6);
  }
  
  /* Simplication */
  
  Expression * Addition::shallowReduce(Context& context, AngleUnit angleUnit) {
    Expression * e = Expression::shallowReduce(context, angleUnit);
    if (e != this) {
      return e;
    }
    /* Step 1: Addition is associative, so let's start by merging children which
     * also are additions themselves. */
    int i = 0;
    int initialNumberOfOperands = numberOfOperands();
    while (i < initialNumberOfOperands) {
      Expression * o = editableOperand(i);
      if (o->type() == Type::Addition) {
        mergeOperands(static_cast<Addition *>(o));
        continue;
      }
      i++;
    }
  
    // Step 2: Sort the operands
    sortOperands(Expression::SimplificationOrder, true);
  
  #if MATRIX_EXACT_REDUCING
    /* Step 2bis: get rid of matrix */
    int n = 1;
    int m = 1;
    /* All operands have been simplified so if any operand contains a matrix, it
     * is at the root node of the operand. Moreover, thanks to the simplification
     * order, all matrix operands (if any) are the last operands. */
    Expression * lastOperand = editableOperand(numberOfOperands()-1);
    if (lastOperand->type() == Type::Matrix) {
      // Create in-place the matrix of addition M (in place of the last operand)
      Matrix * resultMatrix = static_cast<Matrix *>(lastOperand);
      n = resultMatrix->numberOfRows();
      m = resultMatrix->numberOfColumns();
      removeOperand(resultMatrix, false);
      /* Scan (starting at the end) accross the addition operands to find any
       * other matrix */
      int i = numberOfOperands()-1;
      while (i >= 0 && operand(i)->type() == Type::Matrix) {
        Matrix * currentMatrix = static_cast<Matrix *>(editableOperand(i));
        int on = currentMatrix->numberOfRows();
        int om = currentMatrix->numberOfColumns();
        if (on != n || om != m) {
          return replaceWith(new Undefined(), true);
        }
        // Dispatch the current matrix operands in the created additions matrix
        for (int j = 0; j < n*m; j++) {
          Addition * a = new Addition();
          Expression * resultMatrixEntryJ = resultMatrix->editableOperand(j);
          resultMatrix->replaceOperand(resultMatrixEntryJ, a, false);
          a->addOperand(currentMatrix->editableOperand(j));
          a->addOperand(resultMatrixEntryJ);
          a->shallowReduce(context, angleUnit);
        }
        currentMatrix->detachOperands();
        removeOperand(currentMatrix, true);
        i--;
      }
      // Distribute the remaining addition on matrix operands
      for (int i = 0; i < n*m; i++) {
        Addition * a = static_cast<Addition *>(clone());
        Expression * entryI = resultMatrix->editableOperand(i);
        resultMatrix->replaceOperand(entryI, a, false);
        a->addOperand(entryI);
        a->shallowReduce(context, angleUnit);
      }
      return replaceWith(resultMatrix, true)->shallowReduce(context, angleUnit);
    }
  #endif
  
    /* Step 3: Factorize like terms. Thanks to the simplification order, those are
     * next to each other at this point. */
    i = 0;
    while (i < numberOfOperands()-1) {
      Expression * o1 = editableOperand(i);
      Expression * o2 = editableOperand(i+1);
      if (o1->type() == Type::Rational && o2->type() == Type::Rational) {
        Rational r1 = *static_cast<Rational *>(o1);
        Rational r2 = *static_cast<Rational *>(o2);
        Rational a = Rational::Addition(r1, r2);
        replaceOperand(o1, new Rational(a), true);
        removeOperand(o2, true);
        continue;
      }
      if (TermsHaveIdenticalNonRationalFactors(o1, o2)) {
        factorizeOperands(o1, o2, context, angleUnit);
        continue;
      }
      i++;
    }
  
    /* Step 4: Let's remove zeroes if there's any. It's important to do this after
     * having factorized because factorization can lead to new zeroes. For example
     * pi+(-1)*pi. We don't remove the last zero if it's the only operand left
     * though. */
    i = 0;
    while (i < numberOfOperands()) {
      Expression * o = editableOperand(i);
      if (o->type() == Type::Rational && static_cast<Rational *>(o)->isZero() && numberOfOperands() > 1) {
        removeOperand(o, true);
        continue;
      }
      i++;
    }
  
    // Step 5: Let's remove the addition altogether if it has a single operand
    Expression * result = squashUnaryHierarchy();
  
    // Step 6: Last but not least, let's put everything under a common denominator
    if (result == this && parent()->type() != Type::Addition) {
      // squashUnaryHierarchy didn't do anything: we're not an unary hierarchy
      result = factorizeOnCommonDenominator(context, angleUnit);
    }
  
    return result;
  }
  
  Expression * Addition::factorizeOnCommonDenominator(Context & context, AngleUnit angleUnit) {
    // We want to turn (a/b+c/d+e/b) into (a*d+b*c+e*d)/(b*d)
  
    // Step 1: We want to compute the common denominator, b*d
    Multiplication * commonDenominator = new Multiplication();
    for (int i = 0; i < numberOfOperands(); i++) {
      Expression * denominator = operand(i)->cloneDenominator(context, angleUnit);
      if (denominator) {
        // Make commonDenominator = LeastCommonMultiple(commonDenominator, denominator);
        commonDenominator->addMissingFactors(denominator, context, angleUnit);
        delete denominator;
      }
    }
    if (commonDenominator->numberOfOperands() == 0) {
      delete commonDenominator;
      // If commonDenominator is empty this means that no operand was a fraction.
      return this;
    }
  
    // Step 2: Create the numerator. We start with this being a/b+c/d+e/b and we
    // want to create numerator = a/b*b*d + c/d*b*d + e/b*b*d
    Addition * numerator = new Addition();
    for (int i=0; i < numberOfOperands(); i++) {
      Multiplication * m = new Multiplication(operand(i), commonDenominator, true);
      numerator->addOperand(m);
      m->privateShallowReduce(context, angleUnit, true, false);
    }
    // Step 3: Add the denominator
    Power * inverseDenominator = new Power(commonDenominator, new Rational(-1), false);
    Multiplication * result = new Multiplication(numerator, inverseDenominator, false);
  
    // Step 4: Simplify the numerator to a*d + c*b + e*d
    numerator->shallowReduce(context, angleUnit);
  
    // Step 5: Simplify the denominator (in case it's a rational number)
    commonDenominator->deepReduce(context, angleUnit);
    inverseDenominator->shallowReduce(context, angleUnit);
  
    /* Step 6: We simplify the resulting multiplication forbidding any
     * distribution of multiplication on additions (to avoid an infinite loop). */
    return static_cast<Multiplication *>(replaceWith(result, true))->privateShallowReduce(context, angleUnit, false, true);
  }
  
  void Addition::factorizeOperands(Expression * e1, Expression * e2, Context & context, AngleUnit angleUnit) {
    /* This function factorizes two operands which only differ by a rational
     * factor. For example, if this is Addition(2*pi, 3*pi), then 2*pi and 3*pi
     * could be merged, and this turned into Addition(5*pi). */
    assert(e1->parent() == this && e2->parent() == this);
  
    // Step 1: Find the new rational factor
    Rational * r = new Rational(Rational::Addition(RationalFactor(e1), RationalFactor(e2)));
  
    // Step 2: Get rid of one of the operands
    removeOperand(e2, true);
  
    // Step 3: Use the new rational factor. Create a multiplication if needed
    Multiplication * m = nullptr;
    if (e1->type() == Type::Multiplication) {
      m = static_cast<Multiplication *>(e1);
    } else {
      m = new Multiplication();
      e1->replaceWith(m, false);
      m->addOperand(e1);
    }
    if (m->operand(0)->type() == Type::Rational) {
      m->replaceOperand(m->operand(0), r, true);
    } else {
      m->addOperand(r);
    }
  
    // Step 4: Reduce the multiplication (in case the new rational factor is zero)
    m->shallowReduce(context, angleUnit);
  }
  
  const Rational Addition::RationalFactor(Expression * e) {
    if (e->type() == Type::Multiplication && e->operand(0)->type() == Type::Rational) {
      return *(static_cast<const Rational *>(e->operand(0)));
    }
    return Rational(1);
  }
  
  static inline int NumberOfNonRationalFactors(const Expression * e) {
    if (e->type() != Expression::Type::Multiplication) {
      return 1; // Or (e->type() != Type::Rational);
    }
    int result = e->numberOfOperands();
    if (e->operand(0)->type() == Expression::Type::Rational) {
      result--;
    }
    return result;
  }
  
  static inline const Expression * FirstNonRationalFactor(const Expression * e) {
    if (e->type() != Expression::Type::Multiplication) {
      return e;
    }
    if (e->operand(0)->type() == Expression::Type::Rational) {
      return e->numberOfOperands() > 1 ? e->operand(1) : nullptr;
    }
    return e->operand(0);
  }
  
  bool Addition::TermsHaveIdenticalNonRationalFactors(const Expression * e1, const Expression * e2) {
    /* Given two expressions, say wether they only differ by a rational factor.
     * For example, 2*pi and pi do, 2*pi and 2*ln(2) don't. */
  
    int numberOfNonRationalFactorsInE1 = NumberOfNonRationalFactors(e1);
    int numberOfNonRationalFactorsInE2 = NumberOfNonRationalFactors(e2);
  
    if (numberOfNonRationalFactorsInE1 != numberOfNonRationalFactorsInE2) {
      return false;
    }
  
    int numberOfNonRationalFactors = numberOfNonRationalFactorsInE1;
    if (numberOfNonRationalFactors == 1) {
      return FirstNonRationalFactor(e1)->isIdenticalTo(FirstNonRationalFactor(e2));
    } else {
      assert(numberOfNonRationalFactors > 1);
      return Multiplication::HaveSameNonRationalFactors(e1, e2);
    }
  }
  
  Expression * Addition::shallowBeautify(Context & context, AngleUnit angleUnit) {
    /* Beautifying Addition essentially consists in adding Subtractions if needed.
     * In practice, we want to turn "a+(-1)*b" into "a-b". Or, more precisely, any
     * "a+(-r)*b" into "a-r*b" where r is a positive Rational.
     * Note: the process will slightly differ if the negative product occurs on
     * the first term: we want to turn "Addition(Multiplication(-1,b))" into
     * "Opposite(b)".
     * Last but not least, special care must be taken when iterating over operands
     * since we may remove some during the process. */
  
    for (int i=0; i<numberOfOperands(); i++) {
      if (operand(i)->type() != Type::Multiplication || operand(i)->operand(0)->type() != Type::Rational || operand(i)->operand(0)->sign() != Sign::Negative) {
        // Ignore terms which are not like "(-r)*a"
        continue;
      }
  
      Multiplication * m = static_cast<Multiplication *>(editableOperand(i));
  
      if (static_cast<const Rational *>(m->operand(0))->isMinusOne()) {
        m->removeOperand(m->operand(0), true);
      } else {
        m->editableOperand(0)->setSign(Sign::Positive, context, angleUnit);
      }
      Expression * subtractant = m->squashUnaryHierarchy();
  
      if (i == 0) {
        Opposite * o = new Opposite(subtractant, true);
        replaceOperand(subtractant, o, true);
      } else {
        const Expression * op1 = operand(i-1);
        removeOperand(op1, false);
        Subtraction * s = new Subtraction(op1, subtractant, false);
        replaceOperand(subtractant, s, false);
        /* CAUTION: we removed an operand. So we need to decrement i to make sure
         * the next iteration is actually on the next operand. */
        i--;
      }
    }
  
    return squashUnaryHierarchy();
  }
  
  /* Evaluation */
  
  template<typename T>
  std::complex<T> Addition::compute(const std::complex<T> c, const std::complex<T> d) {
    return c+d;
  }
  
  template std::complex<float> Poincare::Addition::compute<float>(std::complex<float>, std::complex<float>);
  template std::complex<double> Poincare::Addition::compute<double>(std::complex<double>, std::complex<double>);
  
  template MatrixComplex<float> Addition::computeOnMatrices<float>(const MatrixComplex<float>,const MatrixComplex<float>);
  template MatrixComplex<double> Addition::computeOnMatrices<double>(const MatrixComplex<double>,const MatrixComplex<double>);
  
  template MatrixComplex<float> Addition::computeOnComplexAndMatrix<float>(std::complex<float> const, const MatrixComplex<float>);
  template MatrixComplex<double> Addition::computeOnComplexAndMatrix<double>(std::complex<double> const, const MatrixComplex<double>);
  
  }