[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

@@ -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<typename T> static Complex<T> computeRealRootOfRationalPow(const std::complex<T> c, T p, T q);
   template<typename T> static Complex<T> compute(const std::complex<T> c, const std::complex<T> d, Preferences::ComplexFormat complexFormat);
 
 private:
@@ -59,11 +60,12 @@ private:
   template<typename T> static MatrixComplex<T> computeOnMatrixAndComplex(const MatrixComplex<T> m, const std::complex<T> d, Preferences::ComplexFormat complexFormat);
   template<typename T> static MatrixComplex<T> computeOnMatrices(const MatrixComplex<T> m, const MatrixComplex<T> n, Preferences::ComplexFormat complexFormat);
   Evaluation<float> approximate(SinglePrecision p, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const override {
-    return ApproximationHelper::MapReduce<float>(this, context, complexFormat, angleUnit, compute<float>, computeOnComplexAndMatrix<float>, computeOnMatrixAndComplex<float>, computeOnMatrices<float>);
+    return templatedApproximate<float>(context, complexFormat, angleUnit);
   }
   Evaluation<double> approximate(DoublePrecision p, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const override {
-    return ApproximationHelper::MapReduce<double>(this, context, complexFormat, angleUnit, compute<double>, computeOnComplexAndMatrix<double>, computeOnMatrixAndComplex<double>, computeOnMatrices<double>);
+    return templatedApproximate<double>(context, complexFormat, angleUnit);
   }
+ template<typename T> Evaluation<T> templatedApproximate(Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const;
 };
 
 class Power final : public Expression {

poincare/src/nth_root.cpp

@@ -46,18 +46,12 @@ Evaluation<T> 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<T> absBasec = basec;
-        absBasec.real(std::fabs(absBasec.real()));
-        // compute root(|x|, q)
-        Complex<T> absBasePowIndex = PowerNode::compute(absBasec, std::complex<T>(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<T>::Builder(-absBasePowIndex.stdComplex()) : absBasePowIndex;
-        }
-      }
+      Complex<T> result = PowerNode::computeRealRootOfRationalPow(basec, (T)1.0, indexc.real());
+       if (!result.isUndefined()) {
+         return result;
+       }
     }
     result = PowerNode::compute(basec, std::complex<T>(1.0)/(indexc), complexFormat);
   }

poincare/src/power.cpp

@@ -127,6 +127,22 @@ bool PowerNode::childAtIndexNeedsUserParentheses(const Expression & child, int c
 
 // Private
 
+template<typename T>
+Complex<T> PowerNode::computeRealRootOfRationalPow(const std::complex<T> 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<T> absc = c;
+    absc.real(std::fabs(absc.real()));
+    // compute |c|^(p/q) which is a real
+    Complex<T> absCPowD = PowerNode::compute(absc, std::complex<T>(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<T>::Builder(-absCPowD.stdComplex()) : absCPowD;
+  }
+  return Complex<T>::Undefined();
+}
+
 template<typename T>
 Complex<T> PowerNode::compute(const std::complex<T> c, const std::complex<T> d, Preferences::ComplexFormat complexFormat) {
   std::complex<T> result;
@@ -262,6 +278,45 @@ template<typename T> MatrixComplex<T> PowerNode::computeOnMatrices(const MatrixC
   return MatrixComplex<T>::Undefined();
 }
 
+template<typename T> Evaluation<T> 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<T> base = childAtIndex(0)->approximate(T(), context, complexFormat, angleUnit);
+    if (base.type() != EvaluationNode<T>::Type::Complex) {
+      goto defaultApproximation;
+    }
+    std::complex<T> c = static_cast<Complex<T> &>(base).stdComplex();
+    T p = NAN;
+    T q = NAN;
+    if (childAtIndex(1)->type() == ExpressionNode::Type::Rational) {
+      const RationalNode * r = static_cast<const RationalNode *>(childAtIndex(1));
+      p = r->signedNumerator().approximate<T>();
+      q = r->denominator().approximate<T>();
+    }
+    // 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<const RationalNode *>(childAtIndex(1)->childAtIndex(0));
+      const RationalNode * qRat = static_cast<const RationalNode *>(childAtIndex(1)->childAtIndex(1));
+      if (!pRat->denominator().isOne() || !qRat->denominator().isOne()) {
+        goto defaultApproximation;
+      }
+      p = pRat->signedNumerator().approximate<T>();
+      q = qRat->signedNumerator().approximate<T>();
+    }
+    if (std::isnan(p) || std::isnan(q)) {
+      goto defaultApproximation;
+    }
+    Complex<T> result = computeRealRootOfRationalPow(c, p, q);
+    if (!result.isUndefined()) {
+      return result;
+    }
+  }
+defaultApproximation:
+  return ApproximationHelper::MapReduce<T>(this, context, complexFormat, angleUnit, compute<T>, computeOnComplexAndMatrix<T>, computeOnMatrixAndComplex<T>, computeOnMatrices<T>);
+}
+
 // 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<Rational>().signedIntegerNumerator();
-    Integer q = childAtIndex(1).convert<Rational>().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<double>::Builder(1.0));

poincare/test/approximation.cpp

@@ -837,11 +837,12 @@ QUIZ_CASE(poincare_approximation_complex_format) {
   assert_expression_approximates_to<double>("√(-1)", "unreal", Radian, Real);
   assert_expression_approximates_to<double>("√(-1)×√(-1)", "unreal", Radian, Real);
   assert_expression_approximates_to<double>("ln(-2)", "unreal", Radian, Real);
-  assert_expression_approximates_to<double>("(-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<double>("(-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<double>("root(-8,3)", "-2", Radian, Real); // Root approximates to the first REAL root in Real mode
   assert_expression_approximates_to<double>("8^(1/3)", "2", Radian, Real);
-  assert_expression_approximates_to<float>("(-8)^(2/3)", "unreal", Radian, Real); // Power always approximates to the principal root (even if unreal)
-  assert_expression_approximates_to<float>("root(-8, 3)^2", "4", Radian, Real); // Root approximates to the first REAL root in Real mode
+  assert_expression_simplifies_approximates_to<float>("(-8)^(2/3)", "4", Radian, Real); // Power have to be simplified first (cf previous comment)
+  assert_expression_approximates_to<float>("root(-8, 3)^2", "4", Radian, Real);
   assert_expression_approximates_to<double>("root(-8,3)", "-2", Radian, Real);
 
   // Cartesian

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);