[poincare] Update rowCanonize pivot selection for more consistent ref results

Change-Id: Id7e856f57ccd3d990077b0f6753089bc6edcc03b

poincare/src/matrix.cpp

@@ -1,4 +1,5 @@
 #include <poincare/matrix.h>
+#include <poincare/absolute_value.h>
 #include <poincare/addition.h>
 #include <poincare/division.h>
 #include <poincare/exception_checkpoint.h>
@@ -204,12 +205,35 @@ Matrix Matrix::rowCanonize(ExpressionNode::ReductionContext reductionContext, Ex
   int k = 0; // column pivot
 
   while (h < m && k < n) {
-    // Find the first non-null pivot
+    /* In non-reduced form, the pivot selection method will affect the output.
+     * Here we prioritize the biggest pivot (in value) to get an output that
+     * does not depends on the order of the rows of the matrix.
+     * We could also take lowest non null pivots, or just first non null as we
+     * already do with reduced forms. Output would be different, but correct. */
+    int iPivot_temp = h;
     int iPivot = h;
-    while (iPivot < m && matrixChild(iPivot, k).isRationalZero()) {
-      iPivot++;
+    float bestPivot = 0.0;
+    while (iPivot_temp < m) {
+      // Using float to find the biggest pivot is sufficient.
+      float pivot = AbsoluteValue::Builder(matrixChild(iPivot_temp, k).clone()).approximateToScalar<float>(reductionContext.context(), reductionContext.complexFormat(), reductionContext.angleUnit());
+      // Handle very low pivots
+      if (pivot == 0.0f && !matrixChild(iPivot_temp, k).isRationalZero()) {
+        pivot = FLT_MIN;
+      }
+
+      if (pivot > bestPivot) {
+        // Update best pivot
+        bestPivot = pivot;
+        iPivot = iPivot_temp;
+        if (reduced) {
+          /* In reduced form, taking the first non null pivot is enough, and
+           * more efficient. */
+          break;
+        }
+      }
+      iPivot_temp++;
     }
-    if (iPivot == m) {
+    if (matrixChild(iPivot, k).isRationalZero()) {
       // No non-null coefficient in this column, skip
       k++;
       if (determinant) {
@@ -273,12 +297,26 @@ void Matrix::ArrayRowCanonize(T * array, int numberOfRows, int numberOfColumns,
   int k = 0; // column pivot
 
   while (h < numberOfRows && k < numberOfColumns) {
-    // Find the first non-null pivot
+    // Find the biggest pivot (in absolute value). See comment on rowCanonize.
+    int iPivot_temp = h;
     int iPivot = h;
-    while (iPivot < numberOfRows && std::abs(array[iPivot*numberOfColumns+k]) < Expression::Epsilon<double>()) {
-      iPivot++;
+    // Using double to stay accurate with any type T
+    double bestPivot = 0.0;
+    while (iPivot_temp < numberOfRows) {
+      double pivot = std::abs(array[iPivot_temp*numberOfColumns+k]);
+      if (pivot > bestPivot) {
+        // Update best pivot
+        bestPivot = pivot;
+        iPivot = iPivot_temp;
+        if (reduced) {
+          /* In reduced form, taking the first non null pivot is enough, and
+           * more efficient. */
+          break;
+        }
+      }
+      iPivot_temp++;
     }
-    if (iPivot == numberOfRows) {
+    if (bestPivot < DBL_MIN) {
       // No non-null coefficient in this column, skip
       k++;
       // Update determinant: det *= 0
@@ -522,7 +560,6 @@ Expression Matrix::computeInverseOrDeterminant(bool computeDeterminant, Expressi
 }
 
 
-template int Matrix::ArrayInverse<float>(float *, int, int);
 template int Matrix::ArrayInverse<double>(double *, int, int);
 template int Matrix::ArrayInverse<std::complex<float>>(std::complex<float> *, int, int);
 template int Matrix::ArrayInverse<std::complex<double>>(std::complex<double> *, int, int);

poincare/test/approximation.cpp

@@ -191,6 +191,9 @@ QUIZ_CASE(poincare_approximation_division) {
   assert_expression_approximates_to<double>("[[1,2][3,4]]/[[3,4][6,9]]", "[[-1,6.6666666666667ᴇ-1][1,0]]");
   assert_expression_approximates_to<double>("3/[[3,4][5,6]]", "[[-9,6][7.5,-4.5]]");
   assert_expression_approximates_to<double>("(3+4𝐢)/[[1,𝐢][3,4]]", "[[4×𝐢,1][-3×𝐢,𝐢]]");
+  // TODO: get rid of the neglectable real or imaginary parts
+  assert_expression_approximates_to<double>("(3+4𝐢)/[[3,4][1,𝐢]]", "[[1+5.5511151231258ᴇ-17×𝐢,-2.2204460492503ᴇ-16+4×𝐢][𝐢,-3×𝐢]]");
+  // [[1,4×𝐢][𝐢,-3×𝐢]] is expected
   assert_expression_approximates_to<float>("1ᴇ20/(1ᴇ20+1ᴇ20𝐢)", "0.5-0.5×𝐢");
   assert_expression_approximates_to<double>("1ᴇ155/(1ᴇ155+1ᴇ155𝐢)", "0.5-0.5×𝐢");
 
@@ -409,20 +412,22 @@ QUIZ_CASE(poincare_approximation_function) {
    * - Rows with only zeros must be at the bottom.
    * - Leading coefficient of other rows must be to the right (strictly) of the
    * - one above.
-   * - (Optional, but sometimes recommended) Leading coefficients must be 1.
-   * NOTE : It would be better if results for ref matched the one commented
-   *  bellow. */
-  assert_expression_approximates_to<double>("ref([[1,0,3,4][5,7,6,8][0,10,11,12]])", "[[1,0,3,4][0,1,-1.2857142857143,-1.7142857142857][0,0,1,1.2215568862275]]");
-  // --> "[[1,1.4,1.2,1.6][0,1,1.1,1.2][0,0,1,1.221557]]"
+   * - (Optional, but sometimes recommended) Leading coefficients must be 1. */
+  assert_expression_approximates_to<double>("ref([[1,0,3,4][5,7,6,8][0,10,11,12]])", "[[1,1.4,1.2,1.6][0,1,1.1,1.2][0,0,1,1.2215568862275]]");
   assert_expression_approximates_to<double>("rref([[1,0,3,4][5,7,6,8][0,10,11,12]])", "[[1,0,0,3.3532934131737ᴇ-1][0,1,0,-0.1437125748503][0,0,1,1.2215568862275]]");
-  assert_expression_approximates_to<double>("ref([[1,0][5,6][0,10]])", "[[1,0][0,1][0,0]]");
-  // --> "[[1,1.2][0,1][0,0]]"
+  assert_expression_approximates_to<double>("ref([[1,0][5,6][0,10]])", "[[1,1.2][0,1][0,0]]");
   assert_expression_approximates_to<double>("rref([[1,0][5,6][0,10]])", "[[1,0][0,1][0,0]]");
   assert_expression_approximates_to<double>("ref([[0,0][0,0][0,0]])", "[[0,0][0,0][0,0]]");
   assert_expression_approximates_to<double>("rref([[0,0][0,0][0,0]])", "[[0,0][0,0][0,0]]");
-  assert_expression_approximates_to<double>("ref([[0,2,-1][5,6,7][12,11,10]])", "[[1,1.2,1.4][0,1,-0.5][0,0,1]]");
-  // --> "[[1,0.9166667,0.8333333][0,1,-0.5][0,0,1]]"
+  assert_expression_approximates_to<double>("ref([[0,2,-1][5,6,7][12,11,10]])", "[[1,9.1666666666667ᴇ-1,8.3333333333333ᴇ-1][0,1,-0.5][0,0,1]]");
   assert_expression_approximates_to<double>("rref([[0,2,-1][5,6,7][12,11,10]])", "[[1,0,0][0,1,0][0,0,1]]");
+  assert_expression_approximates_to<double>("ref([[3,9][2,5]])", "[[1,3][0,1]]");
+  assert_expression_approximates_to<double>("ref([[3,2][5,7]])", "[[1,1.4][0,1]]");
+  assert_expression_approximates_to<double>("ref([[3,11][5,7]])", "[[1,1.4][0,1]]");
+  assert_expression_approximates_to<double>("ref([[2,5][2,7]])", "[[1,2.5][0,1]]");
+  assert_expression_approximates_to<double>("ref([[3,12][-4,1]])", "[[1,-0.25][0,1]]");
+  assert_expression_approximates_to<double>("ref([[0,1][1ᴇ-100,1]])", "[[1,1ᴇ100][0,1]]");
+  assert_expression_approximates_to<double>("rref([[0,1][1ᴇ-100,1]])", "[[1,0][0,1]]");
 
   assert_expression_approximates_to<float>("round(2.3246,3)", "2.325");
   assert_expression_approximates_to<double>("round(2.3245,3)", "2.325");

poincare/test/simplification.cpp

@@ -934,6 +934,7 @@ QUIZ_CASE(poincare_simplification_matrix) {
   assert_parsed_expression_simplify_to("rref([[1,1/√(2),√(4)]])", "[[1,√(2)/2,2]]");
   assert_parsed_expression_simplify_to("ref([[1,0,√(4)][0,1,1/√(2)][0,0,1]])", "[[1,0,2][0,1,√(2)/2][0,0,1]]");
   assert_parsed_expression_simplify_to("rref([[1,0,√(4)][0,1,1/√(2)][0,0,0]])", "[[1,0,2][0,1,√(2)/2][0,0,0]]");
+  assert_parsed_expression_simplify_to("ref([[1,0,3,4][5,7,6,8][0,10,11,12]])", "[[1,7/5,6/5,8/5][0,1,11/10,6/5][0,0,1,204/167]]");
 
   // Expressions with unreduced matrix
   assert_reduce("confidence(cos(2)/25,3)→a");