[poincare] Add ref and rref matrix functions

Change-Id: Id0e57a4f85d551ca5320c4b6e3c0baadae70946d

apps/math_toolbox.cpp

@@ -71,7 +71,9 @@ const ToolboxMessageTree matricesChildren[] = {
   ToolboxMessageTree::Leaf(I18n::Message::DeterminantCommandWithArg, I18n::Message::Determinant),
   ToolboxMessageTree::Leaf(I18n::Message::TransposeCommandWithArg, I18n::Message::Transpose),
   ToolboxMessageTree::Leaf(I18n::Message::TraceCommandWithArg, I18n::Message::Trace),
-  ToolboxMessageTree::Leaf(I18n::Message::DimensionCommandWithArg, I18n::Message::Dimension)
+  ToolboxMessageTree::Leaf(I18n::Message::DimensionCommandWithArg, I18n::Message::Dimension),
+  ToolboxMessageTree::Leaf(I18n::Message::RowEchelonFormCommandWithArg, I18n::Message::RowEchelonForm),
+  ToolboxMessageTree::Leaf(I18n::Message::ReducedRowEchelonFormCommandWithArg, I18n::Message::ReducedRowEchelonForm)
 };
 
 #if LIST_ARE_DEFINED

apps/shared.universal.i18n

@@ -140,9 +140,11 @@ QuoCommandWithArg = "quo(p,q)"
 RandintCommandWithArg = "randint(a,b)"
 RandomCommandWithArg = "random()"
 ReCommandWithArg = "re(z)"
+ReducedRowEchelonFormCommandWithArg = "rref(M)"
 RemCommandWithArg = "rem(p,q)"
 RootCommandWithArg = "root(x,n)"
 RoundCommandWithArg = "round(x,n)"
+RowEchelonFormCommandWithArg = "ref(M)"
 R = "r"
 Shift = "shift"
 Sigma = "σ"

apps/toolbox.de.i18n

@@ -134,6 +134,8 @@ Determinant = "Determinante"
 Transpose = "Transponierte"
 Trace = "Spur"
 Dimension = "Größe"
+RowEchelonForm = "Stufenform"
+ReducedRowEchelonForm = "Reduzierte Stufenform"
 Sort = "Sortieren aufsteigend"
 InvSort = "Sortieren absteigend"
 Maximum = "Maximalwert"

apps/toolbox.en.i18n

@@ -134,6 +134,8 @@ Determinant = "Determinant"
 Transpose = "Transpose"
 Trace = "Trace"
 Dimension = "Size"
+RowEchelonForm = "Row echelon form"
+ReducedRowEchelonForm = "Reduced row echelon form"
 Sort = "Sort ascending "
 InvSort = "Sort descending"
 Maximum = "Maximum"

apps/toolbox.es.i18n

@@ -134,6 +134,8 @@ Determinant = "Determinante"
 Transpose = "Transpuesta"
 Trace = "Traza"
 Dimension = "Tamaño"
+RowEchelonForm = "Matriz escalonada"
+ReducedRowEchelonForm = "Matriz escalonada reducida"
 Sort = "Clasificación ascendente"
 InvSort = "Clasificación descendente"
 Maximum = "Máximo"

apps/toolbox.fr.i18n

@@ -134,6 +134,8 @@ Determinant = "Déterminant de M"
 Transpose = "Transposée de M"
 Trace = "Trace de M"
 Dimension = "Taille de M"
+RowEchelonForm = "Forme échelonnée"
+ReducedRowEchelonForm = "Forme échelonnée réduite"
 Sort = "Tri croissant"
 InvSort = "Tri décroissant"
 Maximum = "Maximum"

apps/toolbox.it.i18n

@@ -134,6 +134,8 @@ Determinant = "Determinante"
 Transpose = "Trasposta"
 Trace = "Traccia"
 Dimension = "Dimensione"
+RowEchelonForm = "Matrice a scalini"
+ReducedRowEchelonForm = "Matrice ridotta a scalini"
 Sort = "Ordine crescente"
 InvSort = "Ordine decrescente"
 Maximum = "Massimo"

apps/toolbox.nl.i18n

@@ -134,6 +134,8 @@ Determinant = "Determinant"
 Transpose = "Getransponeerde"
 Trace = "Spoor"
 Dimension = "Afmeting"
