From c4018d06486de4cca3dd610a5a3bc878603ffbcf Mon Sep 17 00:00:00 2001 From: Hugo Saint-Vignes Date: Thu, 2 Jul 2020 14:47:54 +0200 Subject: [PATCH] [poincare] Update rowCanonize pivot selection for more consistent ref results Change-Id: Id7e856f57ccd3d990077b0f6753089bc6edcc03b --- poincare/src/matrix.cpp | 55 ++++++++++++++++++++++++++------ poincare/test/approximation.cpp | 23 +++++++------ poincare/test/simplification.cpp | 1 + 3 files changed, 61 insertions(+), 18 deletions(-) diff --git a/poincare/src/matrix.cpp b/poincare/src/matrix.cpp index 0c5f37cd7..cabbb2992 100644 --- a/poincare/src/matrix.cpp +++ b/poincare/src/matrix.cpp @@ -1,4 +1,5 @@ #include +#include #include #include #include @@ -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(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()) { - 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 *, int, int); template int Matrix::ArrayInverse(double *, int, int); template int Matrix::ArrayInverse>(std::complex *, int, int); template int Matrix::ArrayInverse>(std::complex *, int, int); diff --git a/poincare/test/approximation.cpp b/poincare/test/approximation.cpp index 58f2304a0..2e294d1f6 100644 --- a/poincare/test/approximation.cpp +++ b/poincare/test/approximation.cpp @@ -191,6 +191,9 @@ QUIZ_CASE(poincare_approximation_division) { assert_expression_approximates_to("[[1,2][3,4]]/[[3,4][6,9]]", "[[-1,6.6666666666667ᴇ-1][1,0]]"); assert_expression_approximates_to("3/[[3,4][5,6]]", "[[-9,6][7.5,-4.5]]"); assert_expression_approximates_to("(3+4𝐢)/[[1,𝐢][3,4]]", "[[4×𝐢,1][-3×𝐢,𝐢]]"); + // TODO: get rid of the neglectable real or imaginary parts + assert_expression_approximates_to("(3+4𝐢)/[[3,4][1,𝐢]]", "[[1+5.5511151231258ᴇ-17×𝐢,-2.2204460492503ᴇ-16+4×𝐢][𝐢,-3×𝐢]]"); + // [[1,4×𝐢][𝐢,-3×𝐢]] is expected assert_expression_approximates_to("1ᴇ20/(1ᴇ20+1ᴇ20𝐢)", "0.5-0.5×𝐢"); assert_expression_approximates_to("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("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("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("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("ref([[1,0][5,6][0,10]])", "[[1,0][0,1][0,0]]"); - // --> "[[1,1.2][0,1][0,0]]" + assert_expression_approximates_to("ref([[1,0][5,6][0,10]])", "[[1,1.2][0,1][0,0]]"); assert_expression_approximates_to("rref([[1,0][5,6][0,10]])", "[[1,0][0,1][0,0]]"); assert_expression_approximates_to("ref([[0,0][0,0][0,0]])", "[[0,0][0,0][0,0]]"); assert_expression_approximates_to("rref([[0,0][0,0][0,0]])", "[[0,0][0,0][0,0]]"); - assert_expression_approximates_to("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("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("rref([[0,2,-1][5,6,7][12,11,10]])", "[[1,0,0][0,1,0][0,0,1]]"); + assert_expression_approximates_to("ref([[3,9][2,5]])", "[[1,3][0,1]]"); + assert_expression_approximates_to("ref([[3,2][5,7]])", "[[1,1.4][0,1]]"); + assert_expression_approximates_to("ref([[3,11][5,7]])", "[[1,1.4][0,1]]"); + assert_expression_approximates_to("ref([[2,5][2,7]])", "[[1,2.5][0,1]]"); + assert_expression_approximates_to("ref([[3,12][-4,1]])", "[[1,-0.25][0,1]]"); + assert_expression_approximates_to("ref([[0,1][1ᴇ-100,1]])", "[[1,1ᴇ100][0,1]]"); + assert_expression_approximates_to("rref([[0,1][1ᴇ-100,1]])", "[[1,0][0,1]]"); assert_expression_approximates_to("round(2.3246,3)", "2.325"); assert_expression_approximates_to("round(2.3245,3)", "2.325"); diff --git a/poincare/test/simplification.cpp b/poincare/test/simplification.cpp index f9bf1c3ab..ebf104323 100644 --- a/poincare/test/simplification.cpp +++ b/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");