diff --git a/poincare/include/poincare/binomial_coefficient.h b/poincare/include/poincare/binomial_coefficient.h
index 5f7094346..801bde559 100644
--- a/poincare/include/poincare/binomial_coefficient.h
+++ b/poincare/include/poincare/binomial_coefficient.h
@@ -43,7 +43,7 @@ public:
   static constexpr Expression::FunctionHelper s_functionHelper = Expression::FunctionHelper("binomial", 2, &UntypedBuilderTwoChildren<BinomialCoefficient>);
 
   // Expression
-  Expression shallowReduce();
+  Expression shallowReduce(Context & context);
 private:
   constexpr static int k_maxNValue = 300;
 };
diff --git a/poincare/include/poincare/ceiling.h b/poincare/include/poincare/ceiling.h
index 707a1c9ea..943542010 100644
--- a/poincare/include/poincare/ceiling.h
+++ b/poincare/include/poincare/ceiling.h
@@ -46,7 +46,7 @@ public:
 
   static constexpr Expression::FunctionHelper s_functionHelper = Expression::FunctionHelper("ceil", 1, &UntypedBuilderOneChild<Ceiling>);
 
-  Expression shallowReduce();
+  Expression shallowReduce(Context & context, Preferences::AngleUnit angleUnit);
 };
 
 }
diff --git a/poincare/include/poincare/derivative.h b/poincare/include/poincare/derivative.h
index 8d26ecf0d..04f5ab069 100644
--- a/poincare/include/poincare/derivative.h
+++ b/poincare/include/poincare/derivative.h
@@ -54,7 +54,7 @@ public:
   static Expression UntypedBuilder(Expression children);
   static constexpr Expression::FunctionHelper s_functionHelper = Expression::FunctionHelper("diff", 3, &UntypedBuilder);
 
-  Expression shallowReduce();
+  Expression shallowReduce(Context & context);
 };
 
 }
diff --git a/poincare/include/poincare/expression.h b/poincare/include/poincare/expression.h
index 2d3a09fe8..49b4b1f2b 100644
--- a/poincare/include/poincare/expression.h
+++ b/poincare/include/poincare/expression.h
@@ -133,9 +133,11 @@ public:
   typedef bool (*ExpressionTest)(const Expression e, Context & context);
   bool recursivelyMatches(ExpressionTest test, Context & context, bool replaceSymbols = true) const;
   // Set of ExpressionTest that can be used with recursivelyMatches
+  static bool IsNAry(const Expression e, Context & context);
   static bool IsApproximate(const Expression e, Context & context);
   static bool IsRandom(const Expression e, Context & context);
   static bool IsMatrix(const Expression e, Context & context);
+  static bool SortedIsMatrix(const Expression e, Context & context);
   static bool IsInfinity(const Expression e, Context & context);
   /* 'characteristicXRange' tries to assess the range on x where the expression
    * (considered as a function on x) has an interesting evolution. For example,
diff --git a/poincare/include/poincare/factorial.h b/poincare/include/poincare/factorial.h
index b00826768..34c0fbc52 100644
--- a/poincare/include/poincare/factorial.h
+++ b/poincare/include/poincare/factorial.h
@@ -54,7 +54,7 @@ public:
   Factorial(const FactorialNode * n) : Expression(n) {}
   static Factorial Builder(Expression child) { return TreeHandle::FixedArityBuilder<Factorial, FactorialNode>(&child, 1); }
 
-  Expression shallowReduce();
+  Expression shallowReduce(Context & context, Preferences::AngleUnit angleUnit);
   Expression shallowBeautify();
 private:
   constexpr static int k_maxOperandValue = 100;
diff --git a/poincare/include/poincare/floor.h b/poincare/include/poincare/floor.h
index e1d4e3ace..f2db0393a 100644
--- a/poincare/include/poincare/floor.h
+++ b/poincare/include/poincare/floor.h
@@ -48,7 +48,7 @@ public:
 
   static constexpr Expression::FunctionHelper s_functionHelper = Expression::FunctionHelper("floor", 1, &UntypedBuilderOneChild<Floor>);
 
-  Expression shallowReduce();
+  Expression shallowReduce(Context & context, Preferences::AngleUnit angleUnit);
 };
 
 }
diff --git a/poincare/include/poincare/frac_part.h b/poincare/include/poincare/frac_part.h
index 09d5637c8..9c4c8a559 100644
--- a/poincare/include/poincare/frac_part.h
+++ b/poincare/include/poincare/frac_part.h
@@ -48,7 +48,7 @@ public:
 
   static constexpr Expression::FunctionHelper s_functionHelper = Expression::FunctionHelper("frac", 1, &UntypedBuilderOneChild<FracPart>);
 
-  Expression shallowReduce();
+  Expression shallowReduce(Context & context, Preferences::AngleUnit angleUnit);
 };
 
 }
diff --git a/poincare/include/poincare/hyperbolic_trigonometric_function.h b/poincare/include/poincare/hyperbolic_trigonometric_function.h
index 58854c30a..228d62f93 100644
--- a/poincare/include/poincare/hyperbolic_trigonometric_function.h
+++ b/poincare/include/poincare/hyperbolic_trigonometric_function.h
@@ -18,7 +18,7 @@ private:
 class HyperbolicTrigonometricFunction : public Expression {
 public:
   HyperbolicTrigonometricFunction(const HyperbolicTrigonometricFunctionNode * n) : Expression(n) {}
-  Expression shallowReduce();
+  Expression shallowReduce(Context & context, Preferences::AngleUnit angleUnit);
 };
 
 }
diff --git a/poincare/include/poincare/integral.h b/poincare/include/poincare/integral.h
index d848e05b1..1f8d31dd6 100644
--- a/poincare/include/poincare/integral.h
+++ b/poincare/include/poincare/integral.h
@@ -60,7 +60,7 @@ public:
   static constexpr Expression::FunctionHelper s_functionHelper = Expression::FunctionHelper("int", 4, &UntypedBuilder);
 
   // Expression
-  Expression shallowReduce();
+  Expression shallowReduce(Context & context);
 };
 
 }
diff --git a/poincare/include/poincare/n_ary_expression.h b/poincare/include/poincare/n_ary_expression.h
index bc30ee7a1..94008a1c6 100644
--- a/poincare/include/poincare/n_ary_expression.h
+++ b/poincare/include/poincare/n_ary_expression.h
@@ -10,6 +10,7 @@ namespace Poincare {
 class NAryExpressionNode : public ExpressionNode { // TODO: VariableArityExpressionNode?
 public:
   void setChildrenInPlace(Expression other) override { assert(false); }
+  static bool IsMatrix(Expression e, Context & context);
 
   //Tree
   int numberOfChildren() const override { return m_numberOfChildren; }
@@ -58,6 +59,7 @@ public:
    * - 0 if all non real children are ComplexCartesian
    * - -1 if some chidren are non-real and non ComplexCartesian */
   int allChildrenAreReal(Context & context) const;