+RowEchelonForm = "Echelonvorm"
+ReducedRowEchelonForm = "Gereduceerde echelonvorm"
 Sort = "Sorteer oplopend "
 InvSort = "Sort aflopend"
 Maximum = "Maximum"

apps/toolbox.pt.i18n

@@ -134,6 +134,8 @@ Determinant = "Determinante"
 Transpose = "Matriz transposta"
 Trace = "Traço"
 Dimension = "Dimensão"
+RowEchelonForm = "Matriz escalonada"
+ReducedRowEchelonForm = "Matriz escalonada reduzida"
 Sort = "Ordem crescente"
 InvSort = "Ordem decrescente"
 Maximum = "Máximo"

poincare/Makefile

@@ -101,6 +101,8 @@ poincare_src += $(addprefix poincare/src/,\
   matrix_inverse.cpp \
   matrix_trace.cpp \
   matrix_transpose.cpp \
+  matrix_ref.cpp \
+  matrix_rref.cpp \
   multiplication.cpp \
   n_ary_expression.cpp \
   naperian_logarithm.cpp \

poincare/include/poincare/expression.h

@@ -63,6 +63,8 @@ class Expression : public TreeHandle {
   friend class MatrixInverse;
   friend class MatrixTrace;
   friend class MatrixTranspose;
+  friend class MatrixRef;
+  friend class MatrixRref;
   friend class Multiplication;
   friend class MultiplicationNode;
   friend class NaperianLogarithm;

poincare/include/poincare/expression_node.h

@@ -103,6 +103,8 @@ public:
     MatrixIdentity,
     MatrixInverse,
     MatrixTranspose,
+    MatrixRef,
+    MatrixRref,
     PredictionInterval,
     Matrix,
     EmptyExpression

poincare/include/poincare/matrix.h

@@ -79,6 +79,7 @@ public:
   template<typename T> static int ArrayInverse(T * array, int numberOfRows, int numberOfColumns);
   static Matrix CreateIdentity(int dim);
   Matrix createTranspose() const;
+  Expression createRef(ExpressionNode::ReductionContext reductionContext, bool * couldComputeRef, bool reduced) const;
   /* createInverse can be called on any matrix, reduced or not, approximated or
    * not. */
   Expression createInverse(ExpressionNode::ReductionContext reductionContext, bool * couldComputeInverse) const;
@@ -94,10 +95,10 @@ private:
   void setNumberOfRows(int rows) { assert(rows >= 0); node()->setNumberOfRows(rows); }
   void setNumberOfColumns(int columns) { assert(columns >= 0); node()->setNumberOfColumns(columns); }
   Expression computeInverseOrDeterminant(bool computeDeterminant, ExpressionNode::ReductionContext reductionContext, bool * couldCompute) const;
-  // rowCanonize turns a matrix in its reduced row echelon form.
-  Matrix rowCanonize(ExpressionNode::ReductionContext reductionContext, Expression * determinant);
+  // rowCanonize turns a matrix in its row echelon form, reduced or not.
+  Matrix rowCanonize(ExpressionNode::ReductionContext reductionContext, Expression * determinant, bool reduced = true);
   // Row canonize the array in place
-  template<typename T> static void ArrayRowCanonize(T * array, int numberOfRows, int numberOfColumns, T * c = nullptr);
+  template<typename T> static void ArrayRowCanonize(T * array, int numberOfRows, int numberOfColumns, T * c = nullptr, bool reduced = true);
 
 };
 

poincare/include/poincare/matrix_complex.h

@@ -46,6 +46,7 @@ public:
   std::complex<T> determinant() const override;
   MatrixComplex<T> inverse() const;
   MatrixComplex<T> transpose() const;
+  MatrixComplex<T> ref(bool reduced) const;
 private:
   // See comment on Matrix
   uint16_t m_numberOfRows;
@@ -63,6 +64,7 @@ public:
   static MatrixComplex<T> CreateIdentity(int dim);
   MatrixComplex<T> inverse() const { return node()->inverse(); }
   MatrixComplex<T> transpose() const { return node()->transpose(); }
+  MatrixComplex<T> ref(bool reduced) const { return node()->ref(reduced); }
   std::complex<T> complexAtIndex(int index) const {
     return node()->complexAtIndex(index);
   }

poincare/include/poincare/matrix_ref.h

@@ -0,0 +1,48 @@
+#ifndef POINCARE_MATRIX_REF_H
+#define POINCARE_MATRIX_REF_H
+
+#include <poincare/expression.h>
+
+namespace Poincare {
+
+class MatrixRefNode final : public ExpressionNode {
+public:
+
+  // TreeNode
+  size_t size() const override { return sizeof(MatrixRefNode); }
+  int numberOfChildren() const override;
+#if POINCARE_TREE_LOG
+  void logNodeName(std::ostream & stream) const override {
+    stream << "MatrixRef";
+  }
+#endif
+
+  // Properties
+  Type type() const override { return Type::MatrixRef; }
+private:
+  // Layout
+  Layout createLayout(Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const override;
+  int serialize(char * buffer, int bufferSize, Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const override;
+  // Simplification
+  Expression shallowReduce(ReductionContext reductionContext) override;
+  LayoutShape leftLayoutShape() const override { return LayoutShape::MoreLetters; };
+  LayoutShape rightLayoutShape() const override { return LayoutShape::BoundaryPunctuation; }
+  // Evaluation
+  Evaluation<float> approximate(SinglePrecision p, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const override { return templatedApproximate<float>(context, complexFormat, angleUnit); }
+  Evaluation<double> approximate(DoublePrecision p, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const override { return templatedApproximate<double>(context, complexFormat, angleUnit); }
+  template<typename T> Evaluation<T> templatedApproximate(Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const;
+};
+
+class MatrixRef final : public Expression {
+public:
+  MatrixRef(const MatrixRefNode * n) : Expression(n) {}
+  static MatrixRef Builder(Expression child) { return TreeHandle::FixedArityBuilder<MatrixRef, MatrixRefNode>({child}); }
+
+  static constexpr Expression::FunctionHelper s_functionHelper = Expression::FunctionHelper("ref", 1, &UntypedBuilderOneChild<MatrixRef>);
+
+  Expression shallowReduce(ExpressionNode::ReductionContext reductionContext);
+};
+
+}
+
+#endif

poincare/include/poincare/matrix_rref.h

@@ -0,0 +1,48 @@
+#ifndef POINCARE_MATRIX_RREF_H
+#define POINCARE_MATRIX_RREF_H
+
+#include <poincare/expression.h>
+
+namespace Poincare {
+
+class MatrixRrefNode final : public ExpressionNode {
+public:
+
+  // TreeNode
+  size_t size() const override { return sizeof(MatrixRrefNode); }
+  int numberOfChildren() const override;
+#if POINCARE_TREE_LOG
+  void logNodeName(std::ostream & stream) const override {
+    stream << "MatrixRref";
+  }
+#endif
+
+  // Properties
+  Type type() const override { return Type::MatrixRref; }
+private:
+  // Layout
+  Layout createLayout(Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const override;
+  int serialize(char * buffer, int bufferSize, Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const override;
+  // Simplification
+  Expression shallowReduce(ReductionContext reductionContext) override;
+  LayoutShape leftLayoutShape() const override { return LayoutShape::MoreLetters; };
+  LayoutShape rightLayoutShape() const override { return LayoutShape::BoundaryPunctuation; }
+  // Evaluation
+  Evaluation<float> approximate(SinglePrecision p, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const override { return templatedApproximate<float>(context, complexFormat, angleUnit); }
+  Evaluation<double> approximate(DoublePrecision p, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const override { return templatedApproximate<double>(context, complexFormat, angleUnit); }
+  template<typename T> Evaluation<T> templatedApproximate(Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const;
+};
+
+class MatrixRref final : public Expression {
+public:
+  MatrixRref(const MatrixRrefNode * n) : Expression(n) {}
+  static MatrixRref Builder(Expression child) { return TreeHandle::FixedArityBuilder<MatrixRref, MatrixRrefNode>({child}); }
+
+  static constexpr Expression::FunctionHelper s_functionHelper = Expression::FunctionHelper("rref", 1, &UntypedBuilderOneChild<MatrixRref>);
+
+  Expression shallowReduce(ExpressionNode::ReductionContext reductionContext);
+};
+
+}
+
+#endif

poincare/include/poincare_nodes.h

@@ -54,6 +54,8 @@
 #include <poincare/matrix_inverse.h>
 #include <poincare/matrix_trace.h>
 #include <poincare/matrix_transpose.h>
+#include <poincare/matrix_ref.h>
+#include <poincare/matrix_rref.h>
 #include <poincare/multiplication.h>
 #include <poincare/naperian_logarithm.h>
 #include <poincare/norm_cdf.h>

poincare/src/expression.cpp

@@ -185,7 +185,9 @@ bool Expression::IsMatrix(const Expression e, Context * context) {
     || e.type() == ExpressionNode::Type::PredictionInterval
     || e.type() == ExpressionNode::Type::MatrixInverse
     || e.type() == ExpressionNode::Type::MatrixIdentity
-    || e.type() == ExpressionNode::Type::MatrixTranspose;
+    || e.type() == ExpressionNode::Type::MatrixTranspose
+    || e.type() == ExpressionNode::Type::MatrixRef
+    || e.type() == ExpressionNode::Type::MatrixRref;
 }
 
 bool Expression::IsInfinity(const Expression e, Context * context) {

poincare/src/matrix.cpp

@@ -180,7 +180,7 @@ int Matrix::ArrayInverse(T * array, int numberOfRows, int numberOfColumns) {
   return 0;
 }
 
-Matrix Matrix::rowCanonize(ExpressionNode::ReductionContext reductionContext, Expression * determinant) {
+Matrix Matrix::rowCanonize(ExpressionNode::ReductionContext reductionContext, Expression * determinant, bool reduced) {
   Expression::SetInterruption(false);
   // The matrix children have to be reduced to be able to spot 0
   deepReduceChildren(reductionContext);
@@ -231,8 +231,11 @@ Matrix Matrix::rowCanonize(ExpressionNode::ReductionContext reductionContext, Ex
       }
       replaceChildInPlace(divisor, Rational::Builder(1));
 
-      /* Set to 0 all M[i][j] i != h, j > k by linear combination */
-      for (int i = 0; i < m; i++) {
+      int l = reduced ? 0 : h + 1;
+      /* Set to 0 all M[i][j] i != h, j > k by linear combination. If a
+       * non-reduced form is computed (ref), only rows below the pivot are
+       * reduced, i > h as well */
+      for (int i = l; i < m; i++) {
         if (i == h) { continue; }
         Expression factor = matrixChild(i, k);
         for (int j = k+1; j < n; j++) {
@@ -255,7 +258,7 @@ Matrix Matrix::rowCanonize(ExpressionNode::ReductionContext reductionContext, Ex
 }
 
 template<typename T>
-void Matrix::ArrayRowCanonize(T * array, int numberOfRows, int numberOfColumns, T * determinant) {
+void Matrix::ArrayRowCanonize(T * array, int numberOfRows, int numberOfColumns, T * determinant, bool reduced) {
   int h = 0; // row pivot
   int k = 0; // column pivot
 
@@ -291,8 +294,11 @@ void Matrix::ArrayRowCanonize(T * array, int numberOfRows, int numberOfColumns,
       }
       array[h*numberOfColumns+k] = 1;
 
-      /* Set to 0 all M[i][j] i != h, j > k by linear combination */
-      for (int i = 0; i < numberOfRows; i++) {
+      int l = reduced ? 0 : h + 1;
+      /* Set to 0 all M[i][j] i != h, j > k by linear combination. If a
+       * non-reduced form is computed (ref), only rows below the pivot are
+       * reduced, i > h as well */
+      for (int i = l; i < numberOfRows; i++) {
         if (i == h) { continue; }
         T factor = array[i*numberOfColumns+k];
         for (int j = k+1; j < numberOfColumns; j++) {
@@ -330,6 +336,36 @@ Matrix Matrix::createTranspose() const {
   return matrix;
 }
 
+Expression Matrix::createRef(ExpressionNode::ReductionContext reductionContext, bool * couldComputeRef, bool reduced) const {
+  // Compute Matrix Row Echelon Form
+  /* If the matrix is too big, the rowCanonization might not be computed exactly
+   * because of a pool allocation error, but we might still be able to compute
+   * it approximately. We thus encapsulate the ref creation in an exception
+   * checkpoint.
+   * We can safely use an exception checkpoint here because we are sure of not
+   * modifying any pre-existing node in the pool. We are sure there is no Store
+   * in the matrix. */
+  Poincare::ExceptionCheckpoint ecp;
+  if (ExceptionRun(ecp)) {
+    /* We clone the current matrix to extract its children later. We can't clone
+     * its children directly. Indeed, the current matrix node (this->node()) is
+     * located before the exception checkpoint. In order to clone its chidlren,
+     * we would temporary increase the reference counter of each child (also
+     * located before the checkpoint). If an exception is raised before
+     * destroying the child handle, its reference counter would be off by one
+     * after the long jump. */
+    Matrix result = clone().convert<Matrix>();
+    *couldComputeRef = true;
+    /* Reduced row echelon form is also called row canonical form. To compute the
+     * row echelon form (non reduced one), fewer steps are required. */
+    result = result.rowCanonize(reductionContext, nullptr, reduced);
+    return std::move(result);
+  } else {
+    *couldComputeRef = false;
+    return Expression();
+  }
+}
+
 Expression Matrix::createInverse(ExpressionNode::ReductionContext reductionContext, bool * couldComputeInverse) const {
   int dim = numberOfRows();
   if (dim != numberOfColumns()) {
@@ -479,7 +515,7 @@ 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);
-template void Matrix::ArrayRowCanonize<std::complex<float> >(std::complex<float>*, int, int, std::complex<float>*);
-template void Matrix::ArrayRowCanonize<std::complex<double> >(std::complex<double>*, int, int, std::complex<double>*);
+template void Matrix::ArrayRowCanonize<std::complex<float> >(std::complex<float>*, int, int, std::complex<float>*, bool);
+template void Matrix::ArrayRowCanonize<std::complex<double> >(std::complex<double>*, int, int, std::complex<double>*, bool);
 
 }

poincare/src/matrix_complex.cpp

@@ -119,6 +119,22 @@ MatrixComplex<T> MatrixComplexNode<T>::transpose() const {
   return result;
 }
 
+template<typename T>
+MatrixComplex<T> MatrixComplexNode<T>::ref(bool reduced) const {
+  // Compute Matrix Row Echelon Form
+  if (numberOfChildren() == 0 || numberOfChildren() > Matrix::k_maxNumberOfCoefficients) {
+    return MatrixComplex<T>::Undefined();
+  }
+  std::complex<T> operandsCopy[Matrix::k_maxNumberOfCoefficients];
+  for (int i = 0; i < numberOfChildren(); i++) {
+    operandsCopy[i] = complexAtIndex(i); // Returns complex<T>(NAN, NAN) if Node type is not Complex
+  }
+  /* Reduced row echelon form is also called row canonical form. To compute the
+   * row echelon form (non reduced one), fewer steps are required. */
+  Matrix::ArrayRowCanonize(operandsCopy, m_numberOfRows, m_numberOfColumns, static_cast<std::complex<T>*>(nullptr), reduced);
+  return MatrixComplex<T>::Builder(operandsCopy, m_numberOfRows, m_numberOfColumns);
+}
+
 // MATRIX COMPLEX REFERENCE
 
 template<typename T>

poincare/src/matrix_ref.cpp

@@ -0,0 +1,61 @@
+#include <poincare/matrix_ref.h>
+#include <poincare/matrix_complex.h>
+#include <poincare/layout_helper.h>
+#include <poincare/matrix.h>
+#include <poincare/serialization_helper.h>
+
+namespace Poincare {
+
+constexpr Expression::FunctionHelper MatrixRef::s_functionHelper;
+
+int MatrixRefNode::numberOfChildren() const { return MatrixRef::s_functionHelper.numberOfChildren(); }
+
+Expression MatrixRefNode::shallowReduce(ReductionContext reductionContext) {
+  return MatrixRef(this).shallowReduce(reductionContext);
+}
+
+Layout MatrixRefNode::createLayout(Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const {
+  return LayoutHelper::Prefix(MatrixRef(this), floatDisplayMode, numberOfSignificantDigits, MatrixRef::s_functionHelper.name());
+}
+
+int MatrixRefNode::serialize(char * buffer, int bufferSize, Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const {
+  return SerializationHelper::Prefix(this, buffer, bufferSize, floatDisplayMode, numberOfSignificantDigits, MatrixRef::s_functionHelper.name());
+}
+
+template<typename T>
+Evaluation<T> MatrixRefNode::templatedApproximate(Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const {
+  Evaluation<T> input = childAtIndex(0)->approximate(T(), context, complexFormat, angleUnit);
+  Evaluation<T> ref;
+  if (input.type() == EvaluationNode<T>::Type::MatrixComplex) {
+    ref = static_cast<MatrixComplex<T>&>(input).ref(false);
+  } else {
+    ref = Complex<T>::Undefined();
+  }
+  assert(!ref.isUninitialized());
+  return ref;
+}
+
+
+Expression MatrixRef::shallowReduce(ExpressionNode::ReductionContext reductionContext) {
+  {
+    Expression e = Expression::defaultShallowReduce();
+    e = e.defaultHandleUnitsInChildren();
+    if (e.isUndefined()) {
+      return e;
+    }
+  }
+  Expression c = childAtIndex(0);
+  if (c.type() == ExpressionNode::Type::Matrix) {
+    bool couldComputeRef = false;
+    Expression result = static_cast<Matrix&>(c).createRef(reductionContext, &couldComputeRef, false);
+    if (couldComputeRef) {
+      replaceWithInPlace(result);
+      return result;
+    }
+    // The matrix could not be transformed properly
+    return *this;
+  }
+  return replaceWithUndefinedInPlace();
+}
+
+}

poincare/src/matrix_rref.cpp

@@ -0,0 +1,61 @@
+#include <poincare/matrix_rref.h>
+#include <poincare/matrix_complex.h>
+#include <poincare/layout_helper.h>
+#include <poincare/matrix.h>
+#include <poincare/serialization_helper.h>
+
+namespace Poincare {
+
+constexpr Expression::FunctionHelper MatrixRref::s_functionHelper;
+
+int MatrixRrefNode::numberOfChildren() const { return MatrixRref::s_functionHelper.numberOfChildren(); }
+
+Expression MatrixRrefNode::shallowReduce(ReductionContext reductionContext) {
+  return MatrixRref(this).shallowReduce(reductionContext);
+}
+
+Layout MatrixRrefNode::createLayout(Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const {
+  return LayoutHelper::Prefix(MatrixRref(this), floatDisplayMode, numberOfSignificantDigits, MatrixRref::s_functionHelper.name());
+}
+
+int MatrixRrefNode::serialize(char * buffer, int bufferSize, Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const {
+  return SerializationHelper::Prefix(this, buffer, bufferSize, floatDisplayMode, numberOfSignificantDigits, MatrixRref::s_functionHelper.name());
+}
+
+template<typename T>
+Evaluation<T> MatrixRrefNode::templatedApproximate(Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const {
+  Evaluation<T> input = childAtIndex(0)->approximate(T(), context, complexFormat, angleUnit);
+  Evaluation<T> rref;
+  if (input.type() == EvaluationNode<T>::Type::MatrixComplex) {
+    rref = static_cast<MatrixComplex<T>&>(input).ref(true);
+  } else {
+    rref = Complex<T>::Undefined();
+  }
+  assert(!rref.isUninitialized());
+  return rref;
+}
+
+
+Expression MatrixRref::shallowReduce(ExpressionNode::ReductionContext reductionContext) {
+  {
+    Expression e = Expression::defaultShallowReduce();
+    e = e.defaultHandleUnitsInChildren();
+    if (e.isUndefined()) {
+      return e;
+    }
+  }
+  Expression c = childAtIndex(0);
+  if (c.type() == ExpressionNode::Type::Matrix) {
+    bool couldComputeRref = false;
+    Expression result = static_cast<Matrix&>(c).createRef(reductionContext, &couldComputeRref, true);
+    if (couldComputeRref) {
+      replaceWithInPlace(result);
+      return result;
+    }
+    // The matrix could not be transformed properly
+    return *this;
+  }
+  return replaceWithUndefinedInPlace();
+}
+
+}

poincare/src/parsing/parser.h

@@ -138,9 +138,11 @@ private:
     &Randint::s_functionHelper,
     &Random::s_functionHelper,
     &RealPart::s_functionHelper,
+    &MatrixRef::s_functionHelper,
     &DivisionRemainder::s_functionHelper,
     &NthRoot::s_functionHelper,
     &Round::s_functionHelper,
+    &MatrixRref::s_functionHelper,
     &SignFunction::s_functionHelper,
     &Sine::s_functionHelper,
     &HyperbolicSine::s_functionHelper,

poincare/src/tree_handle.cpp

@@ -336,6 +336,8 @@ template MatrixIdentity TreeHandle::FixedArityBuilder<MatrixIdentity, MatrixIden
 template MatrixInverse TreeHandle::FixedArityBuilder<MatrixInverse, MatrixInverseNode>(const Tuple &);
 template MatrixTrace TreeHandle::FixedArityBuilder<MatrixTrace, MatrixTraceNode>(const Tuple &);
 template MatrixTranspose TreeHandle::FixedArityBuilder<MatrixTranspose, MatrixTransposeNode>(const Tuple &);
+template MatrixRef TreeHandle::FixedArityBuilder<MatrixRef, MatrixRefNode>(const Tuple &);
+template MatrixRref TreeHandle::FixedArityBuilder<MatrixRref, MatrixRrefNode>(const Tuple &);
 template Multiplication TreeHandle::NAryBuilder<Multiplication, MultiplicationNode>(const Tuple &);
 template NaperianLogarithm TreeHandle::FixedArityBuilder<NaperianLogarithm, NaperianLogarithmNode>(const Tuple &);
 template NormCDF TreeHandle::FixedArityBuilder<NormCDF, NormCDFNode>(const Tuple &);

poincare/test/approximation.cpp

@@ -405,6 +405,25 @@ QUIZ_CASE(poincare_approximation_function) {
   assert_expression_approximates_to<double>("transpose([[1,7,5][4,2,8]])", "[[1,4][7,2][5,8]]");
   assert_expression_approximates_to<double>("transpose([[1,2][4,5][7,8]])", "[[1,4,7][2,5,8]]");
 
+  /* Results for ref depend on the implementation. In any case :
+   * - 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]]"
+  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>("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>("rref([[0,2,-1][5,6,7][12,11,10]])", "[[1,0,0][0,1,0][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/expression_properties.cpp

@@ -63,6 +63,8 @@ QUIZ_CASE(poincare_properties_is_matrix) {
   assert_expression_has_property("inverse([[1,2][3,4]])", &context, Expression::IsMatrix);
   assert_expression_has_property("3*identity(4)", &context, Expression::IsMatrix);
   assert_expression_has_property("transpose([[1,2][3,4]])", &context, Expression::IsMatrix);
+  assert_expression_has_property("ref([[1,2][3,4]])", &context, Expression::IsMatrix);
+  assert_expression_has_property("rref([[1,2][3,4]])", &context, Expression::IsMatrix);
   assert_expression_has_not_property("2*3+1", &context, Expression::IsMatrix);
 }
 

poincare/test/parsing.cpp

@@ -410,6 +410,8 @@ QUIZ_CASE(poincare_parsing_identifiers) {
   assert_parsed_expression_is("tanh(1)", HyperbolicTangent::Builder(BasedInteger::Builder(1)));
   assert_parsed_expression_is("trace(1)", MatrixTrace::Builder(BasedInteger::Builder(1)));
   assert_parsed_expression_is("transpose(1)", MatrixTranspose::Builder(BasedInteger::Builder(1)));
+  assert_parsed_expression_is("ref(1)", MatrixRef::Builder(BasedInteger::Builder(1)));
+  assert_parsed_expression_is("rref(1)", MatrixRref::Builder(BasedInteger::Builder(1)));
   assert_parsed_expression_is("√(1)", SquareRoot::Builder(BasedInteger::Builder(1)));
   assert_text_not_parsable("cos(1,2)");
   assert_text_not_parsable("log(1,2,3)");

poincare/test/simplification.cpp

@@ -447,6 +447,8 @@ QUIZ_CASE(poincare_simplification_units) {
   assert_parsed_expression_simplify_to("tanh(_s)", "undef");
   assert_parsed_expression_simplify_to("trace(_s)", "undef");
   assert_parsed_expression_simplify_to("transpose(_s)", "undef");
+  assert_parsed_expression_simplify_to("ref(_s)", "undef");
+  assert_parsed_expression_simplify_to("rref(_s)", "undef");
 
   /* Valid expressions */
   assert_parsed_expression_simplify_to("-2×_A", "-2×_A");
@@ -927,6 +929,12 @@ QUIZ_CASE(poincare_simplification_matrix) {
   assert_parsed_expression_simplify_to("transpose([[1/√(2),1/2,3][2,1,-3]])", "[[√(2)/2,2][1/2,1][3,-3]]");
   assert_parsed_expression_simplify_to("transpose(√(4))", "2");
 
+  // Ref and Rref
+  assert_parsed_expression_simplify_to("ref([[1,1/√(2),√(4)]])", "[[1,√(2)/2,2]]");
+  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]]");
+
   // Expressions with unreduced matrix
   assert_reduce("confidence(cos(2)/25,3)→a");
   // Check that matrices are not permuted in multiplication