From 12a5f5499ccf56e180578d19ee3207cb27d5183a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89milie=20Feral?= Date: Wed, 26 Feb 2020 11:06:47 +0100 Subject: [PATCH] [poincare] Handle real root finding in Power and NthRoot in order to remove beautifying x^(p/q) -> root(x,q)^p. This triggered precision loss! --- poincare/include/poincare/power.h | 6 ++- poincare/src/nth_root.cpp | 16 +++---- poincare/src/power.cpp | 69 ++++++++++++++++++++++++------- poincare/test/approximation.cpp | 7 ++-- poincare/test/simplification.cpp | 4 +- 5 files changed, 70 insertions(+), 32 deletions(-) diff --git a/poincare/include/poincare/power.h b/poincare/include/poincare/power.h index c5860e77e..b2e161b85 100644 --- a/poincare/include/poincare/power.h +++ b/poincare/include/poincare/power.h @@ -34,6 +34,7 @@ public: int polynomialDegree(Context * context, const char * symbolName) const override; int getPolynomialCoefficients(Context * context, const char * symbolName, Expression coefficients[], ExpressionNode::SymbolicComputation symbolicComputation) const override; + template static Complex computeRealRootOfRationalPow(const std::complex c, T p, T q); template static Complex compute(const std::complex c, const std::complex d, Preferences::ComplexFormat complexFormat); private: @@ -59,11 +60,12 @@ private: template static MatrixComplex computeOnMatrixAndComplex(const MatrixComplex m, const std::complex d, Preferences::ComplexFormat complexFormat); template static MatrixComplex computeOnMatrices(const MatrixComplex m, const MatrixComplex n, Preferences::ComplexFormat complexFormat); Evaluation approximate(SinglePrecision p, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const override { - return ApproximationHelper::MapReduce(this, context, complexFormat, angleUnit, compute, computeOnComplexAndMatrix, computeOnMatrixAndComplex, computeOnMatrices); + return templatedApproximate(context, complexFormat, angleUnit); } Evaluation approximate(DoublePrecision p, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const override { - return ApproximationHelper::MapReduce(this, context, complexFormat, angleUnit, compute, computeOnComplexAndMatrix, computeOnMatrixAndComplex, computeOnMatrices); + return templatedApproximate(context, complexFormat, angleUnit); } + template Evaluation templatedApproximate(Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const; }; class Power final : public Expression { diff --git a/poincare/src/nth_root.cpp b/poincare/src/nth_root.cpp index 74ee32f52..b5081c1c0 100644 --- a/poincare/src/nth_root.cpp +++ b/poincare/src/nth_root.cpp @@ -46,18 +46,12 @@ Evaluation NthRootNode::templatedApproximate(Context * context, Preferences:: /* If the complexFormat is Real, we look for nthroot of form root(x,q) with * x real and q integer because they might have a real form which does not * correspond to the principale angle. */ - if (complexFormat == Preferences::ComplexFormat::Real) { + if (complexFormat == Preferences::ComplexFormat::Real && indexc.imag() == 0.0 && std::round(indexc.real()) == indexc.real()) { // root(x, q) with q integer and x real - if (basec.imag() == 0.0 && indexc.imag() == 0.0 && std::round(indexc.real()) == indexc.real()) { - std::complex absBasec = basec; - absBasec.real(std::fabs(absBasec.real())); - // compute root(|x|, q) - Complex absBasePowIndex = PowerNode::compute(absBasec, std::complex(1.0)/(indexc), complexFormat); - // q odd if (-1)^q = -1 - if (std::pow((T)-1.0, (T)indexc.real()) < 0.0) { - return basec.real() < 0 ? Complex::Builder(-absBasePowIndex.stdComplex()) : absBasePowIndex; - } - } + Complex result = PowerNode::computeRealRootOfRationalPow(basec, (T)1.0, indexc.real()); + if (!result.isUndefined()) { + return result; + } } result = PowerNode::compute(basec, std::complex(1.0)/(indexc), complexFormat); } diff --git a/poincare/src/power.cpp b/poincare/src/power.cpp index 798de23ae..29d25e7c4 100644 --- a/poincare/src/power.cpp +++ b/poincare/src/power.cpp @@ -127,6 +127,22 @@ bool PowerNode::childAtIndexNeedsUserParentheses(const Expression & child, int c // Private +template +Complex PowerNode::computeRealRootOfRationalPow(const std::complex c, T p, T q) { + // Compute real root of c^(p/q) with p, q integers if it has a real root + if (c.imag() == 0 && std::pow((T)-1.0, q) < 0.0) { + /* If c real and q odd integer (q odd if (-1)^q = -1), a real root does + * exist (which is not necessarily the principal root)! */ + std::complex absc = c; + absc.real(std::fabs(absc.real())); + // compute |c|^(p/q) which is a real + Complex absCPowD = PowerNode::compute(absc, std::complex(p/q), Preferences::ComplexFormat::Real); + // c^(p/q) = (sign(c)^p)*|c|^(p/q) = -|c|^(p/q) iff c < 0 and p odd + return c.real() < 0 && std::pow((T)-1.0, p) < 0.0 ? Complex::Builder(-absCPowD.stdComplex()) : absCPowD; + } + return Complex::Undefined(); +} + template Complex PowerNode::compute(const std::complex c, const std::complex d, Preferences::ComplexFormat complexFormat) { std::complex result; @@ -262,6 +278,45 @@ template MatrixComplex PowerNode::computeOnMatrices(const MatrixC return MatrixComplex::Undefined(); } +template Evaluation PowerNode::templatedApproximate(Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const { + /* Special case: c^(p/q) with p, q integers + * In real mode, c^(p/q) might have a real root which is not the principal + * root. We return this value in that case to avoid returning "unreal". */ + if (complexFormat == Preferences::ComplexFormat::Real) { + Evaluation base = childAtIndex(0)->approximate(T(), context, complexFormat, angleUnit); + if (base.type() != EvaluationNode::Type::Complex) { + goto defaultApproximation; + } + std::complex c = static_cast &>(base).stdComplex(); + T p = NAN; + T q = NAN; + if (childAtIndex(1)->type() == ExpressionNode::Type::Rational) { + const RationalNode * r = static_cast(childAtIndex(1)); + p = r->signedNumerator().approximate(); + q = r->denominator().approximate(); + } + // Spot form p/q with p, q integers + if (childAtIndex(1)->type() == ExpressionNode::Type::Division && childAtIndex(1)->childAtIndex(0)->type() == ExpressionNode::Type::Rational && childAtIndex(1)->childAtIndex(1)->type() == ExpressionNode::Type::Rational) { + const RationalNode * pRat = static_cast(childAtIndex(1)->childAtIndex(0)); + const RationalNode * qRat = static_cast(childAtIndex(1)->childAtIndex(1)); + if (!pRat->denominator().isOne() || !qRat->denominator().isOne()) { + goto defaultApproximation; + } + p = pRat->signedNumerator().approximate(); + q = qRat->signedNumerator().approximate(); + } + if (std::isnan(p) || std::isnan(q)) { + goto defaultApproximation; + } + Complex result = computeRealRootOfRationalPow(c, p, q); + if (!result.isUndefined()) { + return result; + } + } +defaultApproximation: + return ApproximationHelper::MapReduce(this, context, complexFormat, angleUnit, compute, computeOnComplexAndMatrix, computeOnMatrixAndComplex, computeOnMatrices); +} + // Power Expression Power::setSign(ExpressionNode::Sign s, ExpressionNode::ReductionContext reductionContext) { @@ -901,20 +956,6 @@ Expression Power::shallowBeautify(ExpressionNode::ReductionContext reductionCont return result; } - /* Optional Step 3: if the ReductionTarget is the SystemForApproximation or - * SystemForAnalysis, turn a^(p/q) into (root(a, q))^p - * Indeed, root(a, q) can have a real root which is not the principale angle - * but that we want to return in real complex format. This special case is - * handled in NthRoot approximation but not in Power approximation. */ - if (reductionContext.target() != ExpressionNode::ReductionTarget::User && childAtIndex(1).type() == ExpressionNode::Type::Rational) { - Integer p = childAtIndex(1).convert().signedIntegerNumerator(); - Integer q = childAtIndex(1).convert().integerDenominator(); - Expression nthRoot = q.isOne() ? childAtIndex(0) : NthRoot::Builder(childAtIndex(0), Rational::Builder(q)); - Expression result = p.isOne() ? nthRoot : Power::Builder(nthRoot, Rational::Builder(p)); - replaceWithInPlace(result); - return result; - } - // Step 4: Force Float(1) in front of an orphan Power of Unit if (parent().isUninitialized() && childAtIndex(0).type() == ExpressionNode::Type::Unit) { Multiplication m = Multiplication::Builder(Float::Builder(1.0)); diff --git a/poincare/test/approximation.cpp b/poincare/test/approximation.cpp index ee6c1580a..362cd621f 100644 --- a/poincare/test/approximation.cpp +++ b/poincare/test/approximation.cpp @@ -837,11 +837,12 @@ QUIZ_CASE(poincare_approximation_complex_format) { assert_expression_approximates_to("√(-1)", "unreal", Radian, Real); assert_expression_approximates_to("√(-1)×√(-1)", "unreal", Radian, Real); assert_expression_approximates_to("ln(-2)", "unreal", Radian, Real); - assert_expression_approximates_to("(-8)^(1/3)", "unreal", Radian, Real); // Power always approximates to the principal root (even if unreal) + // Power/Root approximates to the first REAL root in Real mode + assert_expression_simplifies_approximates_to("(-8)^(1/3)", "-2", Radian, Real); // Power have to be simplified first in order to spot the right form c^(p/q) with p, q integers assert_expression_approximates_to("root(-8,3)", "-2", Radian, Real); // Root approximates to the first REAL root in Real mode assert_expression_approximates_to("8^(1/3)", "2", Radian, Real); - assert_expression_approximates_to("(-8)^(2/3)", "unreal", Radian, Real); // Power always approximates to the principal root (even if unreal) - assert_expression_approximates_to("root(-8, 3)^2", "4", Radian, Real); // Root approximates to the first REAL root in Real mode + assert_expression_simplifies_approximates_to("(-8)^(2/3)", "4", Radian, Real); // Power have to be simplified first (cf previous comment) + assert_expression_approximates_to("root(-8, 3)^2", "4", Radian, Real); assert_expression_approximates_to("root(-8,3)", "-2", Radian, Real); // Cartesian diff --git a/poincare/test/simplification.cpp b/poincare/test/simplification.cpp index 0dd49cc10..cd5bd1c04 100644 --- a/poincare/test/simplification.cpp +++ b/poincare/test/simplification.cpp @@ -1231,8 +1231,8 @@ QUIZ_CASE(poincare_simplification_reduction_target) { assert_parsed_expression_simplify_to("(1+x)/(1+x)", "1", User); // Apply rule x^(2/3) --> root(x,3)^2 for ReductionTarget = System - assert_parsed_expression_simplify_to("x^(2/3)", "root(x,3)^2", SystemForApproximation); - assert_parsed_expression_simplify_to("x^(2/3)", "root(x,3)^2", SystemForAnalysis); + assert_parsed_expression_simplify_to("x^(2/3)", "x^\u00122/3\u0013", SystemForApproximation); + assert_parsed_expression_simplify_to("x^(2/3)", "x^\u00122/3\u0013", SystemForAnalysis); assert_parsed_expression_simplify_to("x^(2/3)", "x^\u00122/3\u0013", User); assert_parsed_expression_simplify_to("x^(1/3)", "root(x,3)", SystemForApproximation); assert_parsed_expression_simplify_to("x^(1/3)", "root(x,3)", SystemForAnalysis);