+  static bool SortedIsMatrix(Expression e, Context & context); // this is supposed to be a sorted
 protected:
   NAryExpressionNode * node() const { return static_cast<NAryExpressionNode *>(Expression::node()); }
 };
diff --git a/poincare/src/binomial_coefficient.cpp b/poincare/src/binomial_coefficient.cpp
index ada1b413f..0dc83df11 100644
--- a/poincare/src/binomial_coefficient.cpp
+++ b/poincare/src/binomial_coefficient.cpp
@@ -15,7 +15,7 @@ constexpr Expression::FunctionHelper BinomialCoefficient::s_functionHelper;
 int BinomialCoefficientNode::numberOfChildren() const { return BinomialCoefficient::s_functionHelper.numberOfChildren(); }
 
 Expression BinomialCoefficientNode::shallowReduce(Context & context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, ReductionTarget target, bool symbolicComputation) {
-  return BinomialCoefficient(this).shallowReduce();
+  return BinomialCoefficient(this).shallowReduce(context);
 }
 
 Layout BinomialCoefficientNode::createLayout(Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const {
@@ -54,7 +54,7 @@ T BinomialCoefficientNode::compute(T k, T n) {
 }
 
 
-Expression BinomialCoefficient::shallowReduce() {
+Expression BinomialCoefficient::shallowReduce(Context & context) {
   {
     Expression e = Expression::defaultShallowReduce();
     if (e.isUndefined()) {
@@ -63,11 +63,11 @@ Expression BinomialCoefficient::shallowReduce() {
   }
   Expression c0 = childAtIndex(0);
   Expression c1 = childAtIndex(1);
-#if MATRIX_EXACT_REDUCING
-  if (c0.type() == ExpressionNode::Type::Matrix || c1.type() == ExpressionNode::Type::Matrix) {
+
+  if (SortedIsMatrix(c0, context) || SortedIsMatrix(c1, context)) {
     return Undefined::Builder();
   }
-#endif
+
   if (c0.type() == ExpressionNode::Type::Rational) {
     Rational r0 = static_cast<Rational&>(c0);
     if (!r0.integerDenominator().isOne() || r0.isNegative()) {
diff --git a/poincare/src/ceiling.cpp b/poincare/src/ceiling.cpp
index 907d48767..d99f60d20 100644
--- a/poincare/src/ceiling.cpp
+++ b/poincare/src/ceiling.cpp
@@ -32,11 +32,11 @@ Complex<T> CeilingNode::computeOnComplex(const std::complex<T> c, Preferences::C
 }
 
 Expression CeilingNode::shallowReduce(Context & context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, ReductionTarget target, bool symbolicComputation) {
-  return Ceiling(this).shallowReduce();
+  return Ceiling(this).shallowReduce(context, angleUnit);
 }
 
 
-Expression Ceiling::shallowReduce() {
+Expression Ceiling::shallowReduce(Context & context, Preferences::AngleUnit angleUnit) {
   {
     Expression e = Expression::defaultShallowReduce();
     if (e.isUndefined()) {
@@ -44,11 +44,9 @@ Expression Ceiling::shallowReduce() {
     }
   }
   Expression c = childAtIndex(0);
-#if MATRIX_EXACT_REDUCING
   if (c.type() == ExpressionNode::Type::Matrix) {
     return SimplificationHelper::Map(*this, context, angleUnit);
   }
-#endif
   if (c.type() == ExpressionNode::Type::Constant) {
     Constant s = static_cast<Constant&>(c);
     Expression result;
diff --git a/poincare/src/complex_argument.cpp b/poincare/src/complex_argument.cpp
index 13dccd0e5..34b1dd7bd 100644
--- a/poincare/src/complex_argument.cpp
+++ b/poincare/src/complex_argument.cpp
@@ -42,11 +42,9 @@ Expression ComplexArgument::shallowReduce(Context & context, Preferences::Comple
     }
   }
   Expression c = childAtIndex(0);
-#if MATRIX_EXACT_REDUCING
   if (c.type() == ExpressionNode::Type::Matrix) {
     return SimplificationHelper::Map(*this, context, angleUnit);
   }
-#endif
   bool real = c.isReal(context);
   if (real) {
     float app = c.node()->approximate(float(), context, complexFormat, angleUnit).toScalar();
diff --git a/poincare/src/complex_cartesian.cpp b/poincare/src/complex_cartesian.cpp
index 000156c78..ecee23825 100644
--- a/poincare/src/complex_cartesian.cpp
+++ b/poincare/src/complex_cartesian.cpp
@@ -265,7 +265,7 @@ ComplexCartesian ComplexCartesian::powerInteger(int n, Context & context, Prefer
     Power apow = Power::Builder(aclone, Rational::Builder(n-i));
     Power bpow = Power::Builder(bclone, Rational::Builder(i));
     Multiplication m = Multiplication::Builder(binom, apow, bpow);
-    binom.shallowReduce();
+    binom.shallowReduce(context);
     apow.shallowReduce(context, complexFormat, angleUnit, target);
     bpow.shallowReduce(context, complexFormat, angleUnit, target);
     if (i/2%2 == 1) {
diff --git a/poincare/src/confidence_interval.cpp b/poincare/src/confidence_interval.cpp
index 68bfc3d87..bb843ae95 100644
--- a/poincare/src/confidence_interval.cpp
+++ b/poincare/src/confidence_interval.cpp
@@ -62,11 +62,9 @@ Expression ConfidenceInterval::shallowReduce(Context & context, Preferences::Com
   }
   Expression c0 = childAtIndex(0);
   Expression c1 = childAtIndex(1);
-#if MATRIX_EXACT_REDUCING
-  if (c0.type() == ExpressionNode::Type::Matrix || c1.type() == ExpressionNode::Type::Matrix) {
+  if (SortedIsMatrix(c0, context) || SortedIsMatrix(c1, context)) {
     return Undefined::Builder();
   }
-#endif
   if (c0.type() == ExpressionNode::Type::Rational) {
     Rational r0 = static_cast<Rational&>(c0);
     if (r0.signedIntegerNumerator().isNegative() || Integer::NaturalOrder(r0.signedIntegerNumerator(), r0.integerDenominator()) > 0) {
diff --git a/poincare/src/conjugate.cpp b/poincare/src/conjugate.cpp
index 9d3f579e3..686c4e293 100644
--- a/poincare/src/conjugate.cpp
+++ b/poincare/src/conjugate.cpp
@@ -40,11 +40,9 @@ Expression Conjugate::shallowReduce(Context & context, Preferences::ComplexForma
     }
   }
   Expression c = childAtIndex(0);
-#if MATRIX_EXACT_REDUCING
   if (c.type() == ExpressionNode::Type::Matrix) {
     return SimplificationHelper::Map(*this, context, angleUnit);
   }
-#endif
   if (c.isReal(context)) {
     replaceWithInPlace(c);
     return c;
diff --git a/poincare/src/cosine.cpp b/poincare/src/cosine.cpp
index a7063d9dc..880fcc863 100644
--- a/poincare/src/cosine.cpp
+++ b/poincare/src/cosine.cpp
@@ -41,12 +41,10 @@ Expression Cosine::shallowReduce(Context & context, Preferences::ComplexFormat c
       return e;
     }
   }
-#if MATRIX_EXACT_REDUCING
-  Expression op = childAtIndex(0);
-  if (op.type() == ExpressionNode::Type::Matrix) {
+  Expression c = childAtIndex(0);
+  if (c.type() == ExpressionNode::Type::Matrix) {
     return SimplificationHelper::Map(*this, context, angleUnit);
   }
-#endif
   return Trigonometry::shallowReduceDirectFunction(*this, context, complexFormat, angleUnit, target);
 }
 
diff --git a/poincare/src/derivative.cpp b/poincare/src/derivative.cpp
index 03115763c..ea8a60542 100644
--- a/poincare/src/derivative.cpp
+++ b/poincare/src/derivative.cpp
@@ -35,7 +35,7 @@ int DerivativeNode::serialize(char * buffer, int bufferSize, Preferences::PrintF
 }
 
 Expression DerivativeNode::shallowReduce(Context & context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, ReductionTarget target, bool symbolicComputation) {
-  return Derivative(this).shallowReduce();
+  return Derivative(this).shallowReduce(context);
 }
 
 template<typename T>
@@ -135,18 +135,19 @@ T DerivativeNode::riddersApproximation(Context & context, Preferences::ComplexFo
   return ans;
 }
 
-Expression Derivative::shallowReduce() {
+Expression Derivative::shallowReduce(Context & context) {
   {
     Expression e = Expression::defaultShallowReduce();
     if (e.isUndefined()) {
       return e;
     }
   }
-#if MATRIX_EXACT_REDUCING
-  if (childAtIndex(0).type() == ExpressionNode::Type::Matrix || || childAtIndex(1).type() == ExpressionNode::Type::Matrix || childAtIndex(2).type() == ExpressionNode::Type::Matrix) {
+  if (SortedIsMatrix(childAtIndex(0), context)
+      || SortedIsMatrix(childAtIndex(1), context)
+      || SortedIsMatrix(childAtIndex(2), context))
+  {
     return Undefined::Builder();
   }
-#endif
   // TODO: to be implemented diff(+) -> +diff() etc
   return *this;
 }
diff --git a/poincare/src/determinant.cpp b/poincare/src/determinant.cpp
index e41a2c6ea..aac6f4e2d 100644
--- a/poincare/src/determinant.cpp
+++ b/poincare/src/determinant.cpp
@@ -41,19 +41,12 @@ Expression Determinant::shallowReduce(Context & context) {
     }
   }
   Expression c0 = childAtIndex(0);
-#if MATRIX_EXACT_REDUCING
-#if 0
-  if (!op.recursivelyMatches(Expression::IsMatrix)) {
-    return replaceWith(op, true);
-  }
-  return this;
-#endif
-#endif
   // det(A) = A if A is not a matrix
-  if (!c0.recursivelyMatches(Expression::IsMatrix, context, true)) {
+  if (!SortedIsMatrix(c0, context)) {
     replaceWithInPlace(c0);
     return c0;
   }
+  //TODO LEA for matrix
   return *this;
 }
 
diff --git a/poincare/src/expression.cpp b/poincare/src/expression.cpp
index 1b3618bb0..d711cf32b 100644
--- a/poincare/src/expression.cpp
+++ b/poincare/src/expression.cpp
@@ -106,6 +106,10 @@ bool Expression::IsRandom(const Expression e, Context & context) {
   return e.isRandom();
 }
 
+bool Expression::IsNAry(const Expression e, Context & context) {
+  return e.type() == ExpressionNode::Type::Addition || e.type() == ExpressionNode::Type::Multiplication;
+}
+
 bool Expression::IsMatrix(const Expression e, Context & context) {
   return e.type() == ExpressionNode::Type::Matrix
     || e.type() == ExpressionNode::Type::ConfidenceInterval
@@ -116,6 +120,17 @@ bool Expression::IsMatrix(const Expression e, Context & context) {
     || e.type() == ExpressionNode::Type::MatrixTranspose;
 }
 
+bool Expression::SortedIsMatrix(const Expression e, Context & context) {
+  // TODO should we replace symbols?
+  if (IsMatrix(e, context)) {
+    return true;
+  }
+  if (IsNAry(e, context)) {
+    return NAryExpression::SortedIsMatrix(e, context);
+  }
+  return false;
+}
+
 bool Expression::IsInfinity(const Expression e, Context & context) {
   return e.type() == ExpressionNode::Type::Infinity;
 }
diff --git a/poincare/src/factor.cpp b/poincare/src/factor.cpp
index 910ce935c..f58d7944b 100644
--- a/poincare/src/factor.cpp
+++ b/poincare/src/factor.cpp
@@ -67,6 +67,7 @@ Expression Factor::shallowBeautify(Context & context, Preferences::ComplexFormat
 }
 
 Multiplication Factor::createMultiplicationOfIntegerPrimeDecomposition(Integer i, Context & context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const {
+  //TODO LEA Matrices?
   assert(!i.isZero());
   assert(!i.isNegative());
   Multiplication m = Multiplication::Builder();
diff --git a/poincare/src/factorial.cpp b/poincare/src/factorial.cpp
index 402b54bd8..bfa3afa88 100644
--- a/poincare/src/factorial.cpp
+++ b/poincare/src/factorial.cpp
@@ -35,7 +35,7 @@ bool FactorialNode::childNeedsParenthesis(const TreeNode * child) const {
 // Simplification
 
 Expression FactorialNode::shallowReduce(Context & context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, ReductionTarget target, bool symbolicComputation) {
-  return Factorial(this).shallowReduce();
+  return Factorial(this).shallowReduce(context, angleUnit);
 }
 
 Expression FactorialNode::shallowBeautify(Context & context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, ReductionTarget target) {
@@ -77,20 +77,19 @@ int FactorialNode::serialize(char * buffer, int bufferSize, Preferences::PrintFl
 }
 
 
-Expression Factorial::shallowReduce() {
+Expression Factorial::shallowReduce(Context & context, Preferences::AngleUnit angleUnit) {
   {
     Expression e = Expression::defaultShallowReduce();
     if (e.isUndefined()) {
       return e;
     }
   }
-#if MATRIX_EXACT_REDUCING
-  if (childAtIndex(0).type() == ExpressionNode::Type::Matrix) {
+  Expression c = childAtIndex(0);
+  if (c.type() == ExpressionNode::Type::Matrix) {
     return SimplificationHelper::Map(*this, context, angleUnit);
   }
-#endif
-  if (childAtIndex(0).type() == ExpressionNode::Type::Rational) {
-    Rational r = childAtIndex(0).convert<Rational>();
+  if (c.type() == ExpressionNode::Type::Rational) {
+    Rational r = c.convert<Rational>();
     if (!r.integerDenominator().isOne() || r.sign() == ExpressionNode::Sign::Negative) {
       Expression result = Undefined::Builder();
       replaceWithInPlace(result);
@@ -104,7 +103,7 @@ Expression Factorial::shallowReduce() {
     replaceWithInPlace(fact);
     return fact;
   }
-  if (childAtIndex(0).type() == ExpressionNode::Type::Constant) {
+  if (c.type() == ExpressionNode::Type::Constant) {
     // e! = undef, i! = undef, pi! = undef
     Expression result = Undefined::Builder();
     replaceWithInPlace(result);
diff --git a/poincare/src/floor.cpp b/poincare/src/floor.cpp
index b5dd5a21e..058f69cca 100644
--- a/poincare/src/floor.cpp
+++ b/poincare/src/floor.cpp
@@ -32,11 +32,10 @@ Complex<T> FloorNode::computeOnComplex(const std::complex<T> c, Preferences::Com
 }
 
 Expression FloorNode::shallowReduce(Context & context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, ReductionTarget target, bool symbolicComputation) {
-  return Floor(this).shallowReduce();
+  return Floor(this).shallowReduce(context, angleUnit);
 }
 
-
-Expression Floor::shallowReduce() {
+Expression Floor::shallowReduce(Context & context, Preferences::AngleUnit angleUnit) {
   {
     Expression e = Expression::defaultShallowReduce();
     if (e.isUndefined()) {
@@ -44,11 +43,9 @@ Expression Floor::shallowReduce() {
     }
   }
   Expression c = childAtIndex(0);
-#if MATRIX_EXACT_REDUCING
   if (c.type() == ExpressionNode::Type::Matrix) {
     return SimplificationHelper::Map(*this, context, angleUnit);
   }
-#endif
   if (c.type() == ExpressionNode::Type::Constant) {
     Constant s = static_cast<Constant &>(c);
     Expression result;
diff --git a/poincare/src/frac_part.cpp b/poincare/src/frac_part.cpp
index 4d804aa77..dd2e7f9c9 100644
--- a/poincare/src/frac_part.cpp
+++ b/poincare/src/frac_part.cpp
@@ -20,7 +20,7 @@ int FracPartNode::serialize(char * buffer, int bufferSize, Preferences::PrintFlo
 }
 
 Expression FracPartNode::shallowReduce(Context & context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, ReductionTarget target, bool symbolicComputation) {
-  return FracPart(this).shallowReduce();
+  return FracPart(this).shallowReduce(context, angleUnit);
 }
 
 template<typename T>
@@ -32,7 +32,7 @@ Complex<T> FracPartNode::computeOnComplex(const std::complex<T> c, Preferences::
 }
 
 
-Expression FracPart::shallowReduce() {
+Expression FracPart::shallowReduce(Context & context, Preferences::AngleUnit angleUnit) {
   {
     Expression e = Expression::defaultShallowReduce();
     if (e.isUndefined()) {
@@ -40,11 +40,9 @@ Expression FracPart::shallowReduce() {
     }
   }
   Expression c = childAtIndex(0);
-#if MATRIX_EXACT_REDUCING
   if (c.type() == ExpressionNode::Type::Matrix) {
     return SimplificationHelper::Map(*this, context, angleUnit);
   }
-#endif
   if (c.type() != ExpressionNode::Type::Rational) {
     return *this;
   }
diff --git a/poincare/src/hyperbolic_arc_cosine.cpp b/poincare/src/hyperbolic_arc_cosine.cpp
index b52aa2383..8ea38dd6c 100644
--- a/poincare/src/hyperbolic_arc_cosine.cpp
+++ b/poincare/src/hyperbolic_arc_cosine.cpp
@@ -26,7 +26,6 @@ Complex<T> HyperbolicArcCosineNode::computeOnComplex(const std::complex<T> c, Pr
   return Complex<T>::Builder(Trigonometry::RoundToMeaningfulDigits(result, c));
 }
 
-
 template Complex<float> Poincare::HyperbolicArcCosineNode::computeOnComplex<float>(std::complex<float>, Preferences::ComplexFormat, Preferences::AngleUnit);
 template Complex<double> Poincare::HyperbolicArcCosineNode::computeOnComplex<double>(std::complex<double>, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit);
 
diff --git a/poincare/src/hyperbolic_arc_sine.cpp b/poincare/src/hyperbolic_arc_sine.cpp
index a2c95af34..7ee51048c 100644
--- a/poincare/src/hyperbolic_arc_sine.cpp
+++ b/poincare/src/hyperbolic_arc_sine.cpp
@@ -11,6 +11,7 @@ constexpr Expression::FunctionHelper HyperbolicArcSine::s_functionHelper;
 Layout HyperbolicArcSineNode::createLayout(Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const {
   return LayoutHelper::Prefix(HyperbolicArcSine(this), floatDisplayMode, numberOfSignificantDigits, HyperbolicArcSine::s_functionHelper.name());
 }
+
 int HyperbolicArcSineNode::serialize(char * buffer, int bufferSize, Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const {
   return SerializationHelper::Prefix(this, buffer, bufferSize, floatDisplayMode, numberOfSignificantDigits, HyperbolicArcSine::s_functionHelper.name());
 }
@@ -29,7 +30,6 @@ Complex<T> HyperbolicArcSineNode::computeOnComplex(const std::complex<T> c, Pref
   return Complex<T>::Builder(Trigonometry::RoundToMeaningfulDigits(result, c));
 }
 
-
 template Complex<float> Poincare::HyperbolicArcSineNode::computeOnComplex<float>(std::complex<float>, Preferences::ComplexFormat, Preferences::AngleUnit);
 template Complex<double> Poincare::HyperbolicArcSineNode::computeOnComplex<double>(std::complex<double>, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit);
 
diff --git a/poincare/src/hyperbolic_arc_tangent.cpp b/poincare/src/hyperbolic_arc_tangent.cpp
index 099835978..13c63bba5 100644
--- a/poincare/src/hyperbolic_arc_tangent.cpp
+++ b/poincare/src/hyperbolic_arc_tangent.cpp
@@ -30,7 +30,6 @@ Complex<T> HyperbolicArcTangentNode::computeOnComplex(const std::complex<T> c, P
   return Complex<T>::Builder(Trigonometry::RoundToMeaningfulDigits(result, c));
 }
 
-
 template Complex<float> Poincare::HyperbolicArcTangentNode::computeOnComplex<float>(std::complex<float>, Preferences::ComplexFormat, Preferences::AngleUnit);
 template Complex<double> Poincare::HyperbolicArcTangentNode::computeOnComplex<double>(std::complex<double>, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit);
 
diff --git a/poincare/src/hyperbolic_cosine.cpp b/poincare/src/hyperbolic_cosine.cpp
index 30f564714..c722ec209 100644
--- a/poincare/src/hyperbolic_cosine.cpp
+++ b/poincare/src/hyperbolic_cosine.cpp
@@ -9,6 +9,7 @@ constexpr Expression::FunctionHelper HyperbolicCosine::s_functionHelper;
 Layout HyperbolicCosineNode::createLayout(Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const {
   return LayoutHelper::Prefix(HyperbolicCosine(this), floatDisplayMode, numberOfSignificantDigits, HyperbolicCosine::s_functionHelper.name());
 }
+
 int HyperbolicCosineNode::serialize(char * buffer, int bufferSize, Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const {
   return SerializationHelper::Prefix(this, buffer, bufferSize, floatDisplayMode, numberOfSignificantDigits, HyperbolicCosine::s_functionHelper.name());
 }
@@ -18,7 +19,6 @@ Complex<T> HyperbolicCosineNode::computeOnComplex(const std::complex<T> c, Prefe
   return Complex<T>::Builder(Trigonometry::RoundToMeaningfulDigits(std::cosh(c), c));
 }
 
-
 template Complex<float> Poincare::HyperbolicCosineNode::computeOnComplex<float>(std::complex<float>, Preferences::ComplexFormat, Preferences::AngleUnit);
 template Complex<double> Poincare::HyperbolicCosineNode::computeOnComplex<double>(std::complex<double>, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit);
 
diff --git a/poincare/src/hyperbolic_sine.cpp b/poincare/src/hyperbolic_sine.cpp
index d30b18b82..232bd8fcb 100644
--- a/poincare/src/hyperbolic_sine.cpp
+++ b/poincare/src/hyperbolic_sine.cpp
@@ -9,6 +9,7 @@ constexpr Expression::FunctionHelper HyperbolicSine::s_functionHelper;
 Layout HyperbolicSineNode::createLayout(Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const {
   return LayoutHelper::Prefix(HyperbolicSine(this), floatDisplayMode, numberOfSignificantDigits, HyperbolicSine::s_functionHelper.name());
 }
+
 int HyperbolicSineNode::serialize(char * buffer, int bufferSize, Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const {
   return SerializationHelper::Prefix(this, buffer, bufferSize, floatDisplayMode, numberOfSignificantDigits, HyperbolicSine::s_functionHelper.name());
 }
@@ -18,7 +19,6 @@ Complex<T> HyperbolicSineNode::computeOnComplex(const std::complex<T> c, Prefere
   return Complex<T>::Builder(Trigonometry::RoundToMeaningfulDigits(std::sinh(c), c));
 }
 
-
 template Complex<float> Poincare::HyperbolicSineNode::computeOnComplex<float>(std::complex<float>, Preferences::ComplexFormat, Preferences::AngleUnit);
 template Complex<double> Poincare::HyperbolicSineNode::computeOnComplex<double>(std::complex<double>, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit);
 
diff --git a/poincare/src/hyperbolic_tangent.cpp b/poincare/src/hyperbolic_tangent.cpp
index df19fe5e3..7c4ad545d 100644
--- a/poincare/src/hyperbolic_tangent.cpp
+++ b/poincare/src/hyperbolic_tangent.cpp
@@ -9,6 +9,7 @@ constexpr Expression::FunctionHelper HyperbolicTangent::s_functionHelper;
 Layout HyperbolicTangentNode::createLayout(Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const {
   return LayoutHelper::Prefix(HyperbolicTangent(this), floatDisplayMode, numberOfSignificantDigits, HyperbolicTangent::s_functionHelper.name());
 }
+
 int HyperbolicTangentNode::serialize(char * buffer, int bufferSize, Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const {
   return SerializationHelper::Prefix(this, buffer, bufferSize, floatDisplayMode, numberOfSignificantDigits, HyperbolicTangent::s_functionHelper.name());
 }
@@ -18,7 +19,6 @@ Complex<T> HyperbolicTangentNode::computeOnComplex(const std::complex<T> c, Pref
   return Complex<T>::Builder(Trigonometry::RoundToMeaningfulDigits(std::tanh(c), c));
 }
 
-
 template Complex<float> Poincare::HyperbolicTangentNode::computeOnComplex<float>(std::complex<float>, Preferences::ComplexFormat, Preferences::AngleUnit);
 template Complex<double> Poincare::HyperbolicTangentNode::computeOnComplex<double>(std::complex<double>, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit);
 
diff --git a/poincare/src/hyperbolic_trigonometric_function.cpp b/poincare/src/hyperbolic_trigonometric_function.cpp
index 11f779c30..ea464f6a8 100644
--- a/poincare/src/hyperbolic_trigonometric_function.cpp
+++ b/poincare/src/hyperbolic_trigonometric_function.cpp
@@ -1,24 +1,23 @@
 #include <poincare/hyperbolic_trigonometric_function.h>
+#include <poincare/simplification_helper.h>
 
 namespace Poincare {
 
 Expression HyperbolicTrigonometricFunctionNode::shallowReduce(Context & context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, ReductionTarget target, bool symbolicComputation) {
-  return HyperbolicTrigonometricFunction(this).shallowReduce();
+  return HyperbolicTrigonometricFunction(this).shallowReduce(context, angleUnit);
 }
 
-Expression HyperbolicTrigonometricFunction::shallowReduce() {
+Expression HyperbolicTrigonometricFunction::shallowReduce(Context & context, Preferences::AngleUnit angleUnit) {
   {
     Expression e = Expression::defaultShallowReduce();
     if (e.isUndefined()) {
       return e;
     }
   }
-#if MATRIX_EXACT_REDUCING
   Expression c = childAtIndex(0);
-  if (c.type() == ExpressionNode::Type::Matrix) {
+  if (childAtIndex(0).type() == ExpressionNode::Type::Matrix) {
     return SimplificationHelper::Map(*this, context, angleUnit);
   }
-#endif
   return *this;
 }
 
diff --git a/poincare/src/imaginary_part.cpp b/poincare/src/imaginary_part.cpp
index e9975fa91..6886fb054 100644
--- a/poincare/src/imaginary_part.cpp
+++ b/poincare/src/imaginary_part.cpp
@@ -32,11 +32,9 @@ Expression ImaginaryPart::shallowReduce(Context & context, Preferences::ComplexF
     }
   }
   Expression c = childAtIndex(0);
-#if MATRIX_EXACT_REDUCING
   if (c.type() == ExpressionNode::Type::Matrix) {
     return SimplificationHelper::Map(*this, context, angleUnit);
   }
-#endif
   if (c.isReal(context)) {
     Expression result = Rational::Builder(0);
     replaceWithInPlace(result);
diff --git a/poincare/src/integral.cpp b/poincare/src/integral.cpp
index d91a5e72f..6f5523a61 100644
--- a/poincare/src/integral.cpp
+++ b/poincare/src/integral.cpp
@@ -40,7 +40,7 @@ int IntegralNode::serialize(char * buffer, int bufferSize, Preferences::PrintFlo
 }
 
 Expression IntegralNode::shallowReduce(Context & context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, ReductionTarget target, bool symbolicComputation) {
-  return Integral(this).shallowReduce();
+  return Integral(this).shallowReduce(context);
 }
 
 template<typename T>
@@ -213,22 +213,20 @@ Expression Integral::UntypedBuilder(Expression children) {
   return Builder(children.childAtIndex(0), children.childAtIndex(1).convert<Symbol>(), children.childAtIndex(2), children.childAtIndex(3));
 }
 
-Expression Integral::shallowReduce() {
+Expression Integral::shallowReduce(Context & context) {
   {
     Expression e = Expression::defaultShallowReduce();
     if (e.isUndefined()) {
       return e;
     }
   }
-#if MATRIX_EXACT_REDUCING
-  if (childAtIndex(0).type() == ExpressionNode::Type::Matrix
-      || childAtIndex(1).type() == ExpressionNode::Type::Matrix
-      || childAtIndex(2).type() == ExpressionNode::Type::Matrix
-      || childAtIndex(3).type() == ExpressionNode::Type::Matrix)
+  if (SortedIsMatrix(childAtIndex(0), context)
+      || SortedIsMatrix(childAtIndex(1), context)
+      || SortedIsMatrix(childAtIndex(2), context)
+      || SortedIsMatrix(childAtIndex(3), context))
   {
     return Undefined::Builder();
   }
-#endif
   return *this;
 }
 
diff --git a/poincare/src/logarithm.cpp b/poincare/src/logarithm.cpp
index ea281ebf8..079d767a5 100644
--- a/poincare/src/logarithm.cpp
+++ b/poincare/src/logarithm.cpp
@@ -94,13 +94,9 @@ Expression CommonLogarithm::shallowReduce(Context & context, Preferences::Comple
     }
   }
   Expression c = childAtIndex(0);
-#if MATRIX_EXACT_REDUCING
-#if 0
   if (c.type() == ExpressionNode::Type::Matrix) {
-    return SimplificationHelper::Map(this, context, angleUnit);
+    return SimplificationHelper::Map(*this, context, angleUnit);
   }
-#endif
-#endif
   Logarithm log = Logarithm::Builder(childAtIndex(0), Rational::Builder(10));
   replaceWithInPlace(log);
   return log.shallowReduce(context, complexFormat, angleUnit, target);
@@ -113,14 +109,12 @@ Expression Logarithm::shallowReduce(Context & context, Preferences::ComplexForma
       return e;
     }
   }
-  Expression c = childAtIndex(0);
-#if MATRIX_EXACT_REDUCING
-#if 0
-  if (c.type() == ExpressionNode::Type::Matrix || childAtIndex(1).type() == ExpressionNode::Type::Matrix) {
+
+  if (SortedIsMatrix(childAtIndex(1), context)) {
     return Undefined::Builder();
   }
-#endif
-#endif
+
+  Expression c = childAtIndex(0);
   if (c.sign(&context) == ExpressionNode::Sign::Negative || childAtIndex(1).sign(&context) == ExpressionNode::Sign::Negative) {
     return *this;
   }
@@ -136,9 +130,12 @@ Expression Logarithm::shallowReduce(Context & context, Preferences::ComplexForma
    *   information on the parent which could later be a power of b.
    */
   bool letLogAtRoot = parentIsAPowerOfSameBase();
+  if (letLogAtRoot) {
+    return *this;
+  }
 
   // log(+inf, a) ?
-  if (!letLogAtRoot && c.type() == ExpressionNode::Type::Infinity && c.sign(&context) == ExpressionNode::Sign::Positive) {
+  if (c.type() == ExpressionNode::Type::Infinity && c.sign(&context) == ExpressionNode::Sign::Positive) {
     Expression base = childAtIndex(1);
     // log(+inf, a) --> ±inf with a rational and a > 0
     if (base.type() == ExpressionNode::Type::Rational && !static_cast<Rational&>(base).isNegative() && !static_cast<Rational&>(base).isZero()) {
@@ -156,7 +153,7 @@ Expression Logarithm::shallowReduce(Context & context, Preferences::ComplexForma
   }
 
   // log(x^y, b)->y*log(x, b) if x>0
-  if (!letLogAtRoot && c.type() == ExpressionNode::Type::Power && c.childAtIndex(0).sign(&context) == ExpressionNode::Sign::Positive) {
+  if (c.type() == ExpressionNode::Type::Power && c.childAtIndex(0).sign(&context) == ExpressionNode::Sign::Positive) {
     Power p = static_cast<Power &>(c);
     Expression x = p.childAtIndex(0);
     Expression y = p.childAtIndex(1);
@@ -168,7 +165,7 @@ Expression Logarithm::shallowReduce(Context & context, Preferences::ComplexForma
     return mult.shallowReduce(context, complexFormat, angleUnit, target);
   }
   // log(x*y, b)->log(x,b)+log(y, b) if x,y>0
-  if (!letLogAtRoot && c.type() == ExpressionNode::Type::Multiplication) {
+  if (c.type() == ExpressionNode::Type::Multiplication) {
     Addition a = Addition::Builder();
     for (int i = 0; i < c.numberOfChildren()-1; i++) {
       Expression factor = c.childAtIndex(i);
@@ -189,7 +186,7 @@ Expression Logarithm::shallowReduce(Context & context, Preferences::ComplexForma
     }
   }
   // log(r) with r Rational
-  if (!letLogAtRoot && c.type() == ExpressionNode::Type::Rational) {
+  if (c.type() == ExpressionNode::Type::Rational) {
     Rational r = static_cast<Rational &>(c);
     Addition a = Addition::Builder();
     // if the log base is Integer: log_b(r) = c + log_b(r') with r = b^c*r'
@@ -205,6 +202,12 @@ Expression Logarithm::shallowReduce(Context & context, Preferences::ComplexForma
     replaceWithInPlace(a);
     return a.shallowReduce(context, complexFormat, angleUnit, target);
   }
+
+  // log(m) with m Matrix
+  if (c.type() == ExpressionNode::Type::Matrix) {
+    return SimplificationHelper::Map(*this, context, angleUnit);
+  }
+
   return *this;
 }
 
diff --git a/poincare/src/n_ary_expression.cpp b/poincare/src/n_ary_expression.cpp
index e1d5b4c14..d2c328b65 100644
--- a/poincare/src/n_ary_expression.cpp
+++ b/poincare/src/n_ary_expression.cpp
@@ -95,4 +95,13 @@ int NAryExpression::allChildrenAreReal(Context & context) const {
   return result;
 }
 
+bool NAryExpression::SortedIsMatrix(Expression e, Context & context) {
+  assert(IsNAry(e, context));
+  int childrenCount = e.numberOfChildren();
+  if (childrenCount > 0) {
+     return SortedIsMatrix(e.childAtIndex(childrenCount - 1), context);
+  }
+  return false;
+}
+
 }
diff --git a/poincare/src/naperian_logarithm.cpp b/poincare/src/naperian_logarithm.cpp
index 5f62307dd..5d359fb28 100644
--- a/poincare/src/naperian_logarithm.cpp
+++ b/poincare/src/naperian_logarithm.cpp
@@ -30,12 +30,11 @@ Expression NaperianLogarithm::shallowReduce(Context & context, Preferences::Comp
       return e;
     }
   }
-#if MATRIX_EXACT_REDUCING
-  if (childAtIndex(0).type() == ExpressionNode::Type::Matrix) {
+  Expression c = childAtIndex(0);
+  if (c.type() == ExpressionNode::Type::Matrix) {
     return SimplificationHelper::Map(*this, context, angleUnit);
   }
-#endif
-  Logarithm l = Logarithm::Builder(childAtIndex(0), Constant::Builder(UCodePointScriptSmallE));
+  Logarithm l = Logarithm::Builder(c, Constant::Builder(UCodePointScriptSmallE));
   replaceWithInPlace(l);
   return l.shallowReduce(context, complexFormat, angleUnit, target);
 }
diff --git a/poincare/src/opposite.cpp b/poincare/src/opposite.cpp
index 70348262d..6e5361593 100644
--- a/poincare/src/opposite.cpp
+++ b/poincare/src/opposite.cpp
@@ -79,8 +79,6 @@ Expression Opposite::shallowReduce(Context & context, Preferences::ComplexFormat
     return result;
   }
   Expression child = result.childAtIndex(0);
-#if MATRIX_EXACT_REDUCING
-#endif
   result = Multiplication::Builder(Rational::Builder(-1), child);
   replaceWithInPlace(result);
   return result.shallowReduce(context, complexFormat, angleUnit, target);
diff --git a/poincare/src/real_part.cpp b/poincare/src/real_part.cpp
index 219fa6f87..d105b607f 100644
--- a/poincare/src/real_part.cpp
+++ b/poincare/src/real_part.cpp
@@ -33,11 +33,9 @@ Expression RealPart::shallowReduce(Context & context, Preferences::ComplexFormat
     }
   }
   Expression c = childAtIndex(0);
-#if MATRIX_EXACT_REDUCING
   if (c.type() == ExpressionNode::Type::Matrix) {
     return SimplificationHelper::Map(*this, context, angleUnit);
   }
-#endif
   if (c.isReal(context)) {
     replaceWithInPlace(c);
     return c;
diff --git a/poincare/src/sign_function.cpp b/poincare/src/sign_function.cpp
index e24218e72..1b448d438 100644
--- a/poincare/src/sign_function.cpp
+++ b/poincare/src/sign_function.cpp
@@ -3,6 +3,7 @@
 #include <poincare/rational.h>
 #include <poincare/layout_helper.h>
 #include <poincare/serialization_helper.h>
+#include <poincare/simplification_helper.h>
 #include <ion.h>
 #include <assert.h>
 #include <math.h>
@@ -59,13 +60,11 @@ Expression SignFunction::shallowReduce(Context & context, Preferences::ComplexFo
       return e;
     }
   }
-#if MATRIX_EXACT_REDUCING
-  if (c.type() == ExpressionNode::Type::Matrix) {
+  Expression child = childAtIndex(0);
+  if (child.type() == ExpressionNode::Type::Matrix) {
     return SimplificationHelper::Map(*this, context, angleUnit);
   }
-#endif
   Rational resultSign = Rational::Builder(1);
-  Expression child = childAtIndex(0);
   ExpressionNode::Sign s = child.sign(&context);
   if (s == ExpressionNode::Sign::Negative) {
     resultSign = Rational::Builder(-1);
diff --git a/poincare/src/sine.cpp b/poincare/src/sine.cpp
index 2d38addfe..7f94392d7 100644
--- a/poincare/src/sine.cpp
+++ b/poincare/src/sine.cpp
@@ -42,12 +42,9 @@ Expression Sine::shallowReduce(Context & context, Preferences::ComplexFormat com
       return e;
     }
   }
-#if MATRIX_EXACT_REDUCING
-  Expression op = childAtIndex(0);
-  if (op.type() == ExpressionNode::Type::Matrix) {
+  if (childAtIndex(0).type() == ExpressionNode::Type::Matrix) {
     return SimplificationHelper::Map(*this, context, angleUnit);
   }
-#endif
   return Trigonometry::shallowReduceDirectFunction(*this, context, complexFormat, angleUnit, target);
 }
 
diff --git a/poincare/src/tangent.cpp b/poincare/src/tangent.cpp
index 9c6b6186b..3923aef14 100644
--- a/poincare/src/tangent.cpp
+++ b/poincare/src/tangent.cpp
@@ -46,8 +46,7 @@ Expression Tangent::shallowReduce(Context & context, Preferences::ComplexFormat
     }
   }
 
-  Expression op = childAtIndex(0);
-  if (op.type() == ExpressionNode::Type::Matrix) {
+  if (childAtIndex(0).type() == ExpressionNode::Type::Matrix) {
     return SimplificationHelper::Map(*this, context, angleUnit);
   }