#include <poincare/power.h>

#include <poincare/addition.h>
#include <poincare/arithmetic.h>
#include <poincare/binomial_coefficient.h>
#include <poincare/cosine.h>
#include <poincare/division.h>
#include <poincare/matrix.h>
#include <poincare/matrix_inverse.h>
#include <poincare/nth_root.h>
#include <poincare/opposite.h>
#include <poincare/parenthesis.h>
#include <poincare/simplification_root.h>
#include <poincare/sine.h>
#include <poincare/square_root.h>
#include <poincare/symbol.h>
#include <poincare/subtraction.h>
#include <poincare/undefined.h>

#include "layout/horizontal_layout.h"
#include "layout/vertical_offset_layout.h"

#include <cmath>
#include <math.h>
#include <ion.h>

extern "C" {
#include <assert.h>
#include <stdlib.h>
}


namespace Poincare {

Expression::Type Power::type() const {
  return Type::Power;
}

Expression * Power::clone() const {
  return new Power(m_operands, true);
}

Expression::Sign Power::sign() const {
  if (shouldStopProcessing()) {
    return Sign::Unknown;
  }
  if (operand(0)->sign() == Sign::Positive && operand(1)->sign() != Sign::Unknown) {
    return Sign::Positive;
  }
  if (operand(0)->sign() == Sign::Negative && operand(1)->type() == Type::Rational) {
    const Rational * r = static_cast<const Rational *>(operand(1));
    if (r->denominator().isOne()) {
      if (Integer::Division(r->numerator(), Integer(2)).remainder.isZero()) {
        return Sign::Positive;
      } else {
        return Sign::Negative;
      }
    }
  }
  return Sign::Unknown;
}

int Power::polynomialDegree(char symbolName) const {
  int deg = Expression::polynomialDegree(symbolName);
  if (deg == 0) {
    return deg;
  }
  int op0Deg = operand(0)->polynomialDegree(symbolName);
  if (op0Deg < 0) {
    return -1;
  }
  if (operand(1)->type() == Type::Rational) {
    const Rational * r = static_cast<const Rational *>(operand(1));
    if (!r->denominator().isOne() || r->sign() == Sign::Negative) {
      return -1;
    }
    if (Integer::NaturalOrder(r->numerator(), Integer(Integer::k_maxExtractableInteger)) > 0) {
      return -1;
    }
    op0Deg *= r->numerator().extractedInt();
    return op0Deg;
  }
  return -1;
}

int Power::privateGetPolynomialCoefficients(char symbolName, Expression * coefficients[]) const {
  int deg = polynomialDegree(symbolName);
  if (deg <= 0) {
    return Expression::privateGetPolynomialCoefficients(symbolName, coefficients);
  }
  /* Here we only consider the case x^4 as privateGetPolynomialCoefficients is
   * supposed to be called after reducing the expression. */
  if (operand(0)->type() == Type::Symbol && static_cast<const Symbol *>(operand(0))->name() == symbolName && operand(1)->type() == Type::Rational) {
    const Rational * r = static_cast<const Rational *>(operand(1));
    if (!r->denominator().isOne() || r->sign() == Sign::Negative) {
      return -1;
    }
    if (Integer::NaturalOrder(r->numerator(), Integer(Integer::k_maxExtractableInteger)) > 0) {
      return -1;
    }
    int n = r->numerator().extractedInt();
    if (n <= k_maxPolynomialDegree) {
      for (int i = 0; i < n; i++) {
        coefficients[i] = new Rational(0);
      }
      coefficients[n] = new Rational(1);
      return n;
    }
  }
  return -1;
}

Expression * Power::setSign(Sign s, Context & context, AngleUnit angleUnit) {
  assert(s == Sign::Positive);
  assert(operand(0)->sign() == Sign::Negative);
  editableOperand(0)->setSign(Sign::Positive, context, angleUnit);
  return this;
}

template<typename T>
std::complex<T> Power::compute(const std::complex<T> c, const std::complex<T> d) {
  /* Openbsd trigonometric functions are numerical implementation and thus are
   * approximative.
   * The error epsilon is ~1E-7 on float and ~1E-15 on double. In order to
   * avoid weird results as e(i*pi) = -1+6E-17*i, we compute the argument of
   * the result of c^d and if arg ~ 0 [Pi], we discard the residual imaginary
   * part and if arg ~ Pi/2 [Pi], we discard the residual real part. */
  std::complex<T> result = std::pow(c, d);
  return ApproximationEngine::truncateRealOrImaginaryPartAccordingToArgument(result);
}

template<typename T> MatrixComplex<T> Power::computeOnComplexAndMatrix(const std::complex<T> c, const MatrixComplex<T> n) {
  return MatrixComplex<T>::Undefined();
}

template<typename T> MatrixComplex<T> Power::computeOnMatrixAndComplex(const MatrixComplex<T> m, const std::complex<T> d) {
 if (m.numberOfRows() != m.numberOfColumns()) {
    return MatrixComplex<T>::Undefined();
  }
  T power = Complex<T>(d).toScalar();
  if (std::isnan(power) || std::isinf(power) || power != (int)power || std::fabs(power) > k_maxApproximatePowerMatrix) {
    return MatrixComplex<T>::Undefined();
  }
  if (power < 0) {
    MatrixComplex<T> * inverse = m.createInverse();
    if (inverse == nullptr) {
      return MatrixComplex<T>::Undefined();
    }
    Complex<T> minusC = Complex<T>(-d);
    MatrixComplex<T> result = Power::computeOnMatrixAndComplex(*inverse, minusC);
    delete inverse;
    return result;
  }
  MatrixComplex<T> result = MatrixComplex<T>::createIdentity(m.numberOfRows());
  // TODO: implement a quick exponentiation
  for (int k = 0; k < (int)power; k++) {
    if (shouldStopProcessing()) {
      return MatrixComplex<T>::Undefined();
    }
    result = Multiplication::computeOnMatrices<T>(result, m);
  }
  return result;
}

template<typename T> MatrixComplex<T> Power::computeOnMatrices(const MatrixComplex<T> m, const MatrixComplex<T> n) {
  return MatrixComplex<T>::Undefined();
}

bool Power::needParenthesisWithParent(const Expression * e) const {
  Type types[] = {Type::Power, Type::Factorial};
  return e->isOfType(types, 2);
}

ExpressionLayout * Power::createLayout(PrintFloat::Mode floatDisplayMode, int numberOfSignificantDigits) const {
  const Expression * indiceOperand = m_operands[1];
  // Delete eventual parentheses of the indice in the pretty print
  if (m_operands[1]->type() == Type::Parenthesis) {
    indiceOperand = m_operands[1]->operand(0);
  }
  HorizontalLayout * result = new HorizontalLayout();
  result->addOrMergeChildAtIndex(m_operands[0]->createLayout(floatDisplayMode, numberOfSignificantDigits), 0, false);
  result->addChildAtIndex(new VerticalOffsetLayout(
        indiceOperand->createLayout(floatDisplayMode, numberOfSignificantDigits),
        VerticalOffsetLayout::Type::Superscript,
        false),
      result->numberOfChildren());
  return result;
}

int Power::simplificationOrderSameType(const Expression * e, bool canBeInterrupted) const {
  int baseComparison = SimplificationOrder(operand(0), e->operand(0), canBeInterrupted);
  if (baseComparison != 0) {
    return baseComparison;
  }
  return SimplificationOrder(operand(1), e->operand(1), canBeInterrupted);
}

int Power::simplificationOrderGreaterType(const Expression * e, bool canBeInterrupted) const {
  int baseComparison = SimplificationOrder(operand(0), e, canBeInterrupted);
  if (baseComparison != 0) {
    return baseComparison;
  }
  Rational one(1);
  return SimplificationOrder(operand(1), &one, canBeInterrupted);
}

Expression * Power::shallowReduce(Context& context, AngleUnit angleUnit) {
  Expression * e = Expression::shallowReduce(context, angleUnit);
  if (e != this) {
    return e;
  }
#if MATRIX_EXACT_REDUCING
  /* Step 0: get rid of matrix */
  if (operand(1)->type() == Type::Matrix) {
    return replaceWith(new Undefined(), true);
  }
  if (operand(0)->type() == Type::Matrix) {
    Matrix * mat = static_cast<Matrix *>(editableOperand(0));
    if (operand(1)->type() != Type::Rational || !static_cast<const Rational *>(operand(1))->denominator().isOne()) {
      return replaceWith(new Undefined(), true);
    }
    Integer exponent = static_cast<const Rational *>(operand(1))->numerator();
    if (mat->numberOfRows() != mat->numberOfColumns()) {
      return replaceWith(new Undefined(), true);
    }
    if (exponent.isNegative()) {
      editableOperand(1)->setSign(Sign::Positive, context, angleUnit);
      Expression * newMatrix = shallowReduce(context, angleUnit);
      Expression * parent = newMatrix->parent();
      MatrixInverse * inv = new MatrixInverse(newMatrix, false);
      parent->replaceOperand(newMatrix, inv, false);
      return inv;
    }
    if (Integer::NaturalOrder(exponent, Integer(k_maxExactPowerMatrix)) > 0) {
      return this;
    }
    int exp = exponent.extractedInt(); // Ok, because 0 < exponent < k_maxExactPowerMatrix
    Matrix * id = Matrix::createIdentity(mat->numberOfRows());
    if (exp == 0) {
      return replaceWith(id, true);
    }
    Multiplication * result = new Multiplication(id, mat->clone());
    // TODO: implement a quick exponentiation
    for (int k = 1; k < exp; k++) {
      result->addOperand(mat->clone());
    }
    replaceWith(result, true);
    return result->shallowReduce(context, angleUnit);
  }
#endif

  /* Step 0: if both operands are true complexes, the result is undefined.
   * We can assert that evaluations is a complex as matrix are not simplified */
  Complex<float> * op0 = static_cast<Complex<float> *>(operand(0)->privateApproximate(float(), context, angleUnit));
  Complex<float> * op1 = static_cast<Complex<float> *>(operand(1)->privateApproximate(float(), context, angleUnit));
  bool bothOperandsComplexes = op0->imag() != 0 && op1->imag() != 0;
  bool nonComplexNegativeOperand0 = op0->imag() == 0 && op0->real() < 0;
  delete op0;
  delete op1;
  if (bothOperandsComplexes) {
    return this;
  }

  /* Step 1: We handle simple cases as x^0, x^1, 0^x and 1^x first for 2 reasons:
   * - we can assert this step that there is no division by 0:
   *   for instance, 0^(-2)->undefined
   * - we save computational time by early escaping for these cases. */
  if (operand(1)->type() == Type::Rational) {
    const Rational * b = static_cast<const Rational *>(operand(1));
    // x^0
    if (b->isZero()) {
      // 0^0 = undef
      if (operand(0)->type() == Type::Rational && static_cast<const Rational *>(operand(0))->isZero()) {
        return replaceWith(new Undefined(), true);
      }
      /* Warning: in all other case but 0^0, we replace x^0 by one. This is
       * almost always true except when x = 0. However, not substituting x^0 by
       * one would prevent from simplifying many expressions like x/x->1. */
      return replaceWith(new Rational(1), true);
    }
    // x^1
    if (b->isOne()) {
      return replaceWith(editableOperand(0), true);
    }
  }
  if (operand(0)->type() == Type::Rational) {
    Rational * a = static_cast<Rational *>(editableOperand(0));
    // 0^x
    if (a->isZero()) {
      if (operand(1)->sign() == Sign::Positive) {
        return replaceWith(new Rational(0), true);
      }
      if (operand(1)->sign() == Sign::Negative) {
        return replaceWith(new Undefined(), true);
      }
    }
    // 1^x
    if (a->isOne()) {
      return replaceWith(new Rational(1), true);
    }
  }

  /* Step 2: We look for square root and sum of square roots (two terms maximum
   * so far) at the denominator and move them to the numerator. */
  Expression * r = removeSquareRootsFromDenominator(context, angleUnit);
  if (r) {
    return r;
  }

  if (operand(1)->type() == Type::Rational) {
    const Rational * b = static_cast<const Rational *>(operand(1));
    // i^(p/q)
    if (operand(0)->type() == Type::Symbol && static_cast<const Symbol *>(operand(0))->name() == Ion::Charset::IComplex) {
      Rational r = Rational::Multiplication(*b, Rational(1, 2));
      return replaceWith(CreateNthRootOfUnity(r))->shallowReduce(context, angleUnit);
    }
  }
  bool letPowerAtRoot = parentIsALogarithmOfSameBase();
  if (operand(0)->type() == Type::Rational) {
    Rational * a = static_cast<Rational *>(editableOperand(0));
    // p^q with p, q rationals
    if (!letPowerAtRoot && operand(1)->type() == Type::Rational) {
      Rational * exp = static_cast<Rational *>(editableOperand(1));
      if (RationalExponentShouldNotBeReduced(a, exp)) {
        return this;
      }
      return simplifyRationalRationalPower(this, a, exp, context, angleUnit);
    }
  }
  // (a)^(1/2) with a < 0 --> i*(-a)^(1/2)
  if (!letPowerAtRoot && nonComplexNegativeOperand0 && operand(1)->type() == Type::Rational && static_cast<const Rational *>(operand(1))->numerator().isOne() && static_cast<const Rational *>(operand(1))->denominator().isTwo()) {
    Expression * o0 = editableOperand(0);
    Expression * m0 = new Multiplication(new Rational(-1), o0, false);
    replaceOperand(o0, m0, false);
    m0->shallowReduce(context, angleUnit);
    Multiplication * m1 = new Multiplication();
    replaceWith(m1, false);
    m1->addOperand(new Symbol(Ion::Charset::IComplex));
    m1->addOperand(this);
    shallowReduce(context, angleUnit);
    return m1->shallowReduce(context, angleUnit);
  }
  // e^(i*Pi*r) with r rational
  if (!letPowerAtRoot && isNthRootOfUnity()) {
    Expression * m = editableOperand(1);
    detachOperand(m);
    Expression * i = m->editableOperand(m->numberOfOperands()-1);
    static_cast<Multiplication *>(m)->removeOperand(i, false);
    if (angleUnit == AngleUnit::Degree) {
      const Expression * pi = m->operand(m->numberOfOperands()-1);
      m->replaceOperand(pi, new Rational(180), true);
    }
    Expression * cos = new Cosine(m, false);
    m = m->shallowReduce(context, angleUnit);
    Expression * sin = new Sine(m, true);
    Expression * complexPart = new Multiplication(sin, i, false);
    sin->shallowReduce(context, angleUnit);
    Expression * a = new Addition(cos, complexPart, false);
    cos->shallowReduce(context, angleUnit);
    complexPart->shallowReduce(context, angleUnit);
    return replaceWith(a, true)->shallowReduce(context, angleUnit);
  }
  // x^log(y,x)->y if y > 0
  if (operand(1)->type() == Type::Logarithm) {
    if (operand(1)->numberOfOperands() == 2 && operand(0)->isIdenticalTo(operand(1)->operand(1))) {
      // y > 0
      if (operand(1)->operand(0)->sign() == Sign::Positive) {
        return replaceWith(editableOperand(1)->editableOperand(0), true);
      }
    }
    // 10^log(y)
    if (operand(1)->numberOfOperands() == 1 && operand(0)->type() == Type::Rational && static_cast<const Rational *>(operand(0))->isTen()) {
      return replaceWith(editableOperand(1)->editableOperand(0), true);
    }
  }
  // (a^b)^c -> a^(b*c) if a > 0 or c is integer
  if (operand(0)->type() == Type::Power) {
    Power * p = static_cast<Power *>(editableOperand(0));
    // Check is a > 0 or c is Integer
    if (p->operand(0)->sign() == Sign::Positive ||
        (operand(1)->type() == Type::Rational && static_cast<Rational *>(editableOperand(1))->denominator().isOne())) {
      return simplifyPowerPower(p, editableOperand(1), context, angleUnit);
    }
  }
  // (a*b*c*...)^r ?
  if (!letPowerAtRoot && operand(0)->type() == Type::Multiplication) {
    Multiplication * m = static_cast<Multiplication *>(editableOperand(0));
    // (a*b*c*...)^n = a^n*b^n*c^n*... if n integer
    if (operand(1)->type() == Type::Rational && static_cast<Rational *>(editableOperand(1))->denominator().isOne()) {
      return simplifyPowerMultiplication(m, editableOperand(1), context, angleUnit);
    }
    // (a*b*...)^r -> |a|^r*(sign(a)*b*...)^r if a rational
    for (int i = 0; i < m->numberOfOperands(); i++) {
      // a is signed and a != -1
      if (m->operand(i)->sign() != Sign::Unknown && (m->operand(i)->type() != Type::Rational || !static_cast<const Rational *>(m->operand(i))->isMinusOne())) {
      //if (m->operand(i)->sign() == Sign::Positive || m->operand(i)->type() == Type::Rational) {
        Expression * r = editableOperand(1);
        Expression * rCopy = r->clone();
        Expression * factor = m->editableOperand(i);
        if (factor->sign() == Sign::Negative) {
          m->replaceOperand(factor, new Rational(-1), false);
          factor = factor->setSign(Sign::Positive, context, angleUnit);
        } else {
          m->removeOperand(factor, false);
        }
        m->shallowReduce(context, angleUnit);
        Power * p = new Power(factor, rCopy, false);
        Multiplication * root = new Multiplication(p, clone(), false);
        p->shallowReduce(context, angleUnit);
        root->editableOperand(1)->shallowReduce(context, angleUnit);
        replaceWith(root, true);
        return root->shallowReduce(context, angleUnit);
      }
    }
  }
  // a^(b+c+...) -> Rational(a^b)*a^c with a and b rational and a != 0
  if (!letPowerAtRoot && operand(0)->type() == Type::Rational && !static_cast<const Rational *>(operand(0))->isZero() && operand(1)->type() == Type::Addition) {
    Addition * a = static_cast<Addition *>(editableOperand(1));
    // Check is b is rational
    if (a->operand(0)->type() == Type::Rational) {
      const Rational * rationalBase = static_cast<const Rational *>(operand(0));
      const Rational * rationalIndex = static_cast<const Rational *>(a->operand(0));
      if (RationalExponentShouldNotBeReduced(rationalBase, rationalIndex)) {
        return this;
      }
      Power * p1 = new Power(editableOperand(0), a->editableOperand(0), true); // a^b
      Power * p2 = static_cast<Power *>(clone());
      static_cast<Addition *>(p2->editableOperand(1))->removeOperand(p2->editableOperand(1)->editableOperand(0), true); // p2 = a^(c+...)
      Multiplication * m = new Multiplication(p1, p2, false);
      simplifyRationalRationalPower(p1, static_cast<Rational *>(p1->editableOperand(0)), static_cast<Rational *>(p1->editableOperand(1)), context, angleUnit);
      replaceWith(m, true);
      return m->shallowReduce(context, angleUnit);
    }
  }

  // (a0+a1+...am)^n with n integer -> a^n+?a^(n-1)*b+?a^(n-2)*b^2+...+b^n (Multinome)
  if (!letPowerAtRoot && operand(1)->type() == Type::Rational && static_cast<const Rational *>(operand(1))->denominator().isOne() && operand(0)->type() == Type::Addition) {
    // Exponent n
    Rational * nr = static_cast<Rational *>(editableOperand(1));
    Integer n = nr->numerator();
    n.setNegative(false);
    /* if n is above 25, the resulting sum would have more than
     * k_maxNumberOfTermsInExpandedMultinome terms so we do not expand it. */
    if (Integer(k_maxNumberOfTermsInExpandedMultinome).isLowerThan(n) || n.isOne()) {
      return this;
    }
    int clippedN = n.extractedInt(); // Authorized because n < k_maxNumberOfTermsInExpandedMultinome
    // Number of terms in addition m
    int m = operand(0)->numberOfOperands();
    /* The multinome (a0+a2+...+a(m-1))^n has BinomialCoefficient(n+m-1,n) terms;
     * we expand the multinome only when the number of terms in the resulting
     * sum has less than k_maxNumberOfTermsInExpandedMultinome terms. */
    if (k_maxNumberOfTermsInExpandedMultinome < BinomialCoefficient::compute(static_cast<double>(clippedN), static_cast<double>(clippedN+m-1))) {
      return this;
    }
    Expression * result = editableOperand(0);
    Expression * a = result->clone();
    for (int i = 2; i <= clippedN; i++) {
      if (result->type() == Type::Addition) {
        Expression * a0 = new Addition();
        for (int j = 0; j < a->numberOfOperands(); j++) {
          Multiplication * m = new Multiplication(result, a->editableOperand(j), true);
          SimplificationRoot rootM(m); // m need to have a parent when applying distributeOnOperandAtIndex
          Expression * expandM = m->distributeOnOperandAtIndex(0, context, angleUnit);
          rootM.detachOperands();
          if (a0->type() == Type::Addition) {
            static_cast<Addition *>(a0)->addOperand(expandM);
          } else {
            a0 = new Addition(a0, expandM, false);
          }
          SimplificationRoot rootA0(a0); // a0 need a parent to be reduced.
          a0 = a0->shallowReduce(context, angleUnit);
          rootA0.detachOperands();
        }
        result = result->replaceWith(a0, true);
      } else {
        Multiplication * m = new Multiplication(a, result, true);
        SimplificationRoot root(m);
        result = result->replaceWith(m->distributeOnOperandAtIndex(0, context, angleUnit), true);
        result = result->shallowReduce(context, angleUnit);
        root.detachOperands();
      }
    }
    delete a;
    if (nr->sign() == Sign::Negative) {
      nr->replaceWith(new Rational(-1), true);
      return shallowReduce(context, angleUnit);
    } else {
      return replaceWith(result, true);
    }
  }
#if 0
  /* We could use the Newton formula instead which is quicker but not immediate
   * to implement in the general case (Newton multinome). */
  // (a+b)^n with n integer -> a^n+?a^(n-1)*b+?a^(n-2)*b^2+...+b^n (Newton)
  if (!letPowerAtRoot && operand(1)->type() == Type::Rational && static_cast<const Rational *>(operand(1))->denominator().isOne() && operand(0)->type() == Type::Addition && operand(0)->numberOfOperands() == 2) {
    Rational * nr = static_cast<Rational *>(editableOperand(1));
    Integer n = nr->numerator();
    n.setNegative(false);
    if (Integer(k_maxExpandedBinome).isLowerThan(n) || n.isOne()) {
      return this;
    }
    int clippedN = n.extractedInt(); // Authorized because n < k_maxExpandedBinome < k_maxNValue
    Expression * x0 = editableOperand(0)->editableOperand(0);
    Expression * x1 = editableOperand(0)->editableOperand(1);
    Addition * a = new Addition();
    for (int i = 0; i <= clippedN; i++) {
      Rational * r = new Rational(static_cast<int>(BinomialCoefficient::compute(static_cast<double>(i), static_cast<double>(clippedN))));
      Power * p0 = new Power(x0->clone(), new Rational(i), false);
      Power * p1 = new Power(x1->clone(), new Rational(clippedN-i), false);
      const Expression * operands[3] = {r, p0, p1};
      Multiplication * m = new Multiplication(operands, 3, false);
      p0->shallowReduce(context, angleUnit);
      p1->shallowReduce(context, angleUnit);
      a->addOperand(m);
      m->shallowReduce(context, angleUnit);
    }
    if (nr->sign() == Sign::Negative) {
      nr->replaceWith(new Rational(-1), true);
      editableOperand(0)->replaceWith(a, true)->shallowReduce(context, angleUnit);
      return shallowReduce(context, angleUnit);
    } else {
      return replaceWith(a, true)->shallowReduce(context, angleUnit);
    }
  }
#endif
  return this;
}

bool Power::parentIsALogarithmOfSameBase() const {
  if (parent()->type() == Type::Logarithm && parent()->operand(0) == this) {
    // parent = log(10^x)
    if (parent()->numberOfOperands() == 1) {
      if (operand(0)->type() == Type::Rational && static_cast<const Rational *>(operand(0))->isTen()) {
        return true;
      }
      return false;
    }
    // parent = log(x^y,x)
    if (operand(0)->isIdenticalTo(parent()->operand(1))) {
      return true;
    }
  }
  // parent = ln(e^x)
  if (parent()->type() == Type::NaperianLogarithm && parent()->operand(0) == this && operand(0)->type() == Type::Symbol && static_cast<const Symbol *>(operand(0))->name() == Ion::Charset::Exponential) {
    return true;
  }
  return false;
}

Expression * Power::simplifyPowerPower(Power * p, Expression * e, Context& context, AngleUnit angleUnit) {
  Expression * p0 = p->editableOperand(0);
  Expression * p1 = p->editableOperand(1);
  p->detachOperands();
  Multiplication * m = new Multiplication(p1, e, false);
  replaceOperand(e, m, false);
  replaceOperand(p, p0, true);
  m->shallowReduce(context, angleUnit);
  return shallowReduce(context, angleUnit);
}

Expression * Power::simplifyPowerMultiplication(Multiplication * m, Expression * r, Context& context, AngleUnit angleUnit) {
  for (int index = 0; index < m->numberOfOperands(); index++) {
    Expression * factor = m->editableOperand(index);
    Power * p = new Power(factor, r, true); // We copy r and factor to avoid inheritance issues
    m->replaceOperand(factor, p, true);
    p->shallowReduce(context, angleUnit);
  }
  detachOperand(m);
  return replaceWith(m, true)->shallowReduce(context, angleUnit); // delete r
}

Expression * Power::simplifyRationalRationalPower(Expression * result, Rational * a, Rational * b, Context& context, AngleUnit angleUnit) {
  if (b->denominator().isOne()) {
    Rational r = Rational::Power(*a, b->numerator());
    return result->replaceWith(new Rational(r),true);
  }
  Expression * n = nullptr;
  Expression * d = nullptr;
  if (b->sign() == Sign::Negative) {
    b->setSign(Sign::Positive);
    n = CreateSimplifiedIntegerRationalPower(a->denominator(), b, false, context, angleUnit);
    d = CreateSimplifiedIntegerRationalPower(a->numerator(), b, true, context, angleUnit);
  } else {
    n = CreateSimplifiedIntegerRationalPower(a->numerator(), b, false, context, angleUnit);
    d = CreateSimplifiedIntegerRationalPower(a->denominator(), b, true, context, angleUnit);
  }
  Multiplication * m = new Multiplication(n, d, false);
  result->replaceWith(m, true);
  return m->shallowReduce(context, angleUnit);
}

Expression * Power::CreateSimplifiedIntegerRationalPower(Integer i, Rational * r, bool isDenominator, Context & context, AngleUnit angleUnit) {
  assert(!i.isZero());
  assert(r->sign() == Sign::Positive);
  if (i.isOne()) {
    return new Rational(1);
  }
  Integer absI = i;
  absI.setNegative(false);
  Integer factors[Arithmetic::k_maxNumberOfPrimeFactors];
  Integer coefficients[Arithmetic::k_maxNumberOfPrimeFactors];
  Arithmetic::PrimeFactorization(&i, factors, coefficients, Arithmetic::k_maxNumberOfPrimeFactors);

  if (coefficients[0].isMinusOne()) {
    /* We could not break i in prime factor (either it might take too many
     * factors or too much time). */
    r->setSign(isDenominator ? Sign::Negative : Sign::Positive);
    return new Power(new Rational(i), r->clone(), false);
  }

  Integer r1(1);
  Integer r2(1);
  int index = 0;
  while (!coefficients[index].isZero() && index < Arithmetic::k_maxNumberOfPrimeFactors) {
    Integer n = Integer::Multiplication(coefficients[index], r->numerator());
    IntegerDivision div = Integer::Division(n, r->denominator());
    r1 = Integer::Multiplication(r1, Integer::Power(factors[index], div.quotient));
    r2 = Integer::Multiplication(r2, Integer::Power(factors[index], div.remainder));
    index++;
  }
  Rational * p1 = new Rational(r2);
  Integer one = isDenominator ? Integer(-1) : Integer(1);
  Rational * p2 = new Rational(one, r->denominator());
  Power * p = new Power(p1, p2, false);
  if (r1.isEqualTo(Integer(1)) && !i.isNegative()) {
    return p;
  }
  Rational * r3 = isDenominator ? new Rational(Integer(1), r1) : new Rational(r1);
  Multiplication * m = new Multiplication(r3, p, false);
  if (r2.isOne()) {
    m->removeOperand(p);
  }
  if (i.isNegative()) {
    Expression * nthRootOfUnity = CreateNthRootOfUnity(*r);
    m->addOperand(nthRootOfUnity);
    nthRootOfUnity->shallowReduce(context, angleUnit);

  }
  m->sortOperands(SimplificationOrder, false);
  return m;
}

Expression * Power::CreateNthRootOfUnity(const Rational r) {
  const Symbol * exp = new Symbol(Ion::Charset::Exponential);
  const Symbol * iComplex = new Symbol(Ion::Charset::IComplex);
  const Symbol * pi = new Symbol(Ion::Charset::SmallPi);
  const Expression * multExpOperands[3] = {iComplex, pi, new Rational(r)};
  Multiplication * mExp = new Multiplication(multExpOperands, 3, false);
  mExp->sortOperands(SimplificationOrder, false);
  return new Power(exp, mExp, false);
#if 0
  const Symbol * iComplex = new Symbol(Ion::Charset::IComplex);
  const Symbol * pi = new Symbol(Ion::Charset::SmallPi);
  Expression * op = new Multiplication(pi, r->clone(), false);
  Cosine * cos = new Cosine(op, false);
  op = op->shallowReduce(context, angleUnit);
  Sine * sin = new Sine(op, true);
  Expression * m = new Multiplication(iComplex, sin, false);
  sin->shallowReduce(context, angleUnit);
  Expression * a = new Addition(cos, m, false);
  cos->shallowReduce(context, angleUnit);
  const Expression * multExpOperands[3] = {pi, r->clone()};
#endif
}

Expression * Power::shallowBeautify(Context& context, AngleUnit angleUnit) {
  // X^-y -> 1/(X->shallowBeautify)^y
  if (operand(1)->sign() == Sign::Negative) {
    Expression * p = cloneDenominator(context, angleUnit);
    Division * d = new Division(new Rational(1), p, false);
    p->shallowReduce(context, angleUnit);
    replaceWith(d, true);
    return d->shallowBeautify(context, angleUnit);
  }
  if (operand(1)->type() == Type::Rational && static_cast<const Rational *>(operand(1))->numerator().isOne()) {
    Integer index = static_cast<const Rational *>(operand(1))->denominator();
    if (index.isEqualTo(Integer(2))) {
      const Expression * sqrtOperand[1] = {operand(0)};
      SquareRoot * sqr = new SquareRoot(sqrtOperand, true);
      return replaceWith(sqr, true);
    }
    const Expression * rootOperand[2] = {operand(0)->clone(), new Rational(index)};
    NthRoot * nr = new NthRoot(rootOperand, false);
    return replaceWith(nr, true);
  }
  // +(a,b)^c ->(+(a,b))^c
  if (operand(0)->type() == Type::Addition || operand(0)->type() == Type::Multiplication) {
    const Expression * o[1] = {operand(0)};
    Parenthesis * p = new Parenthesis(o, true);
    replaceOperand(operand(0), p, true);
  }
  return this;
}

Expression * Power::cloneDenominator(Context & context, AngleUnit angleUnit) const {
  if (operand(1)->sign() == Sign::Negative) {
    Expression * denominator = clone();
    Expression * newExponent = denominator->editableOperand(1)->setSign(Sign::Positive, context, angleUnit);
    if (newExponent->type() == Type::Rational && static_cast<Rational *>(newExponent)->isOne()) {
      delete denominator;
      return operand(0)->clone();
    }
    return denominator;
  }
  return nullptr;
}

bool Power::TermIsARationalSquareRootOrRational(const Expression * e) {
  if (e->type() == Type::Rational) {
    return true;
  }
  if (e->type() == Type::Power && e->operand(0)->type() == Type::Rational && e->operand(1)->type() == Type::Rational && static_cast<const Rational *>(e->operand(1))->isHalf()) {
    return true;
  }
  if (e->type() == Type::Multiplication && e->numberOfOperands() == 2 && e->operand(0)->type() == Type::Rational && e->operand(1)->type() == Type::Power && e->operand(1)->operand(0)->type() == Type::Rational && e->operand(1)->operand(1)->type() == Type::Rational && static_cast<const Rational *>(e->operand(1)->operand(1))->isHalf()) {
  return true;
  }
  return false;
}

const Rational * Power::RadicandInExpression(const Expression * e) {
  if (e->type() == Type::Rational) {
    return nullptr;
  } else if (e->type() == Type::Power) {
    assert(e->type() == Type::Power);
    assert(e->operand(0)->type() == Type::Rational);
    return static_cast<const Rational *>(e->operand(0));
  } else {
    assert(e->type() == Type::Multiplication);
    assert(e->operand(1)->type() == Type::Power);
    assert(e->operand(1)->operand(0)->type() == Type::Rational);
    return static_cast<const Rational *>(e->operand(1)->operand(0));
  }
}

const Rational * Power::RationalFactorInExpression(const Expression * e) {
  if (e->type() == Type::Rational) {
    return static_cast<const Rational *>(e);
  } else if (e->type() == Type::Power) {
    return nullptr;
  } else {
    assert(e->type() == Type::Multiplication);
    assert(e->operand(0)->type() == Type::Rational);
    return static_cast<const Rational *>(e->operand(0));
  }
}

Expression * Power::removeSquareRootsFromDenominator(Context & context, AngleUnit angleUnit) {
  Expression * result = nullptr;

  if (operand(0)->type() == Type::Rational && operand(1)->type() == Type::Rational && (static_cast<const Rational *>(operand(1))->isHalf() || static_cast<const Rational *>(operand(1))->isMinusHalf())) {
      /* We're considering a term of the form sqrt(p/q) (or 1/sqrt(p/q)), with
       * p and q integers.
       * We'll turn those into sqrt(p*q)/q (or sqrt(p*q)/p) . */
      Integer p = static_cast<const Rational *>(operand(0))->numerator();
      assert(!p.isZero()); // We eliminated case of form 0^(-1/2) at first step of shallowReduce
      Integer q = static_cast<const Rational *>(operand(0))->denominator();
      // We do nothing for terms of the form sqrt(p)
      if (!q.isOne() || static_cast<const Rational *>(operand(1))->isMinusHalf()) {
        Power * sqrt = new Power(new Rational(Integer::Multiplication(p, q)), new Rational(1, 2), false);
        if (static_cast<const Rational *>(operand(1))->isHalf()) {
          result = new Multiplication(new Rational(Integer(1), q), sqrt, false);
        } else {
          result = new Multiplication(new Rational(Integer(1), p), sqrt, false); // We use here the assertion that p != 0
        }
        sqrt->shallowReduce(context, angleUnit);
      }
  } else if (operand(1)->type() == Type::Rational && static_cast<const Rational *>(operand(1))->isMinusOne() && operand(0)->type() == Type::Addition && operand(0)->numberOfOperands() == 2 && TermIsARationalSquareRootOrRational(operand(0)->operand(0)) && TermIsARationalSquareRootOrRational(operand(0)->operand(1))) {
    /* We're considering a term of the form
     *
     * 1/(n1/d1*sqrt(p1/q1) + n2/d2*sqrt(p2/q2))
     *
     * and we want to turn it into
     *
     *  n1*q2*d1*d2^2*sqrt(p1*q1) - n2*q1*d2*d1^2*sqrt(p2*q2)
     * -------------------------------------------------------
     *          n1^2*d2^2*p1*q2 - n2^2*d1^2*p2*q1
     */
    const Rational * f1 = RationalFactorInExpression(operand(0)->operand(0));
    const Rational * f2 = RationalFactorInExpression(operand(0)->operand(1));
    const Rational * r1 = RadicandInExpression(operand(0)->operand(0));
    const Rational * r2 = RadicandInExpression(operand(0)->operand(1));
    Integer n1 = (f1 ? f1->numerator() : Integer(1));
    Integer d1 = (f1 ? f1->denominator() : Integer(1));
    Integer p1 = (r1 ? r1->numerator() : Integer(1));
    Integer q1 = (r1 ? r1->denominator() : Integer(1));
    Integer n2 = (f2 ? f2->numerator() : Integer(1));
    Integer d2 = (f2 ? f2->denominator() : Integer(1));
    Integer p2 = (r2 ? r2->numerator() : Integer(1));
    Integer q2 = (r2 ? r2->denominator() : Integer(1));

    // Compute the denominator = n1^2*d2^2*p1*q2 - n2^2*d1^2*p2*q1
    Integer denominator = Integer::Subtraction(
        Integer::Multiplication(
          Integer::Multiplication(
            Integer::Power(n1, Integer(2)),
            Integer::Power(d2, Integer(2))),
          Integer::Multiplication(p1, q2)),
        Integer::Multiplication(
          Integer::Multiplication(
            Integer::Power(n2, Integer(2)),
            Integer::Power(d1, Integer(2))),
          Integer::Multiplication(p2, q1)));

    // Compute the numerator
    Power * sqrt1 = new Power(new Rational(Integer::Multiplication(p1, q1)), new Rational(1, 2), false);
    Power * sqrt2 = new Power(new Rational(Integer::Multiplication(p2, q2)), new Rational(1, 2), false);
    Integer factor1 = Integer::Multiplication(
        Integer::Multiplication(n1, d1),
        Integer::Multiplication(Integer::Power(d2, Integer(2)), q2));
    Multiplication * m1 = new Multiplication(new Rational(factor1), sqrt1, false);
    Integer factor2 = Integer::Multiplication(
        Integer::Multiplication(n2, d2),
        Integer::Multiplication(Integer::Power(d1, Integer(2)), q1));
    Multiplication * m2 = new Multiplication(new Rational(factor2), sqrt2, false);
    Subtraction * numerator = nullptr;
    if (denominator.isNegative()) {
      numerator = new Subtraction(m2, m1, false);
      denominator.setNegative(false);
    } else {
      numerator = new Subtraction(m1, m2, false);
    }

    result = new Multiplication(numerator, new Rational(Integer(1), denominator), false);
    numerator->deepReduce(context, angleUnit);
  }

  if (result) {
    replaceWith(result, true);
    result = result->shallowReduce(context, angleUnit);
  }
  return result;
}

bool Power::isNthRootOfUnity() const {
  if (operand(0)->type() != Type::Symbol || static_cast<const Symbol *>(operand(0))->name() != Ion::Charset::Exponential) {
    return false;
  }
  if (operand(1)->type() != Type::Multiplication) {
    return false;
  }
  if (operand(1)->numberOfOperands() < 2 || operand(1)->numberOfOperands() > 3) {
    return false;
  }
  const Expression * i = operand(1)->operand(operand(1)->numberOfOperands()-1);
  if (i->type() != Type::Symbol || static_cast<const Symbol *>(i)->name() != Ion::Charset::IComplex) {
    return false;
  }
  const Expression * pi = operand(1)->operand(operand(1)->numberOfOperands()-2);
  if (pi->type() != Type::Symbol || static_cast<const Symbol *>(pi)->name() != Ion::Charset::SmallPi) {
    return false;
  }
  if (numberOfOperands() == 2) {
    return true;
  }
  if (operand(1)->operand(0)->type() == Type::Rational) {
    return true;
  }
  return false;
}

bool Power::RationalExponentShouldNotBeReduced(const Rational * b, const Rational * r) {
  if (r->isMinusOne()) {
    return false;
  }
  /* We check that the simplification does not involve too complex power of
   * integers (ie 3^999, 120232323232^50) that would take too much time to
   * compute:
   *  - we cap the exponent at k_maxExactPowerMatrix
   *  - we cap the resulting power at DBL_MAX
   * The complexity of computing a power of rational is mainly due to computing
   * the GCD of the resulting numerator and denominator. Euclide algorithm's
   * complexity is apportionned to the number of decimal digits in the smallest
   * integer. */
  Integer maxIntegerExponent = r->numerator();
  maxIntegerExponent.setNegative(false);
  if (Integer::NaturalOrder(maxIntegerExponent, Integer(k_maxExactPowerMatrix)) > 0) {
    return true;
  }
  double index = maxIntegerExponent.approximate<double>();
  double powerNumerator = std::pow(std::fabs(b->numerator().approximate<double>()), index);
  double powerDenominator = std::pow(std::fabs(b->denominator().approximate<double>()), index);
  if (std::isnan(powerNumerator) || std::isnan(powerDenominator) || std::isinf(powerNumerator) || std::isinf(powerDenominator)) {
    return true;
  }
  return false;
}

template std::complex<float> Power::compute<float>(std::complex<float>, std::complex<float>);
template std::complex<double> Power::compute<double>(std::complex<double>, std::complex<double>);
}