#include <poincare/multiplication.h>
  #include <poincare/addition.h>
  #include <poincare/arithmetic.h>
  #include <poincare/division.h>
  #include <poincare/matrix.h>
  #include <poincare/opposite.h>
  #include <poincare/parenthesis.h>
  #include <poincare/power.h>
  #include <poincare/rational.h>
  #include <poincare/simplification_root.h>
  #include <poincare/subtraction.h>
  #include <poincare/tangent.h>
  #include <poincare/undefined.h>
  #include <cmath>
  #include <ion.h>
  extern "C" {
  #include <assert.h>
  #include <stdlib.h>
  }
  
  namespace Poincare {
  
  Expression::Type Multiplication::type() const {
    return Expression::Type::Multiplication;
  }
  
  Expression * Multiplication::clone() const {
    if (numberOfOperands() == 0) {
      return new Multiplication();
    }
    return new Multiplication(operands(), numberOfOperands(), true);
  }
  
  int Multiplication::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;
    }
    return degree;
  }
  
  int Multiplication::privateGetPolynomialCoefficients(char symbolName, Expression * coefficients[]) const {
    int deg = polynomialDegree(symbolName);
    if (deg < 0 || deg > k_maxPolynomialDegree) {
      return -1;
    }
    // Initialization of coefficients
    for (int i = 1; i <= deg; i++) {
      coefficients[i] = new Rational(0);
    }
    coefficients[0] = new Rational(1);
  
    Expression * intermediateCoefficients[k_maxNumberOfPolynomialCoefficients];
    // Let's note result = a(0)+a(1)*X+a(2)*X^2+a(3)*x^3+..
    for (int i = 0; i < numberOfOperands(); i++) {
      // operand(i) = b(0)+b(1)*X+b(2)*X^2+b(3)*x^3+...
      int degI = operand(i)->privateGetPolynomialCoefficients(symbolName, intermediateCoefficients);
      assert(degI <= k_maxPolynomialDegree);
      for (int j = deg; j > 0; j--) {
        // new coefficients[j] = b(0)*a(j)+b(1)*a(j-1)+b(2)*a(j-2)+...
        Addition * a = new Addition();
        int jbis = j > degI ? degI : j;
        for (int l = 0; l <= jbis ; l++) {
          // Always copy the a and b coefficients are they are used multiple times
          a->addOperand(new Multiplication(intermediateCoefficients[l], coefficients[j-l], true));
        }
        /* a(j) and b(j) are used only to compute coefficient at rank >= j, we
         * can delete them as we compute new coefficient by decreasing ranks. */
        delete coefficients[j];
        if (j <= degI) { delete intermediateCoefficients[j]; };
        coefficients[j] = a;
      }
      // new coefficients[0] = a(0)*b(0)
      coefficients[0] = new Multiplication(coefficients[0], intermediateCoefficients[0], false);
    }
    return deg;
  }
  
  bool Multiplication::needParenthesisWithParent(const Expression * e) const {
    Type types[] = {Type::Division, Type::Power, Type::Factorial};
    return e->isOfType(types, 3);
  }
  
  ExpressionLayout * Multiplication::createLayout(PrintFloat::Mode floatDisplayMode, int numberOfSignificantDigits) const {
    const char middleDotString[] = {Ion::Charset::MiddleDot, 0};
    return LayoutEngine::createInfixLayout(this, floatDisplayMode, numberOfSignificantDigits, middleDotString);
  }
  
  int Multiplication::writeTextInBuffer(char * buffer, int bufferSize, PrintFloat::Mode floatDisplayMode, int numberOfSignificantDigits) const {
    const char multiplicationString[] = {Ion::Charset::MultiplicationSign, 0};
    return LayoutEngine::writeInfixExpressionTextInBuffer(this, buffer, bufferSize, floatDisplayMode, numberOfSignificantDigits, multiplicationString);
  }
  
  Expression::Sign Multiplication::sign() const {
    int sign = 1;
    for (int i = 0; i < numberOfOperands(); i++) {
      sign *= (int)operand(i)->sign();
    }
    return (Sign)sign;
  }
  
  Expression * Multiplication::setSign(Sign s, Context & context, AngleUnit angleUnit) {
    assert(s == Sign::Positive);
    for (int i = 0; i < numberOfOperands(); i++) {
      if (operand(i)->sign() == Sign::Negative) {
        editableOperand(i)->setSign(s, context, angleUnit);
      }
    }
    return shallowReduce(context, angleUnit);
  }
  
  template<typename T>
  std::complex<T> Multiplication::compute(const std::complex<T> c, const std::complex<T> d) {
    return c*d;
  }
  
  template<typename T>
  MatrixComplex<T> Multiplication::computeOnMatrices(const MatrixComplex<T> m, const MatrixComplex<T> n) {
    if (m.numberOfColumns() != n.numberOfRows()) {
      return MatrixComplex<T>::Undefined();
    }
    std::complex<T> * operands = new std::complex<T> [m.numberOfRows()*n.numberOfColumns()];
    for (int i = 0; i < m.numberOfRows(); i++) {
      for (int j = 0; j < n.numberOfColumns(); j++) {
        std::complex<T> c(0.0);
        for (int k = 0; k < m.numberOfColumns(); k++) {
          c += m.complexOperand(i*m.numberOfColumns()+k)*n.complexOperand(k*n.numberOfColumns()+j);
        }
        operands[i*n.numberOfColumns()+j] = c;
      }
    }
    MatrixComplex<T> result = MatrixComplex<T>(operands, m.numberOfRows(), n.numberOfColumns());
    delete[] operands;
    return result;
  }
  
  template<typename T>
  void Multiplication::computeOnArrays(T * m, T * n, T * result, int mNumberOfColumns, int mNumberOfRows, int nNumberOfColumns) {
    for (int i = 0; i < mNumberOfRows; i++) {
      for (int j = 0; j < nNumberOfColumns; j++) {
        T res = 0.0f;
        for (int k = 0; k < mNumberOfColumns; k++) {
          res+= m[i*mNumberOfColumns+k]*n[k*nNumberOfColumns+j];
        }
        result[i*nNumberOfColumns+j] = res;
      }
    }
  }
  
  bool Multiplication::HaveSameNonRationalFactors(const Expression * e1, const Expression * e2) {
    int numberOfNonRationalFactors1 = e1->operand(0)->type() == Type::Rational ? e1->numberOfOperands()-1 : e1->numberOfOperands();
    int numberOfNonRationalFactors2 = e2->operand(0)->type() == Type::Rational ? e2->numberOfOperands()-1 : e2->numberOfOperands();
    if (numberOfNonRationalFactors1 != numberOfNonRationalFactors2) {
      return false;
    }
    int firstNonRationalOperand1 = e1->operand(0)->type() == Type::Rational ? 1 : 0;
    int firstNonRationalOperand2 = e2->operand(0)->type() == Type::Rational ? 1 : 0;
    for (int i = 0; i < numberOfNonRationalFactors1; i++) {
      if (!(e1->operand(firstNonRationalOperand1+i)->isIdenticalTo(e2->operand(firstNonRationalOperand2+i)))) {
        return false;
      }
    }
    return true;
  }
  
  static inline const Expression * Base(const Expression * e) {
    if (e->type() == Expression::Type::Power) {
      return e->operand(0);
    }
    return e;
  }
  
  Expression * Multiplication::shallowReduce(Context& context, AngleUnit angleUnit) {
    return privateShallowReduce(context, angleUnit, true, true);
  }
  
  Expression * Multiplication::privateShallowReduce(Context & context, AngleUnit angleUnit, bool shouldExpand, bool canBeInterrupted) {
    Expression * e = Expression::shallowReduce(context, angleUnit);
    if (e != this) {
      return e;
    }
    /* Step 1: Multiplication is associative, so let's start by merging children
     * which also are multiplications themselves. */
    mergeMultiplicationOperands();
  
    /* Step 2: If any of the operand is zero, the multiplication result is zero */
    for (int i = 0; i < numberOfOperands(); i++) {
      const Expression * o = operand(i);
      if (o->type() == Type::Rational && static_cast<const Rational *>(o)->isZero()) {
        return replaceWith(new Rational(0), true);
      }
    }
  
    // Step 3: Sort the operands
    sortOperands(SimplificationOrder, canBeInterrupted);
  
  #if MATRIX_EXACT_REDUCING
    /* Step 3bis: 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) {
      Matrix * resultMatrix = static_cast<Matrix *>(lastOperand);
      // Use the last matrix operand as the final matrix
      n = resultMatrix->numberOfRows();
      m = resultMatrix->numberOfColumns();
      /* Scan accross the multiplication operands to find any other matrix:
       * (the last operand is the result matrix so we start at
       * numberOfOperands()-2)*/
      int k = numberOfOperands()-2;
      while (k >= 0 && operand(k)->type() == Type::Matrix) {
        Matrix * currentMatrix = static_cast<Matrix *>(editableOperand(k));
        int on = currentMatrix->numberOfRows();
        int om = currentMatrix->numberOfColumns();
        if (om != n) {
          return replaceWith(new Undefined(), true);
        }
        // Create the matrix resulting of the multiplication of the current matrix and the result matrix
       /*                        resultMatrix
        *                          i2= 0..m
        *                         +-+-+-+-+-+
        *                         | | | | | |
        *                         +-+-+-+-+-+
        *                  j=0..n | | | | | |
        *                         +-+-+-+-+-+
        *                         | | | | | |
        *                         +-+-+-+-+-+
        *        currentMatrix
        *           j=0..om
        *         +---+---+---+   +-+-+-+-+-+
        *         |   |   |   |   | | | | | |
        *         +---+---+---+   +-+-+-+-+-+
        *i1=0..on |   |   |   |   | |e| | | |
        *         +---+---+---+   +-+-+-+-+-+
        *         |   |   |   |   | | | | | |
        *         +---+---+---+   +-+-+-+-+-+
        * */
        Expression ** newMatrixOperands = new Expression * [on*m];
        for (int e = 0; e < on*m; e++) {
          newMatrixOperands[e] = new Addition();
          int i2 = e%m;
          int i1 = e/m;
          for (int j = 0; j < n; j++) {
            Expression * mult = new Multiplication(currentMatrix->editableOperand(j+om*i1), resultMatrix->editableOperand(j*m+i2), true);
            static_cast<Addition *>(newMatrixOperands[e])->addOperand(mult);
            mult->shallowReduce(context, angleUnit);
          }
          Reduce(&newMatrixOperands[e], context, angleUnit, false);
        }
        n = on;
        removeOperand(currentMatrix, true);
        resultMatrix = static_cast<Matrix *>(resultMatrix->replaceWith(new Matrix(newMatrixOperands, n, m, false), true));
        k--;
      }
      removeOperand(resultMatrix, false);
      // Distribute the remaining multiplication on matrix operands
      for (int i = 0; i < n*m; i++) {
        Multiplication * m = static_cast<Multiplication *>(clone());
        Expression * entryI = resultMatrix->editableOperand(i);
        resultMatrix->replaceOperand(entryI, m, false);
        m->addOperand(entryI);
        m->shallowReduce(context, angleUnit);
      }
      return replaceWith(resultMatrix, true)->shallowReduce(context, angleUnit);
    }
  #endif
  
    /* Step 4: Gather like terms. For example, turn pi^2*pi^3 into pi^5. Thanks to
     * the simplification order, such terms are guaranteed to be next to each
     * other. */
    int i = 0;
    while (i < numberOfOperands()-1) {
      Expression * oi = editableOperand(i);
      Expression * oi1 = editableOperand(i+1);
      if (TermsHaveIdenticalBase(oi, oi1)) {
        bool shouldFactorizeBase = true;
        if (TermHasRationalBase(oi)) {
          /* Combining powers of a given rational isn't straightforward. Indeed,
           * there are two cases we want to deal with:
           *  - 2*2^(1/2) or 2*2^pi, we want to keep as-is
           *  - 2^(1/2)*2^(3/2) we want to combine. */
          shouldFactorizeBase = oi->type() == Type::Power && oi1->type() == Type::Power;
        }
        if (shouldFactorizeBase) {
          factorizeBase(oi, oi1, context, angleUnit);
          continue;
        }
      } else if (TermHasRationalBase(oi) && TermHasRationalBase(oi1) && TermsHaveIdenticalExponent(oi, oi1)) {
        factorizeExponent(oi, oi1, context, angleUnit);
        continue;
      }
      i++;
    }
  
    /* Step 5: We look for terms of form sin(x)^p*cos(x)^q with p, q rational of
     *opposite signs. We replace them by either:
     * - tan(x)^p*cos(x)^(p+q) if |p|<|q|
     * - tan(x)^(-q)*sin(x)^(p+q) otherwise */
    for (int i = 0; i < numberOfOperands(); i++) {
      Expression * o1 = editableOperand(i);
      if (Base(o1)->type() == Type::Sine && TermHasRationalExponent(o1)) {
        const Expression * x = Base(o1)->operand(0);
        /* Thanks to the SimplificationOrder, Cosine-base factors are after
         * Sine-base factors */
        for (int j = i+1; j < numberOfOperands(); j++) {
          Expression * o2 = editableOperand(j);
          if (Base(o2)->type() == Type::Cosine && TermHasRationalExponent(o2) && Base(o2)->operand(0)->isIdenticalTo(x)) {
            factorizeSineAndCosine(o1, o2, context, angleUnit);
            break;
          }
        }
      }
    }
    /* Replacing sin/cos by tan factors may have mixed factors and factors are
     * guaranteed to be sorted (according ot SimplificationOrder) at the end of
     * shallowReduce */
    sortOperands(SimplificationOrder, true);
  
    /* Step 6: We remove rational operands that appeared in the middle of sorted
     * operands. It's important to do this after having factorized because
     * factorization can lead to new ones. Indeed:
     * pi^(-1)*pi-> 1
     * i*i -> -1
     * 2^(1/2)*2^(1/2) -> 2
     * sin(x)*cos(x) -> 1*tan(x)
     * Last, we remove the only rational operand if it is one and not the only
     * operand. */
    i = 1;
    while (i < numberOfOperands()) {
      Expression * o = editableOperand(i);
      if (o->type() == Type::Rational && static_cast<Rational *>(o)->isOne()) {
        removeOperand(o, true);
        continue;
      }
      if (o->type() == Type::Rational) {
        if (operand(0)->type() == Type::Rational) {
          Rational * o0 = static_cast<Rational *>(editableOperand(0));
          Rational m = Rational::Multiplication(*o0, *(static_cast<Rational *>(o)));
          replaceOperand(o0, new Rational(m), true);
          removeOperand(o, true);
        } else {
          removeOperand(o, false);
          addOperandAtIndex(o, 0);
        }
        continue;
      }
      i++;
    }
    if (operand(0)->type() == Type::Rational && static_cast<Rational *>(editableOperand(0))->isOne() && numberOfOperands() > 1) {
      removeOperand(editableOperand(0), true);
    }
  
  
    /* Step 7: Expand multiplication over addition operands if any. For example,
     * turn (a+b)*c into a*c + b*c. We do not want to do this step right now if
     * the parent is a multiplication to avoid missing factorization such as
     * (x+y)^(-1)*((a+b)*(x+y)).
     * Note: This step must be done after Step 4, otherwise we wouldn't be able to
     * reduce expressions such as (x+y)^(-1)*(x+y)(a+b). */
    if (shouldExpand && parent()->type() != Type::Multiplication) {
      for (int i=0; i<numberOfOperands(); i++) {
        if (operand(i)->type() == Type::Addition) {
          return distributeOnOperandAtIndex(i, context, angleUnit);
        }
      }
    }
  
    // Step 8: Let's remove the multiplication altogether if it has one operand
    Expression * result = squashUnaryHierarchy();
  
    return result;
  }
  
  void Multiplication::mergeMultiplicationOperands() {
    // Multiplication is associative: a*(b*c)->a*b*c
    int i = 0;
    int initialNumberOfOperands = numberOfOperands();
    while (i < initialNumberOfOperands) {
      Expression * o = editableOperand(i);
      if (o->type() == Type::Multiplication) {
        mergeOperands(static_cast<Multiplication *>(o)); // TODO: ensure that matrix operands are not swapped to implement MATRIX_EXACT_REDUCING
        continue;
      }
      i++;
    }
  }
  
  void Multiplication::factorizeSineAndCosine(Expression * o1, Expression * o2, Context & context, AngleUnit angleUnit) {
    assert(o1->parent() == this && o2->parent() == this);
    /* This function turn sin(x)^p * cos(x)^q into either:
     * - tan(x)^p*cos(x)^(p+q) if |p|<|q|
     * - tan(x)^(-q)*sin(x)^(p+q) otherwise */
    const Expression * x = Base(o1)->operand(0);
    Rational p = o1->type() == Type::Power ? *(static_cast<Rational *>(o1->editableOperand(1))) : Rational(1);
    Rational q = o2->type() == Type::Power ? *(static_cast<Rational *>(o2->editableOperand(1))) : Rational(1);
    /* If p and q have the same sign, we cannot replace them by a tangent */
    if ((int)p.sign()*(int)q.sign() > 0) {
      return;
    }
    Rational sumPQ = Rational::Addition(p, q);
    Rational absP = p;
    absP.setSign(Sign::Positive);
    Rational absQ = q;
    absQ.setSign(Sign::Positive);
    Expression * tan = new Tangent(x, true);
    if (Rational::NaturalOrder(absP, absQ) < 0) {
      if (o1->type() == Type::Power) {
        o1->replaceOperand(o1->operand(0), tan, true);
      } else {
        replaceOperand(o1, tan, true);
        o1 = tan;
      }
      o1->shallowReduce(context, angleUnit);
      if (o2->type() == Type::Power) {
        o2->replaceOperand(o2->operand(1), new Rational(sumPQ), true);
      } else {
        Expression * newO2 = new Power(o2, new Rational(sumPQ), false);
        replaceOperand(o2, newO2, false);
        o2 = newO2;
      }
      o2->shallowReduce(context, angleUnit);
    } else {
      if (o2->type() == Type::Power) {
        o2->replaceOperand(o2->operand(1), new Rational(Rational::Multiplication(q, Rational(-1))), true);
        o2->replaceOperand(o2->operand(0), tan, true);
      } else {
        Expression * newO2 = new Power(tan, new Rational(-1), false);
        replaceOperand(o2, newO2, true);
        o2 = newO2;
      }
      o2->shallowReduce(context, angleUnit);
      if (o1->type() == Type::Power) {
        o1->replaceOperand(o1->operand(1), new Rational(sumPQ), true);
      } else {
        Expression * newO1 = new Power(o1, new Rational(sumPQ), false);
        replaceOperand(o1, newO1, false);
        o1 = newO1;
      }
      o1->shallowReduce(context, angleUnit);
    }
  }
  
  void Multiplication::factorizeBase(Expression * e1, Expression * e2, Context & context, AngleUnit angleUnit) {
    /* This function factorizes two operands which have a common base. For example
     * if this is Multiplication(pi^2, pi^3), then pi^2 and pi^3 could be merged
     * and this turned into Multiplication(pi^5). */
    assert(TermsHaveIdenticalBase(e1, e2));
  
    // Step 1: Find the new exponent
    Expression * s = new Addition(CreateExponent(e1), CreateExponent(e2), false);
  
    // Step 2: Get rid of one of the operands
    removeOperand(e2, true);
  
    // Step 3: Use the new exponent
    Power * p = nullptr;
    if (e1->type() == Type::Power) {
      // If e1 is a power, replace the initial exponent with the new one
      e1->replaceOperand(e1->operand(1), s, true);
      p = static_cast<Power *>(e1);
    } else {
      // Otherwise, create a new Power node
      p = new Power(e1, s, false);
      replaceOperand(e1, p, false);
    }
  
    // Step 4: Reduce the new power
    s->shallowReduce(context, angleUnit); // pi^2*pi^3 -> pi^(2+3) -> pi^5
    Expression * reducedP = p->shallowReduce(context, angleUnit); // pi^2*pi^-2 -> pi^0 -> 1
    /* Step 5: Reducing the new power might have turned it into a multiplication,
     * ie: 12^(1/2) -> 2*3^(1/2). In that case, we need to merge the multiplication
     * node with this. */
    if (reducedP->type() == Type::Multiplication) {
      mergeMultiplicationOperands();
    }
  }
  
  void Multiplication::factorizeExponent(Expression * e1, Expression * e2, Context & context, AngleUnit angleUnit) {
    /* This function factorizes operands which share a common exponent. For
     * example, it turns Multiplication(2^x,3^x) into Multiplication(6^x). */
    assert(e1->parent() == this && e2->parent() == this);
  
    const Expression * base1 = e1->operand(0)->clone();
    const Expression * base2 = e2->operand(0);
    e2->detachOperand(base2);
    Expression * m = new Multiplication(base1, base2, false);
    removeOperand(e2, true);
    e1->replaceOperand(e1->operand(0), m, true);
  
    m->shallowReduce(context, angleUnit); // 2^x*3^x -> (2*3)^x -> 6^x
    Expression * reducedE1 = e1->shallowReduce(context, angleUnit); // 2^x*(1/2)^x -> (2*1/2)^x -> 1
    /* Reducing the new power might have turned it into a multiplication,
     * ie: 12^(1/2) -> 2*3^(1/2). In that case, we need to merge the multiplication
     * node with this. */
    if (reducedE1->type() == Type::Multiplication) {
      mergeMultiplicationOperands();
    }
  }
  
  Expression * Multiplication::distributeOnOperandAtIndex(int i, Context & context, AngleUnit angleUnit) {
    // This function turns a*(b+c) into a*b + a*c
    // We avoid deleting and creating a new addition
    Addition * a = static_cast<Addition *>(editableOperand(i));
    removeOperand(a, false);
    for (int j = 0; j < a->numberOfOperands(); j++) {
      Multiplication * m = static_cast<Multiplication *>(clone());
      Expression * termJ = a->editableOperand(j);
      a->replaceOperand(termJ, m, false);
      m->addOperand(termJ);
      m->shallowReduce(context, angleUnit); // pi^(-1)*(pi + x) -> pi^(-1)*pi + pi^(-1)*x -> 1 + pi^(-1)*x
    }
    replaceWith(a, true);
    return a->shallowReduce(context, angleUnit); // Order terms, put under a common denominator if needed
  }
  
  const Expression * Multiplication::CreateExponent(Expression * e) {
    return e->type() == Type::Power ? e->operand(1)->clone() : new Rational(1);
  }
  
  bool Multiplication::TermsHaveIdenticalBase(const Expression * e1, const Expression * e2) {
    return Base(e1)->isIdenticalTo(Base(e2));
  }
  
  bool Multiplication::TermsHaveIdenticalExponent(const Expression * e1, const Expression * e2) {
    /* Note: We will return false for e1=2 and e2=Pi, even though one could argue
     * that these have the same exponent whose value is 1. */
    return e1->type() == Type::Power && e2->type() == Type::Power && (e1->operand(1)->isIdenticalTo(e2->operand(1)));
  }
  
  bool Multiplication::TermHasRationalBase(const Expression * e) {
    return Base(e)->type() == Type::Rational;
  }
  
  bool Multiplication::TermHasRationalExponent(const Expression * e) {
    if (e->type() != Type::Power) {
      return true;
    }
    if (e->operand(1)->type() == Type::Rational) {
      return true;
    }
    return false;
  }
  
  Expression * Multiplication::shallowBeautify(Context & context, AngleUnit angleUnit) {
    /* Beautifying a Multiplication consists in several possible operations:
     * - Add Opposite ((-3)*x -> -(3*x), useful when printing fractions)
     * - Adding parenthesis if needed (a*(b+c) is not a*b+c)
     * - Creating a Division if there's either a term with a power of -1 (a.b^(-1)
     *   shall become a/b) or a non-integer rational term (3/2*a -> (3*a)/2). */
  
    // Step 1: Turn -n*A into -(n*A)
    if (operand(0)->type() == Type::Rational && operand(0)->sign() == Sign::Negative) {
      if (static_cast<const Rational *>(operand(0))->isMinusOne()) {
        removeOperand(editableOperand(0), true);
      } else {
        editableOperand(0)->setSign(Sign::Positive, context, angleUnit);
      }
      Expression * e = squashUnaryHierarchy();
      Opposite * o = new Opposite(e, true);
      e->replaceWith(o, true);
      o->editableOperand(0)->shallowBeautify(context, angleUnit);
      return o;
    }
  
    /* Step 2: Merge negative powers: a*b^(-1)*c^(-pi)*d = a*(b*c^pi)^(-1)
     * This also turns 2/3*a into 2*a*3^(-1) */
    Expression * e = mergeNegativePower(context, angleUnit);
    if (e->type() == Type::Power) {
      return e->shallowBeautify(context, angleUnit);
    }
    assert(e == this);
  
    // Step 3: Add Parenthesis if needed
    for (int i = 0; i < numberOfOperands(); i++) {
      const Expression * o = operand(i);
      if (o->type() == Type::Addition ) {
        Parenthesis * p = new Parenthesis(o, false);
        replaceOperand(o, p, false);
      }
    }
  
    // Step 4: Create a Division if needed
    for (int i = 0; i < numberOfOperands(); i++) {
      if (!(operand(i)->type() == Type::Power && operand(i)->operand(1)->type() == Type::Rational && static_cast<const Rational *>(operand(i)->operand(1))->isMinusOne())) {
        continue;
      }
  
      // Let's remove the denominator-to-be from this
      Power * p = static_cast<Power *>(editableOperand(i));
      Expression * denominatorOperand = p->editableOperand(0);
      p->detachOperand(denominatorOperand);
      removeOperand(p, true);
  
      Expression * numeratorOperand = shallowReduce(context, angleUnit);
      // Delete parenthesis unnecessary on numerator
      if (numeratorOperand->type() == Type::Parenthesis) {
        numeratorOperand = numeratorOperand->replaceWith(numeratorOperand->editableOperand(0), true);
      }
      Expression * originalParent = numeratorOperand->parent();
      Division * d = new Division(numeratorOperand, denominatorOperand, false);
      originalParent->replaceOperand(numeratorOperand, d, false);
      return d->shallowBeautify(context, angleUnit);
    }
  
    return this;
  }
  
  Expression * Multiplication::cloneDenominator(Context & context, AngleUnit angleUnit) const {
    // Merge negative power: a*b^-1*c^(-Pi)*d = a*(b*c^Pi)^-1
    // WARNING: we do not want to change the expression but to create a new one.
    SimplificationRoot root(clone());
    Expression * e = ((Multiplication *)root.operand(0))->mergeNegativePower(context, angleUnit);
    Expression * result = nullptr;
    if (e->type() == Type::Power) {
      result = static_cast<Power *>(e)->cloneDenominator(context, angleUnit);
    } else {
      assert(e->type() == Type::Multiplication);
      for (int i = 0; i < e->numberOfOperands(); i++) {
        // a*b^(-1)*... -> a*.../b
        if (e->operand(i)->type() == Type::Power && e->operand(i)->operand(1)->type() == Type::Rational && static_cast<const Rational *>(e->operand(i)->operand(1))->isMinusOne()) {
          Power * p = static_cast<Power *>(e->editableOperand(i));
          result = p->editableOperand(0);
          p->detachOperand((result));
        }
      }
    }
    root.detachOperand(e);
    delete e;
    return result;
  }
  
  Expression * Multiplication::mergeNegativePower(Context & context, AngleUnit angleUnit) {
    Multiplication * m = new Multiplication();
    // Special case for rational p/q: if q != 1, q should be at denominator
    if (operand(0)->type() == Type::Rational && !static_cast<const Rational *>(operand(0))->denominator().isOne()) {
      Rational * r = static_cast<Rational *>(editableOperand(0));
      m->addOperand(new Rational(r->denominator()));
      if (r->numerator().isOne()) {
        removeOperand(r, true);
      } else {
        replaceOperand(r, new Rational(r->numerator()), true);
      }
    }
    int i = 0;
    while (i < numberOfOperands()) {
      if (operand(i)->type() == Type::Power && operand(i)->operand(1)->sign() == Sign::Negative) {
        Expression * e = editableOperand(i);
        e->editableOperand(1)->setSign(Sign::Positive, context, angleUnit);
        removeOperand(e, false);
        m->addOperand(e);
        e->shallowReduce(context, angleUnit);
      } else {
        i++;
      }
    }
    if (m->numberOfOperands() == 0) {
      delete m;
      return this;
    }
    Power * p = new Power(m, new Rational(-1), false);
    m->sortOperands(SimplificationOrder, true);
    m->squashUnaryHierarchy();
    addOperand(p);
    sortOperands(SimplificationOrder, true);
    return squashUnaryHierarchy();
  }
  
  void Multiplication::addMissingFactors(Expression * factor, Context & context, AngleUnit angleUnit) {
    if (factor->type() == Type::Multiplication) {
      for (int j = 0; j < factor->numberOfOperands(); j++) {
        addMissingFactors(factor->editableOperand(j), context, angleUnit);
      }
      return;
    }
    /* Special case when factor is a Rational: if 'this' has already a rational
     * operand, we replace it by its LCM with factor ; otherwise, we simply add
     * factor as an operand. */
    if (numberOfOperands() > 0 && operand(0)->type() == Type::Rational && factor->type() == Type::Rational) {
      Rational * f = static_cast<Rational *>(factor);
      Rational * o = static_cast<Rational *>(editableOperand(0));
      assert(f->denominator().isOne());
      assert(o->denominator().isOne());
      Integer i = f->numerator();
      Integer j = o->numerator();
      return replaceOperand(o, new Rational(Arithmetic::LCM(&i, &j)));
    }
    if (factor->type() != Type::Rational) {
      /* If factor is not a rational, we merge it with the operand of identical
       * base if any. Otherwise, we add it as an new operand. */
      for (int i = 0; i < numberOfOperands(); i++) {
        if (TermsHaveIdenticalBase(operand(i), factor)) {
          Expression * sub = new Subtraction(CreateExponent(editableOperand(i)), CreateExponent(factor), false);
          Reduce((Expression **)&sub, context, angleUnit);
          if (sub->sign() == Sign::Negative) { // index[0] < index[1]
            if (factor->type() == Type::Power) {
              factor->replaceOperand(factor->editableOperand(1), new Opposite(sub, true), true);
            } else {
              factor = new Power(factor, new Opposite(sub, true), false);
            }
            factor->editableOperand(1)->shallowReduce(context, angleUnit);
            factorizeBase(editableOperand(i), factor, context, angleUnit);
            editableOperand(i)->shallowReduce(context, angleUnit);
          } else if (sub->sign() == Sign::Unknown) {
            factorizeBase(editableOperand(i), factor, context, angleUnit);
            editableOperand(i)->shallowReduce(context, angleUnit);
          } else {}
          delete sub;
          /* Reducing the new operand i can lead to creating a new multiplication
           * (ie 2^(1+2*3^(1/2)) -> 2*2^(2*3^(1/2)). We thus have to get rid of
           * nested multiplication: */
          mergeMultiplicationOperands();
          return;
        }
      }
    }
    addOperand(factor->clone());
    sortOperands(SimplificationOrder, false);
  }
  
  template MatrixComplex<float> Multiplication::computeOnComplexAndMatrix<float>(std::complex<float> const, const MatrixComplex<float>);
  template MatrixComplex<double> Multiplication::computeOnComplexAndMatrix<double>(std::complex<double> const, const MatrixComplex<double>);
  template std::complex<float> Multiplication::compute<float>(const std::complex<float>, const std::complex<float>);
  template std::complex<double> Multiplication::compute<double>(const std::complex<double>, const std::complex<double>);
  template void Multiplication::computeOnArrays<double>(double * m, double * n, double * result, int mNumberOfColumns, int mNumberOfRows, int nNumberOfColumns);
  
  }