From f15a43531f52ea428d0dafbad2d4ae57cff77ad3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9a=20Saviot?= Date: Thu, 9 Aug 2018 15:40:44 +0200 Subject: [PATCH] [poincare] Power --- poincare/Makefile | 1 + .../allocation_failure_evaluation_node.h | 2 - poincare/include/poincare/complex.h | 2 - poincare/include/poincare/evaluation.h | 4 - poincare/include/poincare/expression.h | 3 +- poincare/include/poincare/expression_node.h | 3 +- poincare/include/poincare/matrix_complex.h | 5 +- poincare/include/poincare/power.h | 35 +- poincare/src/division.cpp | 16 +- poincare/src/matrix_complex.cpp | 4 +- poincare/src/power.cpp | 1097 +++++++++-------- 11 files changed, 609 insertions(+), 563 deletions(-) diff --git a/poincare/Makefile b/poincare/Makefile index 4c9a8dc0c..9334c8396 100644 --- a/poincare/Makefile +++ b/poincare/Makefile @@ -6,6 +6,7 @@ objs += $(addprefix poincare/src/,\ addition.o\ matrix_complex.o\ multiplication.o\ + power.o\ tree_node.o\ tree_pool.o\ tree_by_reference.o\ diff --git a/poincare/include/poincare/allocation_failure_evaluation_node.h b/poincare/include/poincare/allocation_failure_evaluation_node.h index 27fe387fd..d9a7269e9 100644 --- a/poincare/include/poincare/allocation_failure_evaluation_node.h +++ b/poincare/include/poincare/allocation_failure_evaluation_node.h @@ -17,8 +17,6 @@ public: Expression complexToExpression(Preferences::Preferences::ComplexFormat complexFormat) const override { return Undefined(); } std::complex trace() const override { return std::complex(NAN); } std::complex determinant() const override { return std::complex(NAN); } - Evaluation inverse() const override { return Complex::Undefined(); } - virtual Evaluation transpose() const override { return Complex::Undefined(); } // TreeNode size_t size() const override { return sizeof(AllocationFailureEvaluationNode); } #if TREE_LOG diff --git a/poincare/include/poincare/complex.h b/poincare/include/poincare/complex.h index a60dff243..2fc873ea6 100644 --- a/poincare/include/poincare/complex.h +++ b/poincare/include/poincare/complex.h @@ -27,8 +27,6 @@ public: Expression complexToExpression(Preferences::Preferences::ComplexFormat complexFormat) const override; std::complex trace() const override { return *this; } std::complex determinant() const override { return *this; } - Evaluation inverse() const override; - Evaluation transpose() const override { return Complex(*this); } }; template diff --git a/poincare/include/poincare/evaluation.h b/poincare/include/poincare/evaluation.h index cd21dad3a..5a5df3307 100644 --- a/poincare/include/poincare/evaluation.h +++ b/poincare/include/poincare/evaluation.h @@ -31,8 +31,6 @@ public: virtual Expression complexToExpression(Preferences::ComplexFormat complexFormat) const = 0; virtual std::complex trace() const = 0; virtual std::complex determinant() const = 0; - virtual Evaluation inverse() const = 0; - virtual Evaluation transpose() const = 0; // TreeNode static TreeNode * FailedAllocationStaticNode(); @@ -56,8 +54,6 @@ public: Expression complexToExpression(Preferences::ComplexFormat complexFormat) const; std::complex trace() const { return node()->trace(); } std::complex determinant() const { return node()->determinant(); } - Evaluation inverse() const { return node()->inverse(); } - Evaluation transpose() const { return node()->transpose(); } protected: Evaluation(EvaluationNode * n) : TreeByValue(n) {} //Evaluation() : TreeByValue() {} diff --git a/poincare/include/poincare/expression.h b/poincare/include/poincare/expression.h index 257b0c245..fd257282f 100644 --- a/poincare/include/poincare/expression.h +++ b/poincare/include/poincare/expression.h @@ -23,6 +23,7 @@ class Expression : public TreeByValue { friend class Division; friend class NAryExpressionNode; friend class Parenthesis; + friend class Power; public: /* Constructor & Destructor */ Expression() : Expression(nullptr) {} @@ -133,7 +134,7 @@ protected: private: /* Properties */ - int getPolynomialCoefficients(char symbolName, Expression coefficients[], Context & context, Preferences::AngleUnit angleUnit) const; + int getPolynomialCoefficients(char symbolName, Expression coefficients[]) const; /* Simplification */ Expression shallowReduce(Context & context, Preferences::AngleUnit angleUnit) const; diff --git a/poincare/include/poincare/expression_node.h b/poincare/include/poincare/expression_node.h index 9707df3fd..34efb2a92 100644 --- a/poincare/include/poincare/expression_node.h +++ b/poincare/include/poincare/expression_node.h @@ -14,9 +14,10 @@ namespace Poincare { * 'this' outdated. They should only be called in a wrapper on Expression. */ class ExpressionNode : public TreeNode, public SerializationHelperInterface { - friend class SymbolNode; friend class DivisionNode; friend class NAryExpressionNode; + friend class PowerNode; + friend class SymbolNode; public: enum class Type : uint8_t { AllocationFailure = 0, diff --git a/poincare/include/poincare/matrix_complex.h b/poincare/include/poincare/matrix_complex.h index 94e8d84c7..d09af2520 100644 --- a/poincare/include/poincare/matrix_complex.h +++ b/poincare/include/poincare/matrix_complex.h @@ -41,8 +41,8 @@ public: Expression complexToExpression(Preferences::Preferences::ComplexFormat complexFormat) const override; std::complex trace() const override; std::complex determinant() const override; - Evaluation inverse() const override; - Evaluation transpose() const override; + MatrixComplex inverse() const; + MatrixComplex transpose() const; private: int m_numberOfRows; int m_numberOfColumns; @@ -65,6 +65,7 @@ public: return MatrixComplex((std::complex *)&undef, 1, 1); } static MatrixComplex createIdentity(int dim); + MatrixComplex inverse() const { return node()->inverse(); } Complex complexAtIndex(int index) { return Complex(node()->complexAtIndex(index)); } diff --git a/poincare/include/poincare/power.h b/poincare/include/poincare/power.h index 1e1a621a7..260eff8b7 100644 --- a/poincare/include/poincare/power.h +++ b/poincare/include/poincare/power.h @@ -2,9 +2,9 @@ #define POINCARE_POWER_H #include -#include #include -#include +#include +#include namespace Poincare { @@ -19,14 +19,14 @@ public: // TreeNode size_t size() const override { return sizeof(PowerNode); } #if TREE_LOG - const char * description() const override { return "Division"; } + const char * description() const override { return "Power"; } #endif int numberOfChildren() const override { return 2; } // Properties - virtual Type type() const override { return Type::Power; } - virtual Sign sign() const override; - virtual Expression setSign(Sign s, Context & context, Preferences::AngleUnit angleUnit) const override; + Type type() const override { return Type::Power; } + Sign sign() const override; + Expression setSign(Sign s, Context & context, Preferences::AngleUnit angleUnit) const override; int polynomialDegree(char symbolName) const override; int getPolynomialCoefficients(char symbolName, Expression coefficients[]) const override; @@ -36,23 +36,25 @@ private: constexpr static int k_maxNumberOfTermsInExpandedMultinome = 25; constexpr static int k_maxExactPowerMatrix = 100; - /* Layout */ + // Layout LayoutRef createLayout(Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const override; - /* Serialize */ + // Serialize bool needsParenthesesWithParent(const SerializationHelperInterface * parent) const override; int serialize(char * buffer, int bufferSize, Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const override { return SerializationHelper::Infix(this, buffer, bufferSize, floatDisplayMode, numberOfSignificantDigits, name()); } static const char * name() { return "^"; } + /* Simplify */ Expression shallowReduce(Context& context, Preferences::AngleUnit angleUnit) const override; - Expression * shallowBeautify(Context & context, Preferences::AngleUnit angleUnit) const override; - int simplificationOrderGreaterType(const Expression * e, bool canBeInterrupted) const override; + Expression shallowBeautify(Context & context, Preferences::AngleUnit angleUnit) const override; + int simplificationOrderGreaterType(const ExpressionNode * e, bool canBeInterrupted) const override; int simplificationOrderSameType(const ExpressionNode * e, bool canBeInterrupted) const override; - Expression * simplifyPowerPower(Power * p, Expression * r, Context & context, Preferences::AngleUnit angleUnit); - Expression * denominator(Context & context, Preferences::AngleUnit angleUnit) const override; - Expression * simplifyPowerMultiplication(Multiplication * m, Expression * r, Context & context, Preferences::AngleUnit angleUnit); + Expression denominator(Context & context, Preferences::AngleUnit angleUnit) const override; +#if 0 + Expression * simplifyPowerPower(Power * p, Expression * r, Context & context, Preferences::AngleUnit angleUnit); + Expression * simplifyPowerMultiplication(Multiplication * m, Expression * r, Context & context, Preferences::AngleUnit angleUnit); Expression * simplifyRationalRationalPower(Expression * result, Rational * a, Rational * b, Context & context, Preferences::AngleUnit angleUnit); Expression * removeSquareRootsFromDenominator(Context & context, Preferences::AngleUnit angleUnit); bool parentIsALogarithmOfSameBase() const; @@ -63,6 +65,7 @@ private: static const Rational * RadicandInExpression(const Expression * e); static const Rational * RationalFactorInExpression(const Expression * e); static bool RationalExponentShouldNotBeReduced(const Rational * b, const Rational * r); +#endif /* Evaluation */ constexpr static int k_maxApproximatePowerMatrix = 1000; template static MatrixComplex computeOnComplexAndMatrix(const std::complex c, const MatrixComplex n); @@ -80,10 +83,10 @@ class Power : public Expression { public: Power(Expression base, Expression exponent); Power(const PowerNode * n) : Expression(n) {} - Expression setSign(Sign s, Context & context, Preferences::AngleUnit angleUnit) const; - Expression shallowReduce(Context & context, Preferences::AngleUnit angleUnit) const; - Expression * shallowBeautify(Context & context, Preferences::AngleUnit angleUnit) const; + Expression setSign(ExpressionNode::Sign s, Context & context, Preferences::AngleUnit angleUnit) const; int getPolynomialCoefficients(char symbolName, Expression coefficients[]) const; + Expression shallowReduce(Context & context, Preferences::AngleUnit angleUnit) const; + Expression shallowBeautify(Context & context, Preferences::AngleUnit angleUnit) const; }; diff --git a/poincare/src/division.cpp b/poincare/src/division.cpp index 57dd79dda..fa5ea3b3b 100644 --- a/poincare/src/division.cpp +++ b/poincare/src/division.cpp @@ -47,12 +47,8 @@ template Complex DivisionNode::compute(const std::complex c, c } template MatrixComplex DivisionNode::computeOnComplexAndMatrix(const std::complex c, const MatrixComplex n) { - Evaluation inverse = n.inverse(); - if (inverse.isUndefined()) { - return MatrixComplex::Undefined(); - } - assert(inverse.node()->type() == EvaluationNode::Type::MatrixComplex || inverse.node()->type() == EvaluationNode::Type::AllocationFailure); - MatrixComplex result = MultiplicationNode::computeOnComplexAndMatrix(c, MatrixComplex(static_cast *>(inverse.node()))); + MatrixComplex inverse = n.inverse(); + MatrixComplex result = MultiplicationNode::computeOnComplexAndMatrix(c, inverse); return result; } @@ -60,12 +56,8 @@ template MatrixComplex DivisionNode::computeOnMatrices(const Matr if (m.numberOfColumns() != n.numberOfColumns()) { return MatrixComplex::Undefined(); } - Evaluation inverse = n.inverse(); - if (inverse.isUndefined()) { - return MatrixComplex::Undefined(); - } - assert(inverse.node()->type() == EvaluationNode::Type::MatrixComplex || inverse.node()->type() == EvaluationNode::Type::AllocationFailure); - MatrixComplex result = MultiplicationNode::computeOnMatrices(m, MatrixComplex(static_cast *>(inverse.node()))); + MatrixComplex inverse = n.inverse(); + MatrixComplex result = MultiplicationNode::computeOnMatrices(m, inverse); return result; } diff --git a/poincare/src/matrix_complex.cpp b/poincare/src/matrix_complex.cpp index 74d209fa2..737adad5e 100644 --- a/poincare/src/matrix_complex.cpp +++ b/poincare/src/matrix_complex.cpp @@ -89,7 +89,7 @@ std::complex MatrixComplexNode::determinant() const { } template -Evaluation MatrixComplexNode::inverse() const { +MatrixComplex MatrixComplexNode::inverse() const { /* TODO if (numberOfRows() != numberOfColumns() || numberOfChildren() == 0 || numberOfChildren() > Matrix::k_maxNumberOfCoefficients) { return MatrixComplex::Undefined(); } @@ -110,7 +110,7 @@ Evaluation MatrixComplexNode::inverse() const { } template -Evaluation MatrixComplexNode::transpose() const { +MatrixComplex MatrixComplexNode::transpose() const { // Intentionally swapping dimensions for transpose MatrixComplex result; for (int i = 0; i < numberOfRows(); i++) { diff --git a/poincare/src/power.cpp b/poincare/src/power.cpp index 44cbcbbac..a5919efda 100644 --- a/poincare/src/power.cpp +++ b/poincare/src/power.cpp @@ -1,18 +1,18 @@ #include #include -#include -#include -#include +//#include +//#include +//#include #include -#include -#include -#include +//#include +//#include +//#include #include #include -#include -#include -#include +//#include +//#include +//#include #include #include #include @@ -26,22 +26,27 @@ extern "C" { #include -#include } namespace Poincare { -Expression::Sign Power::sign() const { - if (shouldStopProcessing()) { +PowerNode * PowerNode::FailedAllocationStaticNode() { + static AllocationFailureExpressionNode failure; + return &failure; +} + +ExpressionNode::Sign PowerNode::sign() const { + if (Expression::shouldStopProcessing()) { return Sign::Unknown; } - if (operand(0)->sign() == Sign::Positive && operand(1)->sign() != Sign::Unknown) { + if (childAtIndex(0)->sign() == Sign::Positive && childAtIndex(1)->sign() != Sign::Unknown) { return Sign::Positive; } - if (operand(0)->sign() == Sign::Negative && operand(1)->type() == Type::Rational) { - const Rational * r = static_cast(operand(1)); + if (childAtIndex(0)->sign() == Sign::Negative && childAtIndex(1)->type() == Type::Rational) { + RationalNode * r = static_cast(childAtIndex(1)); if (r->denominator().isOne()) { - if (Integer::Division(r->numerator(), Integer(2)).remainder.isZero()) { + NaturalIntegerPointer nip = r->numerator(); + if (Integer::Division(Integer(&nip), Integer(2)).remainder.isZero()) { return Sign::Positive; } else { return Sign::Negative; @@ -51,68 +56,43 @@ Expression::Sign Power::sign() const { return Sign::Unknown; } -int Power::polynomialDegree(char symbolName) const { - int deg = Expression::polynomialDegree(symbolName); +Expression PowerNode::setSign(Sign s, Context & context, Preferences::AngleUnit angleUnit) const { + return Power(this).setSign(s, context, angleUnit); +} + +int PowerNode::polynomialDegree(char symbolName) const { + int deg = Power(this).polynomialDegree(symbolName); if (deg == 0) { return deg; } - int op0Deg = operand(0)->polynomialDegree(symbolName); + int op0Deg = childAtIndex(0)->polynomialDegree(symbolName); if (op0Deg < 0) { return -1; } - if (operand(1)->type() == Type::Rational) { - const Rational * r = static_cast(operand(1)); + if (childAtIndex(1)->type() == Type::Rational) { + RationalNode * r = static_cast(childAtIndex(1)); if (!r->denominator().isOne() || r->sign() == Sign::Negative) { return -1; } - if (Integer::NaturalOrder(r->numerator(), Integer(Integer::k_maxExtractableInteger)) > 0) { + NaturalIntegerPointer nip = r->numerator(); + Integer numeratorInt = Integer(&nip); + if (Integer::NaturalOrder(numeratorInt, Integer(Integer::k_maxExtractableInteger)) > 0) { return -1; } - op0Deg *= r->numerator().extractedInt(); + op0Deg *= numeratorInt.extractedInt(); return op0Deg; } return -1; } -int Power::getPolynomialCoefficients(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(operand(0))->name() == symbolName && operand(1)->type() == Type::Rational) { - const Rational * r = static_cast(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] = RationalReference(0); - } - coefficients[n] = RationalReference(1); - return n; - } - } - return -1; +int PowerNode::getPolynomialCoefficients(char symbolName, Expression coefficients[]) const { + return Power(this).getPolynomialCoefficients(symbolName, coefficients); } -Expression PowerNode::setSign(Sign s, Context & context, Preferences::AngleUnit angleUnit) const { - return Power(this).setSign(s, context, angleUnit); -} - -Expression Power::setSign(Sign s, Context & context, Preferences::AngleUnit angleUnit) const { - assert(s == Sign::Positive); - assert(childAtIndex(0).sign() == Sign::Negative); - return Power(childAtIndex(0).setSign(Sign::Positive, context, angleUnit), childAtIndex(1)); -} +// Private template -std::complex Power::compute(const std::complex c, const std::complex d) { +Complex PowerNode::compute(const std::complex c, const std::complex 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 @@ -120,59 +100,19 @@ std::complex Power::compute(const std::complex c, const std::complex d) * 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 result = std::pow(c, d); - return ApproximationHelper::TruncateRealOrImaginaryPartAccordingToArgument(result); + return Complex(ApproximationHelper::TruncateRealOrImaginaryPartAccordingToArgument(result)); } -template MatrixComplex Power::computeOnComplexAndMatrix(const std::complex c, const MatrixComplex n) { - return MatrixComplex::Undefined(); -} +// Layout -template MatrixComplex Power::computeOnMatrixAndComplex(const MatrixComplex m, const std::complex d) { - if (m.numberOfRows() != m.numberOfColumns()) { - return MatrixComplex::Undefined(); - } - T power = Complex(d).toScalar(); - if (std::isnan(power) || std::isinf(power) || power != (int)power || std::fabs(power) > k_maxApproximatePowerMatrix) { - return MatrixComplex::Undefined(); - } - if (power < 0) { - MatrixComplex * inverse = m.createInverse(); - if (inverse == nullptr) { - return MatrixComplex::Undefined(); - } - Complex minusC = Complex(-d); - MatrixComplex result = Power::computeOnMatrixAndComplex(*inverse, minusC); - delete inverse; - return result; - } - MatrixComplex result = MatrixComplex::createIdentity(m.numberOfRows()); - // TODO: implement a quick exponentiation - for (int k = 0; k < (int)power; k++) { - if (shouldStopProcessing()) { - return MatrixComplex::Undefined(); - } - result = Multiplication::computeOnMatrices(result, m); - } - return result; -} - -template MatrixComplex Power::computeOnMatrices(const MatrixComplex m, const MatrixComplex n) { - return MatrixComplex::Undefined(); -} - -bool Power::needsParenthesesWithParent(const SerializationHelperInterface * e) const { - Type types[] = {Type::Power, Type::Factorial}; - return e->isOfType(types, 2); -} - -LayoutRef Power::createLayout(Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const { - const Expression * indiceOperand = m_operands[1]; +LayoutRef PowerNode::createLayout(Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const { + ExpressionNode * indiceOperand = childAtIndex(1); // Delete eventual parentheses of the indice in the pretty print - if (m_operands[1]->type() == Type::Parenthesis) { - indiceOperand = m_operands[1]->operand(0); + if (indiceOperand->type() == Type::Parenthesis) { + indiceOperand = indiceOperand->childAtIndex(0); } HorizontalLayoutRef result = HorizontalLayoutRef(); - result.addOrMergeChildAtIndex(m_operands[0]->createLayout(floatDisplayMode, numberOfSignificantDigits), 0, false); + result.addOrMergeChildAtIndex(childAtIndex(0)->createLayout(floatDisplayMode, numberOfSignificantDigits), 0, false); result.addChildAtIndex(VerticalOffsetLayoutRef( indiceOperand->createLayout(floatDisplayMode, numberOfSignificantDigits), VerticalOffsetLayoutNode::Type::Superscript), @@ -182,367 +122,59 @@ LayoutRef Power::createLayout(Preferences::PrintFloatMode floatDisplayMode, int return result; } -int Power::simplificationOrderSameType(const ExpressionNode * 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); +// Serialize + +bool PowerNode::needsParenthesesWithParent(const SerializationHelperInterface * e) const { + Type types[] = {Type::Power, Type::Factorial}; + return static_cast(e)->isOfType(types, 2); } -int Power::simplificationOrderGreaterType(const Expression * e, bool canBeInterrupted) const { - int baseComparison = SimplificationOrder(operand(0), e, canBeInterrupted); +// Simplify + +Expression PowerNode::shallowReduce(Context& context, Preferences::AngleUnit angleUnit) const { + return Power(this).shallowReduce(context, angleUnit); +} + +Expression PowerNode::shallowBeautify(Context& context, Preferences::AngleUnit angleUnit) const { + return Power(this).shallowBeautify(context, angleUnit); +} + + +int PowerNode::simplificationOrderGreaterType(const ExpressionNode * e, bool canBeInterrupted) const { + int baseComparison = SimplificationOrder(childAtIndex(0), e, canBeInterrupted); if (baseComparison != 0) { return baseComparison; } Rational one(1); - return SimplificationOrder(operand(1), &one, canBeInterrupted); + return SimplificationOrder(childAtIndex(1), one.node(), canBeInterrupted); } -Expression Power::shallowReduce(Context& context, Preferences::AngleUnit angleUnit) { - Expression * e = Expression::shallowReduce(context, angleUnit); - if (e != this) { - return e; +int PowerNode::simplificationOrderSameType(const ExpressionNode * e, bool canBeInterrupted) const { + assert(e->numberOfChildren() > 0); + int baseComparison = SimplificationOrder(childAtIndex(0), e->childAtIndex(0), canBeInterrupted); + if (baseComparison != 0) { + return baseComparison; } -#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(editableOperand(0)); - if (operand(1)->type() != Type::Rational || !static_cast(operand(1))->denominator().isOne()) { - return replaceWith(new Undefined(), true); - } - Integer exponent = static_cast(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 + assert(e->numberOfChildren() > 1); + return SimplificationOrder(childAtIndex(1), e->childAtIndex(1), canBeInterrupted); +} - /* 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 * op0 = static_cast *>(operand(0)->privateApproximate(float(), context, angleUnit)); - Complex * op1 = static_cast *>(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; +Expression PowerNode::denominator(Context & context, Preferences::AngleUnit angleUnit) const { + if (childAtIndex(1)->sign() == Sign::Negative) { + Expression denominator = Power(this); + Expression newExponent = denominator.childAtIndex(1).setSign(Sign::Positive, context, angleUnit); + if (newExponent.type() == Type::Rational && static_cast(&newExponent)->isOne()) { + Expression result = Power(this).childAtIndex(0); + return result; + } + return denominator; } + return Expression(); +} - /* 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(operand(1)); - // x^0 - if (b->isZero()) { - // 0^0 = undef - if (operand(0)->type() == Type::Rational && static_cast(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(RationalReference(1), true); - } - // x^1 - if (b->isOne()) { - return replaceWith(editableOperand(0), true); - } - } - if (operand(0)->type() == Type::Rational) { - Rational * a = static_cast(editableOperand(0)); - // 0^x - if (a->isZero()) { - if (operand(1)->sign() == Sign::Positive) { - return replaceWith(RationalReference(0), true); - } - if (operand(1)->sign() == Sign::Negative) { - return replaceWith(new Undefined(), true); - } - } - // 1^x - if (a->isOne()) { - return replaceWith(RationalReference(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(operand(1)); - // i^(p/q) - if (operand(0)->type() == Type::Symbol && static_cast(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(editableOperand(0)); - // p^q with p, q rationals - if (!letPowerAtRoot && operand(1)->type() == Type::Rational) { - Rational * exp = static_cast(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(operand(1))->numerator().isOne() && static_cast(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->numberOfChildren()-1); - static_cast(m)->removeOperand(i, false); - if (angleUnit == Preferences::AngleUnit::Degree) { - const Expression * pi = m->operand(m->numberOfChildren()-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)->numberOfChildren() == 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)->numberOfChildren() == 1 && operand(0)->type() == Type::Rational && static_cast(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(editableOperand(0)); - // Check is a > 0 or c is Integer - if (p->operand(0)->sign() == Sign::Positive || - (operand(1)->type() == Type::Rational && static_cast(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(editableOperand(0)); - // (a*b*c*...)^n = a^n*b^n*c^n*... if n integer - if (operand(1)->type() == Type::Rational && static_cast(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->numberOfChildren(); i++) { - // a is signed and a != -1 - if (m->operand(i)->sign() != Sign::Unknown && (m->operand(i)->type() != Type::Rational || !static_cast(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(operand(0))->isZero() && operand(1)->type() == Type::Addition) { - Addition * a = static_cast(editableOperand(1)); - // Check is b is rational - if (a->operand(0)->type() == Type::Rational) { - const Rational * rationalBase = static_cast(operand(0)); - const Rational * rationalIndex = static_cast(a->operand(0)); - if (RationalExponentShouldNotBeReduced(rationalBase, rationalIndex)) { - return this; - } - Power * p1 = static_cast(clone()); - replaceOperand(a, a->editableOperand(1), true); - Power * p2 = static_cast(clone()); - Multiplication * m = new Multiplication(p1, p2, false); - simplifyRationalRationalPower(p1, static_cast(p1->editableOperand(0)), static_cast(p1->editableOperand(1)->editableOperand(0)), 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(operand(1))->denominator().isOne() && operand(0)->type() == Type::Addition) { - // Exponent n - Rational * nr = static_cast(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)->numberOfChildren(); - /* 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(clippedN), static_cast(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->numberOfChildren(); 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(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); - } - } +//TODO #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(operand(1))->denominator().isOne() && operand(0)->type() == Type::Addition && operand(0)->numberOfChildren() == 2) { - Rational * nr = static_cast(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(BinomialCoefficient::compute(static_cast(i), static_cast(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()->numberOfChildren() == 1) { - if (operand(0)->type() == Type::Rational && static_cast(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(operand(0))->name() == Ion::Charset::Exponential) { - return true; - } - return false; -} - Expression * Power::simplifyPowerPower(Power * p, Expression * e, Context& context, Preferences::AngleUnit angleUnit) { Expression * p0 = p->editableOperand(0); Expression * p1 = p->editableOperand(1); @@ -554,6 +186,27 @@ Expression * Power::simplifyPowerPower(Power * p, Expression * e, Context& conte return shallowReduce(context, angleUnit); } +bool Power::parentIsALogarithmOfSameBase() const { + if (parent()->type() == Type::Logarithm && parent()->childAtIndex(0) == this) { + // parent = log(10^x) + if (parent()->numberOfChildren() == 1) { + if (childAtIndex(0)->type() == Type::Rational && static_cast(childAtIndex(0))->isTen()) { + return true; + } + return false; + } + // parent = log(x^y,x) + if (childAtIndex(0)->isIdenticalTo(parent()->childAtIndex(1))) { + return true; + } + } + // parent = ln(e^x) + if (parent()->type() == Type::NaperianLogarithm && parent()->childAtIndex(0) == this && childAtIndex(0)->type() == Type::Symbol && static_cast(childAtIndex(0))->name() == Ion::Charset::Exponential) { + return true; + } + return false; +} + Expression * Power::simplifyPowerMultiplication(Multiplication * m, Expression * r, Context& context, Preferences::AngleUnit angleUnit) { for (int index = 0; index < m->numberOfChildren(); index++) { Expression * factor = m->editableOperand(index); @@ -584,7 +237,10 @@ Expression * Power::simplifyRationalRationalPower(Expression * result, Rational result->replaceWith(m, true); return m->shallowReduce(context, angleUnit); } +#endif +//TODO +#if 0 Expression * Power::CreateSimplifiedIntegerRationalPower(Integer i, Rational * r, bool isDenominator, Context & context, Preferences::AngleUnit angleUnit) { assert(!i.isZero()); assert(r->sign() == Sign::Positive); @@ -635,7 +291,10 @@ Expression * Power::CreateSimplifiedIntegerRationalPower(Integer i, Rational * r m->sortOperands(SimplificationOrder, false); return m; } +#endif +//TODO +#if 0 Expression * Power::CreateNthRootOfUnity(const Rational r) { const Symbol * exp = new Symbol(Ion::Charset::Exponential); const Symbol * iComplex = new Symbol(Ion::Charset::IComplex); @@ -659,56 +318,14 @@ Expression * Power::CreateNthRootOfUnity(const Rational r) { #endif } -Expression Power::shallowBeautify(Context& context, Preferences::AngleUnit angleUnit) { - // X^-y -> 1/(X->shallowBeautify)^y - if (operand(1)->sign() == Sign::Negative) { - Expression * p = denominator(context, angleUnit); - Division * d = new Division(RationalReference(1), p, false); - p->shallowReduce(context, angleUnit); - replaceWith(d, true); - return d->shallowBeautify(context, angleUnit); - } - if (operand(1)->type() == Type::Rational && static_cast(operand(1))->numerator().isOne()) { - Integer index = static_cast(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::denominator(Context & context, Preferences::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(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(e->operand(1))->isHalf()) { + if (e->type() == Type::Power && e->childAtIndex(0)->type() == Type::Rational && e->childAtIndex(1)->type() == Type::Rational && static_cast(e->childAtIndex(1))->isHalf()) { return true; } - if (e->type() == Type::Multiplication && e->numberOfChildren() == 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(e->operand(1)->operand(1))->isHalf()) { + if (e->type() == Type::Multiplication && e->numberOfChildren() == 2 && e->childAtIndex(0)->type() == Type::Rational && e->childAtIndex(1)->type() == Type::Power && e->childAtIndex(1)->childAtIndex(0)->type() == Type::Rational && e->childAtIndex(1)->childAtIndex(1)->type() == Type::Rational && static_cast(e->childAtIndex(1)->childAtIndex(1))->isHalf()) { return true; } return false; @@ -719,13 +336,13 @@ const Rational * Power::RadicandInExpression(const Expression * e) { return nullptr; } else if (e->type() == Type::Power) { assert(e->type() == Type::Power); - assert(e->operand(0)->type() == Type::Rational); - return static_cast(e->operand(0)); + assert(e->childAtIndex(0)->type() == Type::Rational); + return static_cast(e->childAtIndex(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(e->operand(1)->operand(0)); + assert(e->childAtIndex(1)->type() == Type::Power); + assert(e->childAtIndex(1)->childAtIndex(0)->type() == Type::Rational); + return static_cast(e->childAtIndex(1)->childAtIndex(0)); } } @@ -736,32 +353,32 @@ const Rational * Power::RationalFactorInExpression(const Expression * e) { return nullptr; } else { assert(e->type() == Type::Multiplication); - assert(e->operand(0)->type() == Type::Rational); - return static_cast(e->operand(0)); + assert(e->childAtIndex(0)->type() == Type::Rational); + return static_cast(e->childAtIndex(0)); } } Expression * Power::removeSquareRootsFromDenominator(Context & context, Preferences::AngleUnit angleUnit) { Expression * result = nullptr; - if (operand(0)->type() == Type::Rational && operand(1)->type() == Type::Rational && (static_cast(operand(1))->isHalf() || static_cast(operand(1))->isMinusHalf())) { + if (childAtIndex(0)->type() == Type::Rational && childAtIndex(1)->type() == Type::Rational && (static_cast(childAtIndex(1))->isHalf() || static_cast(childAtIndex(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(operand(0))->numerator(); + Integer p = static_cast(childAtIndex(0))->numerator(); assert(!p.isZero()); // We eliminated case of form 0^(-1/2) at first step of shallowReduce - Integer q = static_cast(operand(0))->denominator(); + Integer q = static_cast(childAtIndex(0))->denominator(); // We do nothing for terms of the form sqrt(p) - if (!q.isOne() || static_cast(operand(1))->isMinusHalf()) { + if (!q.isOne() || static_cast(childAtIndex(1))->isMinusHalf()) { Power * sqrt = new Power(new Rational(Integer::Multiplication(p, q)), new Rational(1, 2), false); - if (static_cast(operand(1))->isHalf()) { + if (static_cast(childAtIndex(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(operand(1))->isMinusOne() && operand(0)->type() == Type::Addition && operand(0)->numberOfChildren() == 2 && TermIsARationalSquareRootOrRational(operand(0)->operand(0)) && TermIsARationalSquareRootOrRational(operand(0)->operand(1))) { + } else if (childAtIndex(1)->type() == Type::Rational && static_cast(childAtIndex(1))->isMinusOne() && childAtIndex(0)->type() == Type::Addition && childAtIndex(0)->numberOfChildren() == 2 && TermIsARationalSquareRootOrRational(childAtIndex(0)->childAtIndex(0)) && TermIsARationalSquareRootOrRational(childAtIndex(0)->childAtIndex(1))) { /* We're considering a term of the form * * 1/(n1/d1*sqrt(p1/q1) + n2/d2*sqrt(p2/q2)) @@ -772,10 +389,10 @@ Expression * Power::removeSquareRootsFromDenominator(Context & context, Preferen * ------------------------------------------------------- * 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)); + const Rational * f1 = RationalFactorInExpression(childAtIndex(0)->childAtIndex(0)); + const Rational * f2 = RationalFactorInExpression(childAtIndex(0)->childAtIndex(1)); + const Rational * r1 = RadicandInExpression(childAtIndex(0)->childAtIndex(0)); + const Rational * r2 = RadicandInExpression(childAtIndex(0)->childAtIndex(1)); Integer n1 = (f1 ? f1->numerator() : Integer(1)); Integer d1 = (f1 ? f1->denominator() : Integer(1)); Integer p1 = (r1 ? r1->numerator() : Integer(1)); @@ -829,27 +446,27 @@ Expression * Power::removeSquareRootsFromDenominator(Context & context, Preferen } bool Power::isNthRootOfUnity() const { - if (operand(0)->type() != Type::Symbol || static_cast(operand(0))->name() != Ion::Charset::Exponential) { + if (childAtIndex(0)->type() != Type::Symbol || static_cast(childAtIndex(0))->name() != Ion::Charset::Exponential) { return false; } - if (operand(1)->type() != Type::Multiplication) { + if (childAtIndex(1)->type() != Type::Multiplication) { return false; } - if (operand(1)->numberOfChildren() < 2 || operand(1)->numberOfChildren() > 3) { + if (childAtIndex(1)->numberOfChildren() < 2 || childAtIndex(1)->numberOfChildren() > 3) { return false; } - const Expression * i = operand(1)->operand(operand(1)->numberOfChildren()-1); + const Expression * i = childAtIndex(1)->childAtIndex(childAtIndex(1)->numberOfChildren()-1); if (i->type() != Type::Symbol || static_cast(i)->name() != Ion::Charset::IComplex) { return false; } - const Expression * pi = operand(1)->operand(operand(1)->numberOfChildren()-2); + const Expression * pi = childAtIndex(1)->childAtIndex(childAtIndex(1)->numberOfChildren()-2); if (pi->type() != Type::Symbol || static_cast(pi)->name() != Ion::Charset::SmallPi) { return false; } if (numberOfChildren() == 2) { return true; } - if (operand(1)->operand(0)->type() == Type::Rational) { + if (childAtIndex(1)->childAtIndex(0)->type() == Type::Rational) { return true; } return false; @@ -882,6 +499,444 @@ bool Power::RationalExponentShouldNotBeReduced(const Rational * b, const Rationa return false; } -template std::complex Power::compute(std::complex, std::complex); -template std::complex Power::compute(std::complex, std::complex); +#endif + + +template MatrixComplex PowerNode::computeOnComplexAndMatrix(const std::complex c, const MatrixComplex n) { + return MatrixComplex::Undefined(); +} + +template MatrixComplex PowerNode::computeOnMatrixAndComplex(const MatrixComplex m, const std::complex d) { + if (m.numberOfRows() != m.numberOfColumns()) { + return MatrixComplex::Undefined(); + } + T power = Complex(d).toScalar(); + if (std::isnan(power) || std::isinf(power) || power != (int)power || std::fabs(power) > k_maxApproximatePowerMatrix) { + return MatrixComplex::Undefined(); + } + if (power < 0) { + MatrixComplex inverse = m.inverse(); + if (!inverse.isDefined()) { + return MatrixComplex::Undefined(); + } + Complex minusC = Complex(-d); + MatrixComplex result = PowerNode::computeOnMatrixAndComplex(inverse, minusC); + return result; + } + MatrixComplex result = MatrixComplex::createIdentity(m.numberOfRows()); + // TODO: implement a quick exponentiation + for (int k = 0; k < (int)power; k++) { + if (Expression::shouldStopProcessing()) { + return MatrixComplex::Undefined(); + } + result = MultiplicationNode::computeOnMatrices(result, m); + } + return result; +} + +template MatrixComplex PowerNode::computeOnMatrices(const MatrixComplex m, const MatrixComplex n) { + return MatrixComplex::Undefined(); +} + +// Power + +Expression Power::setSign(ExpressionNode::Sign s, Context & context, Preferences::AngleUnit angleUnit) const { + assert(s == Sign::Positive); + assert(childAtIndex(0).sign() == Sign::Negative); + return Power(childAtIndex(0).setSign(ExpressionNode::Sign::Positive, context, angleUnit), childAtIndex(1)); +} + +int Power::getPolynomialCoefficients(char symbolName, Expression coefficients[]) const { + int deg = polynomialDegree(symbolName); + if (deg <= 0) { + return Expression::getPolynomialCoefficients(symbolName, coefficients); + } + /* Here we only consider the case x^4 as privateGetPolynomialCoefficients is + * supposed to be called after reducing the expression. */ + if (childAtIndex(0).type() == ExpressionNode::Type::Symbol && static_cast(childAtIndex(0).node())->name() == symbolName && childAtIndex(1).type() == ExpressionNode::Type::Rational) { + RationalNode * r = static_cast(childAtIndex(1).node()); + if (!r->denominator().isOne() || r->sign() == ExpressionNode::Sign::Negative) { + return -1; + } + NaturalIntegerPointer nip = r->numerator(); + Integer num = Integer(&nip); + if (Integer::NaturalOrder(num, Integer(Integer::k_maxExtractableInteger)) > 0) { + return -1; + } + int n = num.extractedInt(); + if (n <= k_maxPolynomialDegree) { + for (int i = 0; i < n; i++) { + coefficients[i] = Rational(0); + } + coefficients[n] = Rational(1); + return n; + } + } + return -1; +} + +Expression Power::shallowReduce(Context& context, Preferences::AngleUnit angleUnit) const { + Expression result = *this; + return result; + //TODO +#if 0 + Expression * e = Expression::shallowReduce(context, angleUnit); + if (e != this) { + return e; + } +#if MATRIX_EXACT_REDUCING + /* Step 0: get rid of matrix */ + if (childAtIndex(1)->type() == Type::Matrix) { + return replaceWith(new Undefined(), true); + } + if (childAtIndex(0)->type() == Type::Matrix) { + Matrix * mat = static_cast(editableOperand(0)); + if (childAtIndex(1)->type() != Type::Rational || !static_cast(childAtIndex(1))->denominator().isOne()) { + return replaceWith(new Undefined(), true); + } + Integer exponent = static_cast(childAtIndex(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 * op0 = static_cast *>(childAtIndex(0)->privateApproximate(float(), context, angleUnit)); + Complex * op1 = static_cast *>(childAtIndex(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 (childAtIndex(1)->type() == Type::Rational) { + const Rational * b = static_cast(childAtIndex(1)); + // x^0 + if (b->isZero()) { + // 0^0 = undef + if (childAtIndex(0)->type() == Type::Rational && static_cast(childAtIndex(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(RationalReference(1), true); + } + // x^1 + if (b->isOne()) { + return replaceWith(editableOperand(0), true); + } + } + if (childAtIndex(0)->type() == Type::Rational) { + Rational * a = static_cast(editableOperand(0)); + // 0^x + if (a->isZero()) { + if (childAtIndex(1)->sign() == Sign::Positive) { + return replaceWith(RationalReference(0), true); + } + if (childAtIndex(1)->sign() == Sign::Negative) { + return replaceWith(new Undefined(), true); + } + } + // 1^x + if (a->isOne()) { + return replaceWith(RationalReference(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 (childAtIndex(1)->type() == Type::Rational) { + const Rational * b = static_cast(childAtIndex(1)); + // i^(p/q) + if (childAtIndex(0)->type() == Type::Symbol && static_cast(childAtIndex(0))->name() == Ion::Charset::IComplex) { + Rational r = Rational::Multiplication(*b, Rational(1, 2)); + return replaceWith(CreateNthRootOfUnity(r))->shallowReduce(context, angleUnit); + } + } + bool letPowerAtRoot = parentIsALogarithmOfSameBase(); + if (childAtIndex(0)->type() == Type::Rational) { + Rational * a = static_cast(editableOperand(0)); + // p^q with p, q rationals + if (!letPowerAtRoot && childAtIndex(1)->type() == Type::Rational) { + Rational * exp = static_cast(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 && childAtIndex(1)->type() == Type::Rational && static_cast(childAtIndex(1))->numerator().isOne() && static_cast(childAtIndex(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->numberOfChildren()-1); + static_cast(m)->removeOperand(i, false); + if (angleUnit == Preferences::AngleUnit::Degree) { + const Expression * pi = m->childAtIndex(m->numberOfChildren()-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 (childAtIndex(1)->type() == Type::Logarithm) { + if (childAtIndex(1)->numberOfChildren() == 2 && childAtIndex(0)->isIdenticalTo(childAtIndex(1)->childAtIndex(1))) { + // y > 0 + if (childAtIndex(1)->childAtIndex(0)->sign() == Sign::Positive) { + return replaceWith(editableOperand(1)->editableOperand(0), true); + } + } + // 10^log(y) + if (childAtIndex(1)->numberOfChildren() == 1 && childAtIndex(0)->type() == Type::Rational && static_cast(childAtIndex(0))->isTen()) { + return replaceWith(editableOperand(1)->editableOperand(0), true); + } + } + // (a^b)^c -> a^(b*c) if a > 0 or c is integer + if (childAtIndex(0)->type() == Type::Power) { + Power * p = static_cast(editableOperand(0)); + // Check is a > 0 or c is Integer + if (p->childAtIndex(0)->sign() == Sign::Positive || + (childAtIndex(1)->type() == Type::Rational && static_cast(editableOperand(1))->denominator().isOne())) { + return simplifyPowerPower(p, editableOperand(1), context, angleUnit); + } + } + // (a*b*c*...)^r ? + if (!letPowerAtRoot && childAtIndex(0)->type() == Type::Multiplication) { + Multiplication * m = static_cast(editableOperand(0)); + // (a*b*c*...)^n = a^n*b^n*c^n*... if n integer + if (childAtIndex(1)->type() == Type::Rational && static_cast(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->numberOfChildren(); i++) { + // a is signed and a != -1 + if (m->childAtIndex(i)->sign() != Sign::Unknown && (m->childAtIndex(i)->type() != Type::Rational || !static_cast(m->childAtIndex(i))->isMinusOne())) { + //if (m->childAtIndex(i)->sign() == Sign::Positive || m->childAtIndex(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 && childAtIndex(0)->type() == Type::Rational && !static_cast(childAtIndex(0))->isZero() && childAtIndex(1)->type() == Type::Addition) { + Addition * a = static_cast(editableOperand(1)); + // Check is b is rational + if (a->childAtIndex(0)->type() == Type::Rational) { + const Rational * rationalBase = static_cast(childAtIndex(0)); + const Rational * rationalIndex = static_cast(a->childAtIndex(0)); + if (RationalExponentShouldNotBeReduced(rationalBase, rationalIndex)) { + return this; + } + Power * p1 = static_cast(clone()); + replaceOperand(a, a->editableOperand(1), true); + Power * p2 = static_cast(clone()); + Multiplication * m = new Multiplication(p1, p2, false); + simplifyRationalRationalPower(p1, static_cast(p1->editableOperand(0)), static_cast(p1->editableOperand(1)->editableOperand(0)), 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 && childAtIndex(1)->type() == Type::Rational && static_cast(childAtIndex(1))->denominator().isOne() && childAtIndex(0)->type() == Type::Addition) { + // Exponent n + Rational * nr = static_cast(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 = childAtIndex(0)->numberOfChildren(); + /* 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(clippedN), static_cast(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->numberOfChildren(); 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(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 && childAtIndex(1)->type() == Type::Rational && static_cast(childAtIndex(1))->denominator().isOne() && childAtIndex(0)->type() == Type::Addition && childAtIndex(0)->numberOfChildren() == 2) { + Rational * nr = static_cast(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(BinomialCoefficient::compute(static_cast(i), static_cast(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; +#endif +} + +Expression Power::shallowBeautify(Context& context, Preferences::AngleUnit angleUnit) const { + Expression result = *this; + return result; + //TODO +#if 0 + // X^-y -> 1/(X->shallowBeautify)^y + if (childAtIndex(1)->sign() == Sign::Negative) { + Expression * p = denominator(context, angleUnit); + Division * d = new Division(RationalReference(1), p, false); + p->shallowReduce(context, angleUnit); + replaceWith(d, true); + return d->shallowBeautify(context, angleUnit); + } + if (childAtIndex(1)->type() == Type::Rational && static_cast(childAtIndex(1))->numerator().isOne()) { + Integer index = static_cast(childAtIndex(1))->denominator(); + if (index.isEqualTo(Integer(2))) { + const Expression * sqrtOperand[1] = {childAtIndex(0)}; + SquareRoot * sqr = new SquareRoot(sqrtOperand, true); + return replaceWith(sqr, true); + } + const Expression * rootOperand[2] = {childAtIndex(0)->clone(), new Rational(index)}; + NthRoot * nr = new NthRoot(rootOperand, false); + return replaceWith(nr, true); + } + // +(a,b)^c ->(+(a,b))^c + if (childAtIndex(0)->type() == Type::Addition || childAtIndex(0)->type() == Type::Multiplication) { + const Expression * o[1] = {childAtIndex(0)}; + Parenthesis * p = new Parenthesis(o, true); + replaceOperand(childAtIndex(0), p, true); + } + return this; +#endif +} + +template Complex PowerNode::compute(std::complex, std::complex); +template Complex PowerNode::compute(std::complex, std::complex); }