[poincare] Replace Poincare::Complex by std::complex

2026-01-19 00:37:25 +01:00 · 2018-07-05 13:29:55 +02:00
parent e7619633bd
commit dd0bc491f6
154 changed files with 1122 additions and 1558 deletions
--- a/apps/calculation/calculation.cpp
+++ b/apps/calculation/calculation.cpp
@@ -172,7 +172,7 @@ Expression * Calculation::approximateOutput(Context * context) {
      m_approximateOutput = exp->approximate<double>(*context);
      delete exp;
    } else {
-      m_approximateOutput = new Complex<double>(Complex<double>::Float(NAN));
+      m_approximateOutput = new Undefined();
    }
  }
  return m_approximateOutput;
--- a/apps/graph/cartesian_function.cpp
+++ b/apps/graph/cartesian_function.cpp
@@ -21,9 +21,9 @@ void CartesianFunction::setDisplayDerivative(bool display) {
 }

 double CartesianFunction::approximateDerivative(double x, Poincare::Context * context) const {
-  Poincare::Complex<double> abscissa = Poincare::Complex<double>::Float(x);
-  Poincare::Expression * args[2] = {expression(context), &abscissa};
-  Poincare::Derivative derivative(args, true);
+  Poincare::Expression * abscissa = Expression::CreateDecimal(x);
+  Poincare::Expression * args[2] = {expression(context)->clone(), abscissa};
+  Poincare::Derivative derivative(args, false); // derivative takes ownership of abscissa and the clone of expression
  /* TODO: when we will simplify derivative, we might want to simplify the
   * derivative here. However, we might want to do it once for all x (to avoid
   * lagging in the derivative table. */
@@ -31,10 +31,10 @@ double CartesianFunction::approximateDerivative(double x, Poincare::Context * co
 }

 double CartesianFunction::sumBetweenBounds(double start, double end, Poincare::Context * context) const {
-  Poincare::Complex<double> x = Poincare::Complex<double>::Float(start);
-  Poincare::Complex<double> y = Poincare::Complex<double>::Float(end);
-  Poincare::Expression * args[3] = {expression(context), &x, &y};
-  Poincare::Integral integral(args, true);
+  Poincare::Expression * x = Expression::CreateDecimal(start);
+  Poincare::Expression * y = Expression::CreateDecimal(end);
+  Poincare::Expression * args[3] = {expression(context)->clone(), x, y};
+  Poincare::Integral integral(args, false); // Integral takes ownership of args
  /* TODO: when we will simplify integral, we might want to simplify the
   * integral here. However, we might want to do it once for all x (to avoid
   * lagging in the derivative table. */
--- a/apps/regression/model/model.cpp
+++ b/apps/regression/model/model.cpp
@@ -1,6 +1,5 @@
 #include "model.h"
 #include "../store.h"
-#include <poincare/complex.h>
 #include <poincare/decimal.h>
 #include <poincare/matrix.h>
 #include <poincare/multiplication.h>
--- a/apps/regression/regression_context.cpp
+++ b/apps/regression/regression_context.cpp
@@ -20,8 +20,11 @@ const Expression * RegressionContext::expressionForSymbol(const Symbol * symbol)

    assert(m_seriesPairIndex >= 0);
    assert(m_seriesPairIndex < m_store->numberOfPairsOfSeries(series));
-    m_value = Complex<double>::Float(m_store->get(series, storeI, m_seriesPairIndex));
-    return &m_value;
+    //m_value = Complex<double>::Float(m_store->get(series, storeI, m_seriesPairIndex));
+    // FIXME
+    //m_value = Complex<double>::Undefined();
+    //return &m_value;
+    return nullptr;
  } else {
    return m_parentContext->expressionForSymbol(symbol);
  }
--- a/apps/sequence/cache_context.cpp
+++ b/apps/sequence/cache_context.cpp
@@ -5,30 +5,42 @@ using namespace Poincare;

 namespace Sequence {

-template<typename T>
-CacheContext<T>::CacheContext(Context * parentContext) :
-  VariableContext<T>('n', parentContext),
-  m_values{{Complex<T>::Float(NAN), Complex<T>::Float(NAN)},
-    {Complex<T>::Float(NAN), Complex<T>::Float(NAN)}}
+CacheContext::CacheContext(Context * parentContext) :
+  VariableContext('n', parentContext),
+  m_values{{new Undefined(), new Undefined()},
+    {new Undefined(), new Undefined()}}
 {
 }

-template<typename T>
-const Expression * CacheContext<T>::expressionForSymbol(const Symbol * symbol) {
+CacheContext::~CacheContext() {
+  Poincare::Expression * m_values[MaxNumberOfSequences][MaxRecurrenceDepth];
+  for (int i = 0; i < MaxNumberOfSequences; i++) {
+    for (int j = 0; j < MaxRecurrenceDepth; j++) {
+      if (m_values[i][j]) {
+        delete m_values[i][j];
+      }
+    }
+  }
+}
+
+const Expression * CacheContext::expressionForSymbol(const Symbol * symbol) {
  if (symbol->name() == Symbol::SpecialSymbols::un || symbol->name() == Symbol::SpecialSymbols::un1 ||
      symbol->name() == Symbol::SpecialSymbols::vn || symbol->name() == Symbol::SpecialSymbols::vn1) {
-    return &m_values[nameIndexForSymbol(symbol)][rankIndexForSymbol(symbol)];
+    return m_values[nameIndexForSymbol(symbol)][rankIndexForSymbol(symbol)];
  }
-  return VariableContext<T>::expressionForSymbol(symbol);
+  return VariableContext::expressionForSymbol(symbol);
 }

 template<typename T>
-void CacheContext<T>::setValueForSymbol(T value, const Poincare::Symbol * symbol) {
-  m_values[nameIndexForSymbol(symbol)][rankIndexForSymbol(symbol)] = Complex<T>::Float(value);
+void CacheContext::setValueForSymbol(T value, const Poincare::Symbol * symbol) {
+  if (m_values[nameIndexForSymbol(symbol)][rankIndexForSymbol(symbol)]) {
+    delete m_values[nameIndexForSymbol(symbol)][rankIndexForSymbol(symbol)];
+    m_values[nameIndexForSymbol(symbol)][rankIndexForSymbol(symbol)] = nullptr;
+  }
+  m_values[nameIndexForSymbol(symbol)][rankIndexForSymbol(symbol)] = Expression::CreateDecimal(value);
 }

-template<typename T>
-int CacheContext<T>::nameIndexForSymbol(const Poincare::Symbol * symbol) {
+int CacheContext::nameIndexForSymbol(const Poincare::Symbol * symbol) {
  switch (symbol->name()) {
    case Symbol::SpecialSymbols::un:
      return 0;
@@ -43,8 +55,7 @@ int CacheContext<T>::nameIndexForSymbol(const Poincare::Symbol * symbol) {
  }
 }

-template<typename T>
-int CacheContext<T>::rankIndexForSymbol(const Poincare::Symbol * symbol) {
+int CacheContext::rankIndexForSymbol(const Poincare::Symbol * symbol) {
  switch (symbol->name()) {
    case Symbol::SpecialSymbols::un:
      return 0;
@@ -59,7 +70,7 @@ int CacheContext<T>::rankIndexForSymbol(const Poincare::Symbol * symbol) {
  }
 }

-template class CacheContext<float>;
-template class CacheContext<double>;
+template void CacheContext::setValueForSymbol<double>(double, Poincare::Symbol const*);
+template void CacheContext::setValueForSymbol<float>(float, Poincare::Symbol const*);

 }
--- a/apps/sequence/cache_context.h
+++ b/apps/sequence/cache_context.h
@@ -6,16 +6,16 @@

 namespace Sequence {

-template<typename T>
-class CacheContext : public Poincare::VariableContext<T> {
+class CacheContext : public Poincare::VariableContext {
 public:
  CacheContext(Poincare::Context * parentContext);
+  ~CacheContext();
  const Poincare::Expression * expressionForSymbol(const Poincare::Symbol * symbol) override;
-  void setValueForSymbol(T value, const Poincare::Symbol * symbol);
+  template<typename T> void setValueForSymbol(T value, const Poincare::Symbol * symbol);
 private:
  int nameIndexForSymbol(const Poincare::Symbol * symbol);
  int rankIndexForSymbol(const Poincare::Symbol * symbol);
-  Poincare::Complex<T> m_values[MaxNumberOfSequences][MaxRecurrenceDepth];
+  Poincare::Expression * m_values[MaxNumberOfSequences][MaxRecurrenceDepth];
 };

 }
--- a/apps/sequence/sequence.cpp
+++ b/apps/sequence/sequence.cpp
@@ -319,14 +319,13 @@ T Sequence::approximateToNextRank(int n, SequenceContext * sqctx) const {
  if (n < m_initialRank || n < 0) {
    return NAN;
  }
-  CacheContext<T> ctx = CacheContext<T>(sqctx);
+  CacheContext ctx = CacheContext(sqctx);
  T un = sqctx->valueOfSequenceAtPreviousRank<T>(0, 0);
  T unm1 = sqctx->valueOfSequenceAtPreviousRank<T>(0, 1);
  T unm2 = sqctx->valueOfSequenceAtPreviousRank<T>(0, 2);
  T vn = sqctx->valueOfSequenceAtPreviousRank<T>(1, 0);
  T vnm1 = sqctx->valueOfSequenceAtPreviousRank<T>(1, 1);
  T vnm2 = sqctx->valueOfSequenceAtPreviousRank<T>(1, 2);
-  Poincare::Symbol nSymbol(symbol());
  Poincare::Symbol vnSymbol(Symbol::SpecialSymbols::vn);
  Poincare::Symbol vn1Symbol(Symbol::SpecialSymbols::vn1);
  Poincare::Symbol unSymbol(Symbol::SpecialSymbols::un);
--- a/apps/shared/store_context.h
+++ b/apps/shared/store_context.h
@@ -2,7 +2,6 @@
 #define SHARED_STORE_CONTEXT_H

 #include <poincare/context.h>
-#include <poincare/complex.h>
 #include "double_pair_store.h"
 #include <cmath>

@@ -14,8 +13,8 @@ public:
    Poincare::Context(),
    m_store(store),
    m_seriesPairIndex(-1),
-    m_parentContext(nullptr),
-    m_value(Poincare::Complex<double>::Float(NAN))
+    m_parentContext(nullptr)//,
+    //m_value(Poincare::Complex<double>::Float(NAN))
  {}
  void setParentContext(Poincare::Context * parentContext) { m_parentContext = parentContext; }
  void setSeriesPairIndex(int j) { m_seriesPairIndex = j; }
@@ -24,7 +23,7 @@ protected:
  Shared::DoublePairStore * m_store;
  int m_seriesPairIndex;
  Poincare::Context * m_parentContext;
-  Poincare::Complex<double> m_value;
+  //Poincare::Complex<double> m_value; // FIXME
 };

 }
--- a/apps/statistics/statistics_context.cpp
+++ b/apps/statistics/statistics_context.cpp
@@ -20,8 +20,11 @@ const Expression * StatisticsContext::expressionForSymbol(const Symbol * symbol)

    assert(m_seriesPairIndex >= 0);
    assert(m_seriesPairIndex < m_store->numberOfPairsOfSeries(series));
-    m_value = Complex<double>::Float(m_store->get(series, storeI, m_seriesPairIndex));
-    return &m_value;
+    //m_value = Complex<double>::Float(m_store->get(series, storeI, m_seriesPairIndex));
+    // FIXME
+    //m_value = Complex<double>::Undefined();
+    //return &m_value;
+    return nullptr;
  } else {
    return m_parentContext->expressionForSymbol(symbol);
  }
--- a/poincare/Makefile
+++ b/poincare/Makefile
@@ -14,7 +14,6 @@ objs += $(addprefix poincare/src/,\
  binomial_coefficient.o\
  bounded_static_hierarchy.o\
  ceiling.o\
-  complex.o\
  complex_argument.o\
  confidence_interval.o\
  conjugate.o\
@@ -29,6 +28,7 @@ objs += $(addprefix poincare/src/,\
  empty_expression.o\
  equal.o\
  expression_layout_cursor.o\
+  evaluation.o\
  expression_lexer.o\
  expression_parser.o\
  expression.o\
--- a/poincare/include/poincare.h
+++ b/poincare/include/poincare.h
@@ -10,7 +10,6 @@
 #include <poincare/arc_tangent.h>
 #include <poincare/binomial_coefficient.h>
 #include <poincare/ceiling.h>
-#include <poincare/complex.h>
 #include <poincare/complex_argument.h>
 #include <poincare/confidence_interval.h>
 #include <poincare/conjugate.h>
--- a/poincare/include/poincare/absolute_value.h
+++ b/poincare/include/poincare/absolute_value.h
@@ -23,11 +23,11 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  template<typename T> static Complex<T> computeOnComplex(const Complex<T> c, AngleUnit angleUnit);
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  template<typename T> static std::complex<T> computeOnComplex(const std::complex<T> c, AngleUnit angleUnit);
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<float>(this, context, angleUnit, computeOnComplex<float>);
  }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
  return ApproximationEngine::map<double>(this, context, angleUnit, computeOnComplex<double>);
  }
 };
--- a/poincare/include/poincare/addition.h
+++ b/poincare/include/poincare/addition.h
@@ -22,12 +22,12 @@ public:
  int polynomialDegree(char symbolName) const override;
  int privateGetPolynomialCoefficients(char symbolName, Expression * coefficients[]) const override;
  /* Evaluation */
-  template<typename T> static Complex<T> compute(const Complex<T> c, const Complex<T> d);
-  template<typename T> static Matrix * computeOnMatrices(const Matrix * m, const Matrix * n) {
+  template<typename T> static std::complex<T> compute(const std::complex<T> c, const std::complex<T> d);
+  template<typename T> static MatrixComplex<T> computeOnMatrices(const MatrixComplex<T> m, const MatrixComplex<T> n) {
    return ApproximationEngine::elementWiseOnComplexMatrices(m, n, compute<T>);
  }
-  template<typename T> static Matrix * computeOnComplexAndMatrix(const Complex<T> * c, const Matrix * m) {
-    return ApproximationEngine::elementWiseOnComplexAndComplexMatrix(c, m, compute<T>);
+  template<typename T> static MatrixComplex<T> computeOnComplexAndMatrix(const std::complex<T> c, const MatrixComplex<T> m) {
+    return ApproximationEngine::elementWiseOnMatrixComplexAndComplex(m, c, compute<T>);
  }
 private:
  /* Layout */
@@ -48,13 +48,13 @@ private:
  static const Rational RationalFactor(Expression * e);
  static bool TermsHaveIdenticalNonRationalFactors(const Expression * e1, const Expression * e2);
  /* Evaluation */
-  template<typename T> static Matrix * computeOnMatrixAndComplex(const Matrix * m, const Complex<T> * c) {
-    return ApproximationEngine::elementWiseOnComplexAndComplexMatrix(c, m, compute<T>);
+  template<typename T> static MatrixComplex<T> computeOnMatrixAndComplex(const MatrixComplex<T> m, const std::complex<T> c) {
+    return ApproximationEngine::elementWiseOnMatrixComplexAndComplex(m, c, compute<T>);
  }
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::mapReduce<float>(this, context, angleUnit, compute<float>, computeOnComplexAndMatrix<float>, computeOnMatrixAndComplex<float>, computeOnMatrices<float>);
   }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::mapReduce<double>(this, context, angleUnit, compute<double>, computeOnComplexAndMatrix<double>, computeOnMatrixAndComplex<double>, computeOnMatrices<double>);
   }
 };
--- a/poincare/include/poincare/approximation_engine.h
+++ b/poincare/include/poincare/approximation_engine.h
@@ -1,25 +1,25 @@
 #ifndef POINCARE_APPROXIMATION_ENGINE_H
 #define POINCARE_APPROXIMATION_ENGINE_H

+#include <poincare/evaluation.h>
 #include <poincare/expression.h>
-#include <poincare/complex.h>
-#include <poincare/matrix.h>
+#include <complex>

 namespace Poincare {

 class ApproximationEngine {
 public:
-  template <typename T> using ComplexCompute = Complex<T>(*)(const Complex<T>, Expression::AngleUnit angleUnit);
-  template<typename T> static Expression * map(const Expression * expression, Context& context, Expression::AngleUnit angleUnit, ComplexCompute<T> compute);
+  template <typename T> using ComplexCompute = std::complex<T>(*)(const std::complex<T>, Expression::AngleUnit angleUnit);
+  template<typename T> static Evaluation<T> * map(const Expression * expression, Context& context, Expression::AngleUnit angleUnit, ComplexCompute<T> compute);

-  template <typename T> using ComplexAndComplexReduction = Complex<T>(*)(const Complex<T>, const Complex<T>);
-  template <typename T> using ComplexAndMatrixReduction = Matrix * (*)(const Complex<T> * c, const Matrix * m);
-  template <typename T> using MatrixAndComplexReduction = Matrix * (*)(const Matrix * m, const Complex<T> * c);
-  template <typename T> using MatrixAndMatrixReduction = Matrix * (*)(const Matrix * m, const Matrix * n);
-  template<typename T> static Expression * mapReduce(const Expression * expression, Context& context, Expression::AngleUnit angleUnit, ComplexAndComplexReduction<T> computeOnComplexes, ComplexAndMatrixReduction<T> computeOnComplexAndMatrix, MatrixAndComplexReduction<T> computeOnMatrixAndComplex, MatrixAndMatrixReduction<T> computeOnMatrices);
+  template <typename T> using ComplexAndComplexReduction = std::complex<T>(*)(const std::complex<T>, const std::complex<T>);
+  template <typename T> using ComplexAndMatrixReduction = MatrixComplex<T>(*)(const std::complex<T> c, const MatrixComplex<T> m);
+  template <typename T> using MatrixAndComplexReduction = MatrixComplex<T>(*)(const MatrixComplex<T> m, const std::complex<T> c);
+  template <typename T> using MatrixAndMatrixReduction = MatrixComplex<T>(*)(const MatrixComplex<T> m, const MatrixComplex<T> n);
+  template<typename T> static Evaluation<T> * mapReduce(const Expression * expression, Context& context, Expression::AngleUnit angleUnit, ComplexAndComplexReduction<T> computeOnComplexes, ComplexAndMatrixReduction<T> computeOnComplexAndMatrix, MatrixAndComplexReduction<T> computeOnMatrixAndComplex, MatrixAndMatrixReduction<T> computeOnMatrices);

-  template<typename T> static Matrix * elementWiseOnComplexAndComplexMatrix(const Complex<T> * c, const Matrix * n, ComplexAndComplexReduction<T> computeOnComplexes);
-  template<typename T> static Matrix * elementWiseOnComplexMatrices(const Matrix * m, const Matrix * n, ComplexAndComplexReduction<T> computeOnComplexes);
+  template<typename T> static MatrixComplex<T> elementWiseOnMatrixComplexAndComplex(const MatrixComplex<T> n, const std::complex<T> c, ComplexAndComplexReduction<T> computeOnComplexes);
+  template<typename T> static MatrixComplex<T> elementWiseOnComplexMatrices(const MatrixComplex<T> m, const MatrixComplex<T> n, ComplexAndComplexReduction<T> computeOnComplexes);

 };

--- a/poincare/include/poincare/arc_cosine.h
+++ b/poincare/include/poincare/arc_cosine.h
@@ -26,11 +26,11 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  template<typename T> static Complex<T> computeOnComplex(const Complex<T> c, AngleUnit angleUnit);
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  template<typename T> static std::complex<T> computeOnComplex(const std::complex<T> c, AngleUnit angleUnit);
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<float>(this, context, angleUnit,computeOnComplex<float>);
  }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<double>(this, context, angleUnit, computeOnComplex<double>);
  }
 };
--- a/poincare/include/poincare/arc_sine.h
+++ b/poincare/include/poincare/arc_sine.h
@@ -24,11 +24,11 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context & context, AngleUnit angleUnit) override;
  /* Evaluation */
-  template<typename T> static Complex<T> computeOnComplex(const Complex<T> c, AngleUnit angleUnit);
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
-    return ApproximationEngine::map<float>(this, context, angleUnit,computeOnComplex<float>);
+  template<typename T> static std::complex<T> computeOnComplex(const std::complex<T> c, AngleUnit angleUnit);
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+    return ApproximationEngine::map<float>(this, context, angleUnit, computeOnComplex<float>);
  }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<double>(this, context, angleUnit, computeOnComplex<double>);
  }
 };
--- a/poincare/include/poincare/arc_tangent.h
+++ b/poincare/include/poincare/arc_tangent.h
@@ -24,11 +24,11 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  template<typename T> static Complex<T> computeOnComplex(const Complex<T> c, AngleUnit angleUnit);
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  template<typename T> static std::complex<T> computeOnComplex(const std::complex<T> c, AngleUnit angleUnit);
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<float>(this, context, angleUnit,computeOnComplex<float>);
  }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<double>(this, context, angleUnit, computeOnComplex<double>);
  }
 };
--- a/poincare/include/poincare/binomial_coefficient.h
+++ b/poincare/include/poincare/binomial_coefficient.h
@@ -22,9 +22,9 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
-  template<typename T> Expression * templatedApproximate(Context& context, AngleUnit angleUnit) const;
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
+  template<typename T> Complex<T> * templatedApproximate(Context& context, AngleUnit angleUnit) const;
 };

 }
--- a/poincare/include/poincare/ceiling.h
+++ b/poincare/include/poincare/ceiling.h
@@ -22,11 +22,11 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  template<typename T> static Complex<T> computeOnComplex(const Complex<T> c, AngleUnit angleUnit);
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  template<typename T> static std::complex<T> computeOnComplex(const std::complex<T> c, AngleUnit angleUnit);
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<float>(this, context, angleUnit, computeOnComplex<float>);
  }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<double>(this, context, angleUnit, computeOnComplex<double>);
  }
 };
--- a/poincare/include/poincare/complex.h
+++ b/poincare/include/poincare/complex.h
@@ -1,65 +0,0 @@
-#ifndef POINCARE_COMPLEX_H
-#define POINCARE_COMPLEX_H
-
-#include <poincare/preferences.h>
-#include <poincare/static_hierarchy.h>
-#include <poincare/integer.h>
-#include <assert.h>
-
-namespace Poincare {
-
-template<typename T>
-class Complex : public StaticHierarchy<0> {
-public:
-  Complex() : m_a(0), m_b(0) {}
-  static Complex<T> Float(T x);
-  static Complex<T> Cartesian(T a, T b);
-  static Complex<T> Polar(T r, T theta);
-  Complex(const char * integralPart, int integralPartLength, bool integralNegative,
-        const char * fractionalPart, int fractionalPartLength,
-        const char * exponent, int exponentLength, bool exponentNegative);
-  Complex(const Complex & other);
-  Complex& operator=(const Complex& other);
-
-  T a() const;
-  T b() const;
-  T r() const;
-  T th() const;
-  Complex<T> conjugate() const;
-  T toScalar() const;
-
-  /* Expression */
-  Expression::Type type() const override;
-  Complex<T> * clone() const override;
-  int writeTextInBuffer(char * buffer, int bufferSize, int numberOfSignificantDigits = PrintFloat::k_numberOfStoredSignificantDigits) const override;
-
-  /* Properties */
-  bool needParenthesisWithParent(const Expression * e) const override;
-
-  /* Simplification: complex does not implement simplificationOrderSameType
-   * because Complex expressions are always transformed into an addition of
-   * Decimal and I symbol before compared with another Expression. */
-
-private:
-  Complex(T a, T b);
-  /* Layout */
-  ExpressionLayout * privateCreateLayout(PrintFloat::Mode floatDisplayMode, Expression::ComplexFormat complexFormat) const override;
-  /* Simplification */
-  static Expression * CreateDecimal(T f);
-  Expression * shallowReduce(Context & context, AngleUnit angleUnit) override;
-  /* Evaluation */
-  Expression * privateApproximate(Expression::SinglePrecision p, Context& context, Expression::AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
-  Expression * privateApproximate(Expression::DoublePrecision p, Context& context, Expression::AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
- template<typename U> Complex<U> * templatedApproximate(Context& context, Expression::AngleUnit angleUnit) const;
-  /* convertComplexToText and convertFloatToTextPrivate return the string length
-   * of the buffer (does not count the 0 last char)*/
-  int convertComplexToText(char * buffer, int bufferSize, int numberOfSignificantDigits, PrintFloat::Mode floatDisplayMode, Expression::ComplexFormat complexFormat, char multiplicationSign) const;
-  ExpressionLayout * createPolarLayout(PrintFloat::Mode floatDisplayMode) const;
-  ExpressionLayout * createCartesianLayout(PrintFloat::Mode floatDisplayMode) const;
-  T m_a;
-  T m_b;
-};
-
-}
-
-#endif
--- a/poincare/include/poincare/complex_argument.h
+++ b/poincare/include/poincare/complex_argument.h
@@ -24,11 +24,11 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  template<typename T> static Complex<T> computeOnComplex(const Complex<T> c, AngleUnit angleUnit);
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  template<typename T> static std::complex<T> computeOnComplex(const std::complex<T> c, AngleUnit angleUnit);
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<float>(this, context, angleUnit,computeOnComplex<float>);
  }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<double>(this, context, angleUnit, computeOnComplex<double>);
  }
 };
--- a/poincare/include/poincare/confidence_interval.h
+++ b/poincare/include/poincare/confidence_interval.h
@@ -24,9 +24,9 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  Expression * privateApproximate(Expression::SinglePrecision p, Context& context, Expression::AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
-  Expression * privateApproximate(Expression::DoublePrecision p, Context& context, Expression::AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
-  template<typename T> Expression * templatedApproximate(Context& context, AngleUnit angleUnit) const;
+  Evaluation<float> * privateApproximate(Expression::SinglePrecision p, Context& context, Expression::AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
+  Evaluation<double> * privateApproximate(Expression::DoublePrecision p, Context& context, Expression::AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
+  template<typename T> Evaluation<T> * templatedApproximate(Context& context, AngleUnit angleUnit) const;
 };

 class SimplePredictionInterval : public ConfidenceInterval {
--- a/poincare/include/poincare/conjugate.h
+++ b/poincare/include/poincare/conjugate.h
@@ -21,11 +21,11 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  template<typename T> static Complex<T> computeOnComplex(const Complex<T> c, AngleUnit angleUnit);
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  template<typename T> static std::complex<T> computeOnComplex(const std::complex<T> c, AngleUnit angleUnit);
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<float>(this, context, angleUnit,computeOnComplex<float>);
  }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<double>(this, context, angleUnit, computeOnComplex<double>);
  }
 };
--- a/poincare/include/poincare/cosine.h
+++ b/poincare/include/poincare/cosine.h
@@ -16,7 +16,9 @@ public:
  Type type() const override;
  Expression * clone() const override;
  float characteristicXRange(Context & context, AngleUnit angleUnit = AngleUnit::Default) const override;
-  template<typename T> static Complex<T> computeOnComplex(const Complex<T> c, AngleUnit angleUnit = AngleUnit::Radian);
+  template<typename T> static std::complex<T> computeOnComplex(const std::complex<T> c, AngleUnit angleUnit = AngleUnit::Radian) {
+    return Trigonometry::computeOnComplex(c, angleUnit, std::cos);
+  }
 private:
  /* Layout */
  ExpressionLayout * privateCreateLayout(PrintFloat::Mode floatDisplayMode, ComplexFormat complexFormat) const override {
@@ -29,10 +31,10 @@ private:
  /* Simplication */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<float>(this, context, angleUnit,computeOnComplex<float>);
  }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<double>(this, context, angleUnit, computeOnComplex<double>);
  }
 };
--- a/poincare/include/poincare/decimal.h
+++ b/poincare/include/poincare/decimal.h
@@ -36,9 +36,9 @@ private:
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  Expression * shallowBeautify(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
-  template<typename T> Expression * templatedApproximate(Context& context, Expression::AngleUnit angleUnit) const;
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
+  template<typename T> Evaluation<T> * templatedApproximate(Context& context, Expression::AngleUnit angleUnit) const;

  int convertToText(char * buffer, int bufferSize, PrintFloat::Mode mode, int numberOfSignificantDigits) const;
  // Worst case is -1.2345678901234E-1000
--- a/poincare/include/poincare/derivative.h
+++ b/poincare/include/poincare/derivative.h
@@ -25,9 +25,9 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
-  template<typename T> Expression * templatedApproximate(Context& context, AngleUnit angleUnit) const;
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
+  template<typename T> Complex<T> * templatedApproximate(Context& context, AngleUnit angleUnit) const;
  template<typename T> T growthRateAroundAbscissa(T x, T h, Context & context, AngleUnit angleUnit) const;
  template<typename T> T riddersApproximation(Context & context, AngleUnit angleUnit, T x, T h, T * error) const;
  // TODO: Change coefficients?
--- a/poincare/include/poincare/determinant.h
+++ b/poincare/include/poincare/determinant.h
@@ -3,7 +3,6 @@

 #include <poincare/layout_engine.h>
 #include <poincare/static_hierarchy.h>
-#include <poincare/complex.h>

 namespace Poincare {

@@ -24,9 +23,9 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context & context, AngleUnit angleUnit) override;
  /* Evaluation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
- template<typename T> Expression * templatedApproximate(Context& context, AngleUnit angleUnit) const;
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
+ template<typename T> Evaluation<T> * templatedApproximate(Context& context, AngleUnit angleUnit) const;
 };

 }
--- a/poincare/include/poincare/division.h
+++ b/poincare/include/poincare/division.h
@@ -15,7 +15,7 @@ class Division : public StaticHierarchy<2> {
 public:
  Type type() const override;
  Expression * clone() const override;
-  template<typename T> static Complex<T> compute(const Complex<T> c, const Complex<T> d);
+  template<typename T> static std::complex<T> compute(const std::complex<T> c, const std::complex<T> d);
  int polynomialDegree(char symbolName) const override;
 private:
  /* Layout */
@@ -27,16 +27,15 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  template<typename T> static Matrix * computeOnMatrixAndComplex(const Matrix * m, const Complex<T> * c) {
-    return ApproximationEngine::elementWiseOnComplexAndComplexMatrix(c, m, compute<T>);
+  template<typename T> static MatrixComplex<T> computeOnMatrixAndComplex(const MatrixComplex<T> m, const std::complex<T> c) {
+    return ApproximationEngine::elementWiseOnMatrixComplexAndComplex(m, c, compute<T>);
  }
-  template<typename T> static Matrix * computeOnComplexAndMatrix(const Complex<T> * c, const Matrix * n);
-  template<typename T> static Matrix * computeOnMatrices(const Matrix * m, const Matrix * n);
-
-  virtual Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  template<typename T> static MatrixComplex<T> computeOnComplexAndMatrix(const std::complex<T> c, const MatrixComplex<T> n);
+  template<typename T> static MatrixComplex<T> computeOnMatrices(const MatrixComplex<T> m, const MatrixComplex<T> n);
+  virtual Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::mapReduce<float>(this, context, angleUnit, compute<float>, computeOnComplexAndMatrix<float>, computeOnMatrixAndComplex<float>, computeOnMatrices<float>);
  }
-  virtual Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  virtual Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::mapReduce<double>(this, context, angleUnit, compute<double>, computeOnComplexAndMatrix<double>, computeOnMatrixAndComplex<double>, computeOnMatrices<double>);
  }
 };
--- a/poincare/include/poincare/division_quotient.h
+++ b/poincare/include/poincare/division_quotient.h
@@ -3,7 +3,6 @@

 #include <poincare/layout_engine.h>
 #include <poincare/static_hierarchy.h>
-#include <poincare/complex.h>

 namespace Poincare {

@@ -24,8 +23,8 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context & context, AngleUnit angleUnit) override;
  /* Evaluation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
 template<typename T> Complex<T> * templatedApproximate(Context& context, AngleUnit angleUnit) const;
 };

--- a/poincare/include/poincare/division_remainder.h
+++ b/poincare/include/poincare/division_remainder.h
@@ -3,7 +3,6 @@

 #include <poincare/layout_engine.h>
 #include <poincare/static_hierarchy.h>
-#include <poincare/complex.h>

 namespace Poincare {

@@ -24,8 +23,8 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context & context, AngleUnit angleUnit) override;
  /* Evaluation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
 template<typename T> Complex<T> * templatedApproximate(Context& context, AngleUnit angleUnit) const;

 };
--- a/poincare/include/poincare/empty_expression.h
+++ b/poincare/include/poincare/empty_expression.h
@@ -2,7 +2,7 @@
 #define POINCARE_EMPTY_EXPRESSION_H

 #include <poincare/static_hierarchy.h>
-#include <poincare/complex.h>
+#include <poincare/evaluation.h>

 namespace Poincare {

@@ -19,8 +19,8 @@ private:
  /* Layout */
  ExpressionLayout * privateCreateLayout(PrintFloat::Mode floatDisplayMode, ComplexFormat complexFormat) const override;
  /* Evaluation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
  template<typename T> Complex<T> * templatedApproximate(Context& context, AngleUnit angleUnit) const;
 };

--- a/poincare/include/poincare/equal.h
+++ b/poincare/include/poincare/equal.h
@@ -24,9 +24,9 @@ private:
    return LayoutEngine::writeInfixExpressionTextInBuffer(this, buffer, bufferSize, numberOfSignificantDigits, "=");
  }
  /* Evalutation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
-  template<typename T> Expression * templatedApproximate(Context& context, AngleUnit angleUnit) const;
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
+  template<typename T> Evaluation<T> * templatedApproximate(Context& context, AngleUnit angleUnit) const;
 };

 }
--- a/poincare/include/poincare/evaluation.h
+++ b/poincare/include/poincare/evaluation.h
@@ -0,0 +1,78 @@
+#ifndef POINCARE_EVALUATION_H
+#define POINCARE_EVALUATION_H
+
+#include <complex>
+extern "C" {
+#include <stdint.h>
+}
+
+namespace Poincare {
+
+template<typename T>
+class Evaluation {
+public:
+  enum class Type : uint8_t {
+    Complex,
+    MatrixComplex
+  };
+  virtual Type type() const = 0;
+  virtual ~Evaluation() {}
+  virtual bool isUndefined() const = 0;
+  virtual T toScalar() const { return NAN; }
+  virtual std::complex<T> createTrace() const = 0;
+  virtual std::complex<T> createDeterminant() const = 0;
+  virtual Evaluation * createInverse() const = 0;
+  virtual Evaluation * createTranspose() const = 0;
+};
+
+template<typename T>
+class Complex : public std::complex<T>, public Evaluation<T> {
+public:
+  Complex(T a, T b = 0.0) : std::complex<T>(a, b) {}
+  Complex(std::complex<T> c) : std::complex<T>(c) {}
+  static Complex Undefined() {
+    return Complex(NAN, NAN);
+  }
+  virtual ~Complex() {}
+  typename Poincare::Evaluation<T>::Type type() const override { return Poincare::Evaluation<T>::Type::Complex; }
+  bool isUndefined() const override {
+    return (std::isnan(this->real()) && std::isnan(this->imag()));
+  }
+  T toScalar() const override;
+  std::complex<T> createTrace() const override { return *this; }
+  std::complex<T> createDeterminant() const override { return *this; }
+  Complex<T> * createInverse() const override;
+  Complex<T> * createTranspose() const override { return new Complex<T>(*this); }
+};
+
+template<typename T>
+class MatrixComplex : public Evaluation<T> {
+public:
+  MatrixComplex(std::complex<T> * operands, int numberOfRows, int numberOfColumns);
+  static MatrixComplex Undefined() {
+    std::complex<T> undef = std::complex<T>(NAN, NAN);
+    return MatrixComplex<T>(&undef, 1, 1);
+  }
+  virtual ~MatrixComplex() {}
+  typename Poincare::Evaluation<T>::Type type() const override { return Poincare::Evaluation<T>::Type::MatrixComplex; }
+  const std::complex<T> complexOperand(int i) const { return m_operands[i]; }
+  int numberOfComplexOperands() const { return m_numberOfRows*m_numberOfColumns; }
+  int numberOfRows() const { return m_numberOfRows; }
+  int numberOfColumns() const { return m_numberOfColumns; }
+  bool isUndefined() const override {
+    return (numberOfRows() == 1 && numberOfColumns() == 1 && std::isnan(complexOperand(0).real()) && std::isnan(complexOperand(0).imag()));
+  }
+  std::complex<T> createTrace() const override;
+  std::complex<T> createDeterminant() const override;
+  MatrixComplex<T> * createInverse() const override;
+  MatrixComplex<T> * createTranspose() const override;
+  static MatrixComplex<T> createIdentity(int dim);
+private:
+  std::complex<T> * m_operands;
+  int m_numberOfRows;
+  int m_numberOfColumns;
+};
+
+}
+
+#endif
--- a/poincare/include/poincare/evaluation_mode.h
+++ b/poincare/include/poincare/evaluation_mode.h
@@ -0,0 +1,18 @@
+#ifndef POINCARE_EVALUATION_MODE_H
+#define POINCARE_EVALUATION_MODE_H
+
+namespace Poincare {
+
+enum class ComplexFormat {
+  Cartesian = 0,
+  Polar = 1,
+  Default = 2
+};
+
+enum class AngleUnit {
+  Degree = 0,
+  Radian = 1,
+  Default = 2
+};
+
+}
--- a/poincare/include/poincare/expression.h
+++ b/poincare/include/poincare/expression.h
@@ -3,6 +3,8 @@

 #include <poincare/expression_layout.h>
 #include <poincare/print_float.h>
+#include <poincare/evaluation.h>
+#include <complex>
 extern "C" {
 #include <assert.h>
 }
@@ -10,7 +12,6 @@ extern "C" {
 namespace Poincare {

 class Context;
-template<class T> class Complex;
 class Rational;

 class Expression {
@@ -78,10 +79,8 @@ class Expression {
  friend class ApproximationEngine;
  friend class SimplificationEngine;
  friend class LayoutEngine;
-  friend class Complex<float>;
-  friend class Complex<double>;
  friend class EmptyExpression;
-
+  friend class Randint;
 public:
  enum class Type : uint8_t {
    Undefined = 0,
@@ -139,7 +138,6 @@ public:
    Sum,
    Symbol,

-    Complex,
    Matrix,
    ConfidenceInterval,
    MatrixDimension,
@@ -257,12 +255,13 @@ public:
  static void Simplify(Expression ** expressionAddress, Context & context, AngleUnit angleUnit = AngleUnit::Default);
  static void Reduce(Expression ** expressionAddress, Context & context, AngleUnit angleUnit = AngleUnit::Default, bool recursively = true);

-  /* Evaluation Engine
-   * The function evaluate creates a new expression and thus mallocs memory.
+  /* Evaluation Engine */
+  template<typename T> static Expression * CreateDecimal(T f);
+  /* The function approximate creates a new expression and thus mallocs memory.
   * Do not forget to delete the new expression to avoid leaking. */
-  template<typename T> Expression * approximate(Context& context, AngleUnit angleUnit = AngleUnit::Default) const;
-  template<typename T> T approximateToScalar(Context& context, AngleUnit angleUnit = AngleUnit::Default) const;
-  template<typename T> static T approximateToScalar(const char * text, Context& context, AngleUnit angleUnit = AngleUnit::Default);
+  template<typename T> Expression * approximate(Context& context, AngleUnit angleUnit = AngleUnit::Default, ComplexFormat complexFormat = ComplexFormat::Default) const;
+  template<typename T> T approximateToScalar(Context& context, AngleUnit angleUnit = AngleUnit::Default, ComplexFormat complexFormat = ComplexFormat::Default) const;
+  template<typename T> static T approximateToScalar(const char * text, Context& context, AngleUnit angleUnit = AngleUnit::Default, ComplexFormat complexFormat = ComplexFormat::Default);
  template<typename T> T approximateWithValueForSymbol(char symbol, T x, Context & context) const;

  /* Expression roots/extrema solver*/
@@ -330,8 +329,9 @@ private:
    return nullptr;
  }
  /* Evaluation Engine */
-  virtual Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const = 0;
-  virtual Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const = 0;
+  template<typename T> static Expression * complexToExpression(std::complex<T> c, ComplexFormat complexFormat);
+  virtual Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const = 0;
+  virtual Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const = 0;

  /* Expression roots/extrema solver*/
  constexpr static double k_solverPrecision = 1.0E-5;
--- a/poincare/include/poincare/factor.h
+++ b/poincare/include/poincare/factor.h
@@ -5,7 +5,6 @@
 #include <poincare/static_hierarchy.h>
 #include <poincare/multiplication.h>
 #include <poincare/rational.h>
-#include <poincare/complex.h>
 #include <cmath>

 namespace Poincare {
@@ -28,10 +27,10 @@ private:
  Expression * shallowBeautify(Context& context, AngleUnit angleUnit) override;
  Expression * createMultiplicationOfIntegerPrimeDecomposition(Integer i, Context & context, AngleUnit angleUnit);
  /* Evaluation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
-  template<typename T> Expression * templatedApproximate(Context& context, AngleUnit angleUnit) const {
-    return operand(0)->approximate<T>(context, angleUnit);
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
+  template<typename T> Evaluation<T> * templatedApproximate(Context& context, AngleUnit angleUnit) const {
+    return operand(0)->privateApproximate(T(), context, angleUnit);
  }
 };

--- a/poincare/include/poincare/factorial.h
+++ b/poincare/include/poincare/factorial.h
@@ -21,11 +21,11 @@ private:
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  Expression * shallowBeautify(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  template<typename T> static Complex<T> computeOnComplex(const Complex<T> c, AngleUnit angleUnit);
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  template<typename T> static std::complex<T> computeOnComplex(const std::complex<T> c, AngleUnit angleUnit);
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<float>(this, context, angleUnit,computeOnComplex<float>);
  }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<double>(this, context, angleUnit, computeOnComplex<double>);
  }

--- a/poincare/include/poincare/floor.h
+++ b/poincare/include/poincare/floor.h
@@ -22,11 +22,11 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  template<typename T> static Complex<T> computeOnComplex(const Complex<T> c, AngleUnit angleUnit);
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  template<typename T> static std::complex<T> computeOnComplex(const std::complex<T> c, AngleUnit angleUnit);
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<float>(this, context, angleUnit, computeOnComplex<float>);
  }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<double>(this, context, angleUnit, computeOnComplex<double>);
  }
 };
--- a/poincare/include/poincare/frac_part.h
+++ b/poincare/include/poincare/frac_part.h
@@ -24,11 +24,11 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  template<typename T> static Complex<T> computeOnComplex(const Complex<T> c, AngleUnit angleUnit);
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  template<typename T> static std::complex<T> computeOnComplex(const std::complex<T> c, AngleUnit angleUnit);
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<float>(this, context, angleUnit, computeOnComplex<float>);
  }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<double>(this, context, angleUnit, computeOnComplex<double>);
  }
 };
--- a/poincare/include/poincare/global_context.h
+++ b/poincare/include/poincare/global_context.h
@@ -3,7 +3,7 @@

 #include <poincare/context.h>
 #include <poincare/matrix.h>
-#include <poincare/complex.h>
+#include <poincare/decimal.h>

 namespace Poincare {

@@ -28,15 +28,14 @@ public:
  static constexpr uint16_t k_maxNumberOfListExpressions = 10;
  static constexpr uint16_t k_maxNumberOfMatrixExpressions = 10;
 private:
-  static Complex<double> * defaultExpression();
+  static Decimal * defaultExpression();
  int symbolIndex(const Symbol * symbol) const;
-  Complex<double> * m_expressions[k_maxNumberOfScalarExpressions];
+  Expression * m_expressions[k_maxNumberOfScalarExpressions];
  Matrix * m_matrixExpressions[k_maxNumberOfMatrixExpressions];
  /* Matrix layout memoization */
  ExpressionLayout * m_matrixLayout[k_maxNumberOfMatrixExpressions];
-  Complex<double> m_pi;
-  Complex<double> m_e;
-  Complex<double> m_i;
+  Decimal m_pi;
+  Decimal m_e;
 };

 }
--- a/poincare/include/poincare/great_common_divisor.h
+++ b/poincare/include/poincare/great_common_divisor.h
@@ -3,7 +3,6 @@

 #include <poincare/layout_engine.h>
 #include <poincare/static_hierarchy.h>
-#include <poincare/complex.h>

 namespace Poincare {

@@ -24,8 +23,8 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
  template<typename T> Complex<T> * templatedApproximate(Context& context, AngleUnit angleUnit) const;
 };

--- a/poincare/include/poincare/hyperbolic_arc_cosine.h
+++ b/poincare/include/poincare/hyperbolic_arc_cosine.h
@@ -24,11 +24,11 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  template<typename T> static Complex<T> computeOnComplex(const Complex<T> c, AngleUnit angleUnit);
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  template<typename T> static std::complex<T> computeOnComplex(const std::complex<T> c, AngleUnit angleUnit);
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<float>(this, context, angleUnit,computeOnComplex<float>);
  }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<double>(this, context, angleUnit, computeOnComplex<double>);
  }
 };
--- a/poincare/include/poincare/hyperbolic_arc_sine.h
+++ b/poincare/include/poincare/hyperbolic_arc_sine.h
@@ -24,11 +24,11 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  template<typename T> static Complex<T> computeOnComplex(const Complex<T> c, AngleUnit angleUnit);
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  template<typename T> static std::complex<T> computeOnComplex(const std::complex<T> c, AngleUnit angleUnit);
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<float>(this, context, angleUnit,computeOnComplex<float>);
  }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<double>(this, context, angleUnit, computeOnComplex<double>);
  }
 };
--- a/poincare/include/poincare/hyperbolic_arc_tangent.h
+++ b/poincare/include/poincare/hyperbolic_arc_tangent.h
@@ -24,11 +24,11 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  template<typename T> static Complex<T> computeOnComplex(const Complex<T> c, AngleUnit angleUnit);
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  template<typename T> static std::complex<T> computeOnComplex(const std::complex<T> c, AngleUnit angleUnit);
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<float>(this, context, angleUnit,computeOnComplex<float>);
  }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<double>(this, context, angleUnit, computeOnComplex<double>);
  }
 };
--- a/poincare/include/poincare/hyperbolic_cosine.h
+++ b/poincare/include/poincare/hyperbolic_cosine.h
@@ -12,7 +12,7 @@ class HyperbolicCosine : public StaticHierarchy<1>  {
 public:
  Type type() const override;
  Expression * clone() const override;
-  template<typename T> static Complex<T> computeOnComplex(const Complex<T> c, AngleUnit angleUnit);
+  template<typename T> static std::complex<T> computeOnComplex(const std::complex<T> c, AngleUnit angleUnit);
 private:
  /* Layout */
  ExpressionLayout * privateCreateLayout(PrintFloat::Mode floatDisplayMode, ComplexFormat complexFormat) const override {
@@ -25,10 +25,10 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<float>(this, context, angleUnit,computeOnComplex<float>);
  }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<double>(this, context, angleUnit, computeOnComplex<double>);
  }
 };
--- a/poincare/include/poincare/hyperbolic_sine.h
+++ b/poincare/include/poincare/hyperbolic_sine.h
@@ -12,7 +12,7 @@ class HyperbolicSine : public StaticHierarchy<1>  {
 public:
  Type type() const override;
  Expression * clone() const override;
-  template<typename T> static Complex<T> computeOnComplex(const Complex<T> c, AngleUnit angleUnit);
+  template<typename T> static std::complex<T> computeOnComplex(const std::complex<T> c, AngleUnit angleUnit);
 private:
  /* Layout */
  ExpressionLayout * privateCreateLayout(PrintFloat::Mode floatDisplayMode, ComplexFormat complexFormat) const override {
@@ -25,10 +25,10 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<float>(this, context, angleUnit,computeOnComplex<float>);
  }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<double>(this, context, angleUnit, computeOnComplex<double>);
  }
 };
--- a/poincare/include/poincare/hyperbolic_tangent.h
+++ b/poincare/include/poincare/hyperbolic_tangent.h
@@ -12,7 +12,7 @@ class HyperbolicTangent : public StaticHierarchy<1>  {
 public:
  Type type() const override;
  Expression * clone() const override;
-  template<typename T> static Complex<T> computeOnComplex(const Complex<T> c, AngleUnit angleUnit);
+  template<typename T> static std::complex<T> computeOnComplex(const std::complex<T> c, AngleUnit angleUnit);
 private:
  /* Layout */
  ExpressionLayout * privateCreateLayout(PrintFloat::Mode floatDisplayMode, ComplexFormat complexFormat) const override {
@@ -25,10 +25,10 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<float>(this, context, angleUnit,computeOnComplex<float>);
  }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<double>(this, context, angleUnit, computeOnComplex<double>);
  }
 };
--- a/poincare/include/poincare/imaginary_part.h
+++ b/poincare/include/poincare/imaginary_part.h
@@ -24,11 +24,11 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  template<typename T> static Complex<T> computeOnComplex(const Complex<T> c, AngleUnit angleUnit);
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  template<typename T> static std::complex<T> computeOnComplex(const std::complex<T> c, AngleUnit angleUnit);
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<float>(this, context, angleUnit,computeOnComplex<float>);
  }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<double>(this, context, angleUnit, computeOnComplex<double>);
  }
 };
--- a/poincare/include/poincare/integral.h
+++ b/poincare/include/poincare/integral.h
@@ -4,7 +4,6 @@
 #include <poincare/static_hierarchy.h>
 #include <poincare/variable_context.h>
 #include <poincare/layout_engine.h>
-#include <poincare/complex.h>

 namespace Poincare {

@@ -23,8 +22,8 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
 template<typename T> Complex<T> * templatedApproximate(Context& context, AngleUnit angleUnit) const;
  template<typename T>
  struct DetailedResult
--- a/poincare/include/poincare/least_common_multiple.h
+++ b/poincare/include/poincare/least_common_multiple.h
@@ -3,7 +3,6 @@

 #include <poincare/layout_engine.h>
 #include <poincare/static_hierarchy.h>
-#include <poincare/complex.h>

 namespace Poincare {

@@ -24,8 +23,8 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
  template<typename T> Complex<T> * templatedApproximate(Context& context, AngleUnit angleUnit) const;
 };

--- a/poincare/include/poincare/logarithm.h
+++ b/poincare/include/poincare/logarithm.h
@@ -26,10 +26,10 @@ private:
  bool parentIsAPowerOfSameBase() const;
  Expression * splitInteger(Integer i, bool isDenominator, Context & context, AngleUnit angleUnit);
  /* Evaluation */
-  template<typename T> static Complex<T> computeOnComplex(const Complex<T> c, AngleUnit angleUnit);
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
-  template<typename T> Expression * templatedApproximate(Context& context, AngleUnit angleUnit) const;
+  template<typename T> static std::complex<T> computeOnComplex(const std::complex<T> c, AngleUnit angleUnit);
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
+  template<typename T> Evaluation<T> * templatedApproximate(Context& context, AngleUnit angleUnit) const;
 };

 }
--- a/poincare/include/poincare/matrix.h
+++ b/poincare/include/poincare/matrix.h
@@ -9,6 +9,7 @@ namespace Poincare {
 class Multiplication;

 class Matrix : public DynamicHierarchy {
+  template<typename T> friend class MatrixComplex;
 public:
  Matrix(MatrixData * matrixData); // pilfer the operands of matrixData
  Matrix(const Expression * const * operands, int numberOfRows, int numberOfColumns, bool cloneOperands = true);
@@ -24,32 +25,27 @@ public:

  /* Operation on matrix */
  int rank(Context & context, AngleUnit angleUnit, bool inPlace);
-  /* createInverse can be called on any matrix reduce or not, approximate or not. */
-  Expression * createInverse(Context & context, AngleUnit angleUnit) const;
-  /* createDeterminant and createTrace can only be called on a matrix of complex
-   * expressions. createDeterminant and createTrace return a complex
-   * expression. */
-  template<typename T> Complex<T> * createTrace() const;
-  template<typename T> Complex<T> * createDeterminant(Context & context, AngleUnit angleUnit) const;
-  /* createApproximateInverse has to be called on a matrix of complex and will
-   * return a matrix of complex if possible and nullptr otherwise. */
-  // TODO: createApproximateInverse should rather use ArrayInverse with T=complex<T>
-  template<typename T> Matrix * createApproximateInverse() const;
  // Inverse the array in-place. Array has to be given in the form array[row_index][column_index]
  template<typename T> static int ArrayInverse(T * array, int numberOfRows, int numberOfColumns);
+#if MATRIX_EXACT_REDUCING
+  //template<typename T> Expression * createTrace() const;
+  //template<typename T> Expression * createDeterminant() const;
  Matrix * createTranspose() const;
  static Matrix * createIdentity(int dim);
-  template<typename T> static Matrix * createApproximateIdentity(int dim);
+  /* createInverse can be called on any matrix reduce or not, approximate or not. */
+  Expression * createInverse(Context & context, AngleUnit angleUnit) const;
+#endif
 private:
  /* rowCanonize turns a matrix in its reduced row echelon form. */
  void rowCanonize(Context & context, AngleUnit angleUnit, Multiplication * m = nullptr);
-  template<typename T> static void ArrayRowCanonize(T * array, int numberOfRows, int numberOfColumns);
+  // Row canonize the array in place
+  template<typename T> static void ArrayRowCanonize(T * array, int numberOfRows, int numberOfColumns, T * c = nullptr);
  /* Layout */
  ExpressionLayout * privateCreateLayout(PrintFloat::Mode floatDisplayMode, ComplexFormat complexFormat) const override;
  /* Evaluation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
- template<typename T> Expression * templatedApproximate(Context& context, AngleUnit angleUnit) const;
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
+ template<typename T> Evaluation<T> * templatedApproximate(Context& context, AngleUnit angleUnit) const;
  int m_numberOfRows;
 };

--- a/poincare/include/poincare/matrix_dimension.h
+++ b/poincare/include/poincare/matrix_dimension.h
@@ -3,7 +3,6 @@

 #include <poincare/layout_engine.h>
 #include <poincare/static_hierarchy.h>
-#include <poincare/complex.h>

 namespace Poincare {

@@ -24,9 +23,9 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
- template<typename T> Expression * templatedApproximate(Context& context, AngleUnit angleUnit) const;
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
+ template<typename T> Evaluation<T> * templatedApproximate(Context& context, AngleUnit angleUnit) const;
 };

 }
--- a/poincare/include/poincare/matrix_inverse.h
+++ b/poincare/include/poincare/matrix_inverse.h
@@ -3,7 +3,6 @@

 #include <poincare/layout_engine.h>
 #include <poincare/static_hierarchy.h>
-#include <poincare/complex.h>

 namespace Poincare {

@@ -24,9 +23,9 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context & context, AngleUnit angleUnit) override;
  /* Evaluation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
-  template<typename T> Expression * templatedApproximate(Context& context, AngleUnit angleUnit) const;
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
+  template<typename T> Evaluation<T> * templatedApproximate(Context& context, AngleUnit angleUnit) const;
 };

 }
--- a/poincare/include/poincare/matrix_trace.h
+++ b/poincare/include/poincare/matrix_trace.h
@@ -3,7 +3,6 @@

 #include <poincare/layout_engine.h>
 #include <poincare/static_hierarchy.h>
-#include <poincare/complex.h>

 namespace Poincare {

@@ -24,9 +23,9 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context & context, AngleUnit angleUnit) override;
  /* Evaluation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
- template<typename T> Expression * templatedApproximate(Context& context, AngleUnit angleUnit) const;
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
+ template<typename T> Complex<T> * templatedApproximate(Context& context, AngleUnit angleUnit) const;
 };

 }
--- a/poincare/include/poincare/matrix_transpose.h
+++ b/poincare/include/poincare/matrix_transpose.h
@@ -3,7 +3,6 @@

 #include <poincare/layout_engine.h>
 #include <poincare/static_hierarchy.h>
-#include <poincare/complex.h>

 namespace Poincare {

@@ -24,9 +23,9 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context & context, AngleUnit angleUnit) override;
  /* Evaluation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
-  template<typename T> Expression * templatedApproximate(Context& context, AngleUnit angleUnit) const;
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
+  template<typename T> Evaluation<T> * templatedApproximate(Context& context, AngleUnit angleUnit) const;
 };

 }
--- a/poincare/include/poincare/multiplication.h
+++ b/poincare/include/poincare/multiplication.h
@@ -25,12 +25,12 @@ public:
  int polynomialDegree(char symbolName) const override;
  int privateGetPolynomialCoefficients(char symbolName, Expression * coefficients[]) const override;
  /* Evaluation */
-  template<typename T> static Complex<T> compute(const Complex<T> c, const Complex<T> d);
-  template<typename T> static Matrix * computeOnComplexAndMatrix(const Complex<T> * c, const Matrix * m) {
-    return ApproximationEngine::elementWiseOnComplexAndComplexMatrix(c, m, compute<T>);
+  template<typename T> static std::complex<T> compute(const std::complex<T> c, const std::complex<T> d);
+  template<typename T> static MatrixComplex<T> computeOnComplexAndMatrix(const std::complex<T> c, const MatrixComplex<T> m) {
+    return ApproximationEngine::elementWiseOnMatrixComplexAndComplex(m, c, compute<T>);
  }
-  template<typename T> static Matrix * computeOnMatrices(const Matrix * m, const Matrix * n);
  template<typename T> static void computeOnArrays(T * m, T * n, T * result, int mNumberOfColumns, int mNumberOfRows, int nNumberOfColumns);
+  template<typename T> static MatrixComplex<T> computeOnMatrices(const MatrixComplex<T> m, const MatrixComplex<T> n);
 private:
  /* Property */
  Expression * setSign(Sign s, Context & context, AngleUnit angleUnit) override;
@@ -59,13 +59,13 @@ private:
  Expression * mergeNegativePower(Context & context, AngleUnit angleUnit);
  /* Evaluation */

-  template<typename T> static Matrix * computeOnMatrixAndComplex(const Matrix * m, const Complex<T> * c) {
-    return ApproximationEngine::elementWiseOnComplexAndComplexMatrix(c, m, compute<T>);
+  template<typename T> static MatrixComplex<T> computeOnMatrixAndComplex(const MatrixComplex<T> m, const std::complex<T> c) {
+    return ApproximationEngine::elementWiseOnMatrixComplexAndComplex(m, c, compute<T>);
  }
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::mapReduce<float>(this, context, angleUnit, compute<float>, computeOnComplexAndMatrix<float>, computeOnMatrixAndComplex<float>, computeOnMatrices<float>);
  }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::mapReduce<double>(this, context, angleUnit, compute<double>, computeOnComplexAndMatrix<double>, computeOnMatrixAndComplex<double>, computeOnMatrices<double>);
  }
 };
--- a/poincare/include/poincare/naperian_logarithm.h
+++ b/poincare/include/poincare/naperian_logarithm.h
@@ -24,11 +24,11 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  template<typename T> static Complex<T> computeOnComplex(const Complex<T> c, AngleUnit angleUnit);
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  template<typename T> static std::complex<T> computeOnComplex(const std::complex<T> c, AngleUnit angleUnit);
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<float>(this, context, angleUnit,computeOnComplex<float>);
  }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<double>(this, context, angleUnit, computeOnComplex<double>);
  }
 };
--- a/poincare/include/poincare/nth_root.h
+++ b/poincare/include/poincare/nth_root.h
@@ -3,7 +3,6 @@

 #include <poincare/static_hierarchy.h>
 #include <poincare/layout_engine.h>
-#include <poincare/complex.h>

 namespace Poincare {

@@ -21,10 +20,9 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  template<typename T> static Complex<T> compute(const Complex<T> c, const Complex<T> d);
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
- template<typename T> Expression * templatedApproximate(Context& context, AngleUnit angleUnit) const;
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
+ template<typename T> Evaluation<T> * templatedApproximate(Context& context, AngleUnit angleUnit) const;

 };

--- a/poincare/include/poincare/opposite.h
+++ b/poincare/include/poincare/opposite.h
@@ -14,7 +14,7 @@ public:
  Type type() const override;
  int polynomialDegree(char symbolName) const override;
  Sign sign() const override;
-  template<typename T> static Complex<T> compute(const Complex<T> c, AngleUnit angleUnit);
+  template<typename T> static std::complex<T> compute(const std::complex<T> c, AngleUnit angleUnit);
 private:
  /* Layout */
  bool needParenthesisWithParent(const Expression * e) const override;
@@ -23,10 +23,10 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<float>(this, context, angleUnit, compute<float>);
  }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<double>(this, context, angleUnit, compute<double>);
  }
 };
--- a/poincare/include/poincare/parenthesis.h
+++ b/poincare/include/poincare/parenthesis.h
@@ -22,9 +22,9 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
- template<typename T> Expression * templatedApproximate(Context& context, AngleUnit angleUnit) const;
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
+ template<typename T> Evaluation<T> * templatedApproximate(Context& context, AngleUnit angleUnit) const;

 };

--- a/poincare/include/poincare/permute_coefficient.h
+++ b/poincare/include/poincare/permute_coefficient.h
@@ -3,7 +3,6 @@

 #include <poincare/layout_engine.h>
 #include <poincare/static_hierarchy.h>
-#include <poincare/complex.h>

 namespace Poincare {

@@ -25,8 +24,8 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
  template<typename T> Complex<T> * templatedApproximate(Context& context, AngleUnit angleUnit) const;
 };

--- a/poincare/include/poincare/power.h
+++ b/poincare/include/poincare/power.h
@@ -23,7 +23,7 @@ public:
  Sign sign() const override;
  int polynomialDegree(char symbolName) const override;
  int privateGetPolynomialCoefficients(char symbolName, Expression * coefficients[]) const override;
-  template<typename T> static Complex<T> compute(const Complex<T> c, const Complex<T> d);
+  template<typename T> static std::complex<T> compute(const std::complex<T> c, const std::complex<T> d);
 private:
  constexpr static int k_maxNumberOfTermsInExpandedMultinome = 25;
  constexpr static int k_maxExactPowerMatrix = 100;
@@ -56,13 +56,13 @@ private:
  static bool RationalExponentShouldNotBeReduced(const Rational * b, const Rational * r);
  /* Evaluation */
  constexpr static int k_maxApproximatePowerMatrix = 1000;
-  template<typename T> static Matrix * computeOnComplexAndMatrix(const Complex<T> * c, const Matrix * n) { return nullptr; }
-  template<typename T> static Matrix * computeOnMatrixAndComplex(const Matrix * m, const Complex<T> * d);
-  template<typename T> static Matrix * computeOnMatrices(const Matrix * m, const Matrix * n) { return nullptr; }
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  template<typename T> static MatrixComplex<T> computeOnComplexAndMatrix(const std::complex<T> c, const MatrixComplex<T> n);
+  template<typename T> static MatrixComplex<T> computeOnMatrixAndComplex(const MatrixComplex<T> m, const std::complex<T> d);
+  template<typename T> static MatrixComplex<T> computeOnMatrices(const MatrixComplex<T> m, const MatrixComplex<T> n);
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::mapReduce<float>(this, context, angleUnit, compute<float>, computeOnComplexAndMatrix<float>, computeOnMatrixAndComplex<float>, computeOnMatrices<float>);
  }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::mapReduce<double>(this, context, angleUnit, compute<double>, computeOnComplexAndMatrix<double>, computeOnMatrixAndComplex<double>, computeOnMatrices<double>);
  }
 };
--- a/poincare/include/poincare/prediction_interval.h
+++ b/poincare/include/poincare/prediction_interval.h
@@ -24,9 +24,9 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  Expression * privateApproximate(Expression::SinglePrecision p, Context& context, Expression::AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
-  Expression * privateApproximate(Expression::DoublePrecision p, Context& context, Expression::AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
-  template<typename T> Expression * templatedApproximate(Context& context, AngleUnit angleUnit) const;
+  Evaluation<float> * privateApproximate(Expression::SinglePrecision p, Context& context, Expression::AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
+  Evaluation<double> * privateApproximate(Expression::DoublePrecision p, Context& context, Expression::AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
+  template<typename T> Evaluation<T> * templatedApproximate(Context& context, AngleUnit angleUnit) const;
 };

 }
--- a/poincare/include/poincare/product.h
+++ b/poincare/include/poincare/product.h
@@ -14,13 +14,13 @@ private:
  const char * name() const override;
  int emptySequenceValue() const override;
  ExpressionLayout * createSequenceLayoutWithArgumentLayouts(ExpressionLayout * argumentLayout, ExpressionLayout * subscriptLayout, ExpressionLayout * superscriptLayout) const override;
-  Expression * evaluateWithNextTerm(DoublePrecision p, Expression * a, Expression * b) const override {
+  Evaluation<double> * evaluateWithNextTerm(DoublePrecision p, Evaluation<double> * a, Evaluation<double> * b) const override {
    return templatedApproximateWithNextTerm<double>(a, b);
  }
-  Expression * evaluateWithNextTerm(SinglePrecision p, Expression * a, Expression * b) const override {
+  Evaluation<float> * evaluateWithNextTerm(SinglePrecision p, Evaluation<float> * a, Evaluation<float> * b) const override {
    return templatedApproximateWithNextTerm<float>(a, b);
  }
-  template<typename T> Expression * templatedApproximateWithNextTerm(Expression * a, Expression * b) const;
+  template<typename T> Evaluation<T> * templatedApproximateWithNextTerm(Evaluation<T> * a, Evaluation<T> * b) const;
 };

 }
--- a/poincare/include/poincare/randint.h
+++ b/poincare/include/poincare/randint.h
@@ -21,13 +21,13 @@ private:
  }
  const char * name() const { return "randint"; }
  /* Evaluation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    return templateApproximate<float>(context, angleUnit);
  }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    return templateApproximate<double>(context, angleUnit);
  }
-  template <typename T> Expression * templateApproximate(Context& context, AngleUnit angleUnit) const;
+  template <typename T> Evaluation<T> * templateApproximate(Context& context, AngleUnit angleUnit) const;
 };

 }
--- a/poincare/include/poincare/random.h
+++ b/poincare/include/poincare/random.h
@@ -3,7 +3,6 @@

 #include <poincare/static_hierarchy.h>
 #include <poincare/layout_engine.h>
-#include <poincare/complex.h>

 namespace Poincare {

@@ -25,14 +24,14 @@ private:
  }
  const char * name() const { return "random"; }
  /* Evaluation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    return templateApproximate<float>();
  }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    return templateApproximate<double>();
  }
-  template <typename T> Expression * templateApproximate() const {
-    return new Complex<T>(Complex<T>::Float(random<T>()));
+  template <typename T> Evaluation<T> * templateApproximate() const {
+    return new Complex<T>(random<T>());
  }
 };

--- a/poincare/include/poincare/rational.h
+++ b/poincare/include/poincare/rational.h
@@ -3,7 +3,7 @@

 #include <poincare/integer.h>
 #include <poincare/static_hierarchy.h>
-#include <poincare/complex.h>
+#include <poincare/evaluation.h>

 namespace Poincare {

@@ -45,8 +45,8 @@ private:
  bool needParenthesisWithParent(const Expression * e) const override;
  ExpressionLayout * privateCreateLayout(PrintFloat::Mode floatDisplayMode, ComplexFormat complexFormat) const override;
  int writeTextInBuffer(char * buffer, int bufferSize, int numberOfSignificantDigits = PrintFloat::k_numberOfStoredSignificantDigits) const override;
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
  template<typename U> Complex<U> * templatedApproximate(Context& context, Expression::AngleUnit angleUnit) const;
  Expression * shallowBeautify(Context & context, AngleUnit angleUnit) override;
  Expression * setSign(Sign s);
--- a/poincare/include/poincare/real_part.h
+++ b/poincare/include/poincare/real_part.h
@@ -24,11 +24,11 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  template<typename T> static Complex<T> computeOnComplex(const Complex<T> c, AngleUnit angleUnit);
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  template<typename T> static std::complex<T> computeOnComplex(const std::complex<T> c, AngleUnit angleUnit);
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<float>(this, context, angleUnit,computeOnComplex<float>);
  }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<double>(this, context, angleUnit, computeOnComplex<double>);
  }
 };
--- a/poincare/include/poincare/round.h
+++ b/poincare/include/poincare/round.h
@@ -3,7 +3,6 @@

 #include <poincare/layout_engine.h>
 #include <poincare/static_hierarchy.h>
-#include <poincare/complex.h>

 namespace Poincare {

@@ -24,8 +23,8 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Complex */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
  template<typename T> Complex<T> * templatedApproximate(Context& context, AngleUnit angleUnit) const;
 };

--- a/poincare/include/poincare/sequence.h
+++ b/poincare/include/poincare/sequence.h
@@ -17,12 +17,12 @@ private:
  virtual ExpressionLayout * createSequenceLayoutWithArgumentLayouts(ExpressionLayout * subscriptLayout, ExpressionLayout * superscriptLayout, ExpressionLayout * argumentLayout) const = 0;
  virtual const char * name() const = 0;
  /* Evaluation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
- template<typename T> Expression * templatedApproximate(Context& context, AngleUnit angleUnit) const;
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
+ template<typename T> Evaluation<T> * templatedApproximate(Context& context, AngleUnit angleUnit) const;
  virtual int emptySequenceValue() const = 0;
-  virtual Expression * evaluateWithNextTerm(SinglePrecision p, Expression * a, Expression * b) const = 0;
-  virtual Expression * evaluateWithNextTerm(DoublePrecision p, Expression * a, Expression * b) const = 0;
+  virtual Evaluation<float> * evaluateWithNextTerm(SinglePrecision p, Evaluation<float> * a, Evaluation<float> * b) const = 0;
+  virtual Evaluation<double> * evaluateWithNextTerm(DoublePrecision p, Evaluation<double> * a, Evaluation<double> * b) const = 0;
 };

 }
--- a/poincare/include/poincare/simplification_root.h
+++ b/poincare/include/poincare/simplification_root.h
@@ -24,11 +24,11 @@ public:
    return nullptr;
  }
  int writeTextInBuffer(char * buffer, int bufferSize, int numberOfSignificantDigits = PrintFloat::k_numberOfStoredSignificantDigits) const override { return 0; }
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    assert(false);
    return nullptr;
  }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    assert(false);
    return nullptr;
  }
--- a/poincare/include/poincare/sine.h
+++ b/poincare/include/poincare/sine.h
@@ -15,7 +15,9 @@ class Sine : public StaticHierarchy<1> {
 public:
  Type type() const override;
  Expression * clone() const override;
-  template<typename T> static Complex<T> computeOnComplex(const Complex<T> c, AngleUnit angleUnit = AngleUnit::Radian);
+  template<typename T> static std::complex<T> computeOnComplex(const std::complex<T> c, AngleUnit angleUnit = AngleUnit::Radian) {
+    return Trigonometry::computeOnComplex(c, angleUnit, std::sin);
+  }
 private:
  /* Layout */
  ExpressionLayout * privateCreateLayout(PrintFloat::Mode floatDisplayMode, ComplexFormat complexFormat) const override {
@@ -28,10 +30,10 @@ private:
  /* Simplication */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<float>(this, context, angleUnit,computeOnComplex<float>);
  }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<double>(this, context, angleUnit, computeOnComplex<double>);
  }
 };
--- a/poincare/include/poincare/square_root.h
+++ b/poincare/include/poincare/square_root.h
@@ -19,11 +19,11 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  template<typename T> static Complex<T> computeOnComplex(const Complex<T> c, AngleUnit angleUnit);
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  template<typename T> static std::complex<T> computeOnComplex(const std::complex<T> c, AngleUnit angleUnit);
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<float>(this, context, angleUnit,computeOnComplex<float>);
  }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<double>(this, context, angleUnit, computeOnComplex<double>);
  }

--- a/poincare/include/poincare/store.h
+++ b/poincare/include/poincare/store.h
@@ -21,9 +21,9 @@ private:
  ExpressionLayout * privateCreateLayout(PrintFloat::Mode floatDisplayMode, ComplexFormat complexFormat) const override;
  int writeTextInBuffer(char * buffer, int bufferSize, int numberOfSignificantDigits = PrintFloat::k_numberOfStoredSignificantDigits) const override;
  /* Evalutation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
-  template<typename T> Expression * templatedApproximate(Context& context, AngleUnit angleUnit) const;
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
+  template<typename T> Evaluation<T> * templatedApproximate(Context& context, AngleUnit angleUnit) const;

  const Symbol * symbol() const { return static_cast<const Symbol *>(operand(1)); }
  const Expression * value() const { return operand(0); }
--- a/poincare/include/poincare/subtraction.h
+++ b/poincare/include/poincare/subtraction.h
@@ -12,7 +12,7 @@ class Subtraction : public StaticHierarchy<2> {
 public:
  Type type() const override;
  Expression * clone() const override;
-  template<typename T> static Complex<T> compute(const Complex<T> c, const Complex<T> d);
+  template<typename T> static std::complex<T> compute(const std::complex<T> c, const std::complex<T> d);
  int polynomialDegree(char symbolName) const override;
 private:
  /* Layout */
@@ -27,18 +27,18 @@ private:
  /* Simplification */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  template<typename T> static Matrix * computeOnMatrixAndComplex(const Matrix * m, const Complex<T> * c) {
-    return ApproximationEngine::elementWiseOnComplexAndComplexMatrix(c, m, compute<T>);
+  template<typename T> static MatrixComplex<T> computeOnMatrixAndComplex(const MatrixComplex<T> m, const std::complex<T> c) {
+    return ApproximationEngine::elementWiseOnMatrixComplexAndComplex(m, c, compute<T>);
  }
-  template<typename T> static Matrix * computeOnComplexAndMatrix(const Complex<T> * c, const Matrix * n);
-  template<typename T> static Matrix * computeOnMatrices(const Matrix * m, const Matrix * n) {
+  template<typename T> static MatrixComplex<T> computeOnComplexAndMatrix(const std::complex<T> c, const MatrixComplex<T> n);
+  template<typename T> static MatrixComplex<T> computeOnMatrices(const MatrixComplex<T> m, const MatrixComplex<T> n) {
    return ApproximationEngine::elementWiseOnComplexMatrices(m, n, compute<T>);
  }

-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::mapReduce<float>(this, context, angleUnit, compute<float>, computeOnComplexAndMatrix<float>, computeOnMatrixAndComplex<float>, computeOnMatrices<float>);
  }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::mapReduce<double>(this, context, angleUnit, compute<double>, computeOnComplexAndMatrix<double>, computeOnMatrixAndComplex<double>, computeOnMatrices<double>);
  }
 };
--- a/poincare/include/poincare/sum.h
+++ b/poincare/include/poincare/sum.h
@@ -14,13 +14,13 @@ private:
  const char * name() const override;
  int emptySequenceValue() const override;
  ExpressionLayout * createSequenceLayoutWithArgumentLayouts(ExpressionLayout * argumentLayout, ExpressionLayout * subscriptLayout, ExpressionLayout * superscriptLayout) const override;
-  Expression * evaluateWithNextTerm(DoublePrecision p, Expression * a, Expression * b) const override {
+  Evaluation<double> * evaluateWithNextTerm(DoublePrecision p, Evaluation<double> * a, Evaluation<double> * b) const override {
    return templatedApproximateWithNextTerm<double>(a, b);
  }
-  Expression * evaluateWithNextTerm(SinglePrecision p, Expression * a, Expression * b) const override {
+  Evaluation<float> * evaluateWithNextTerm(SinglePrecision p, Evaluation<float> * a, Evaluation<float> * b) const override {
    return templatedApproximateWithNextTerm<float>(a, b);
  }
-  template<typename T> Expression * templatedApproximateWithNextTerm(Expression * a, Expression * b) const;
+  template<typename T> Evaluation<T> * templatedApproximateWithNextTerm(Evaluation<T> * a, Evaluation<T> * b) const;
 };

 }
--- a/poincare/include/poincare/symbol.h
+++ b/poincare/include/poincare/symbol.h
@@ -70,9 +70,9 @@ private:
  ExpressionLayout * privateCreateLayout(PrintFloat::Mode floatDisplayMode, ComplexFormat complexFormat) const override;
  int writeTextInBuffer(char * buffer, int bufferSize, int numberOfSignificantDigits = PrintFloat::k_numberOfStoredSignificantDigits) const override;
  /* Evaluation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
- template<typename T> Expression * templatedApproximate(Context& context, AngleUnit angleUnit) const;
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
+ template<typename T> Evaluation<T> * templatedApproximate(Context& context, AngleUnit angleUnit) const;
  const char m_name;
 };

--- a/poincare/include/poincare/tangent.h
+++ b/poincare/include/poincare/tangent.h
@@ -4,7 +4,7 @@
 #include <poincare/layout_engine.h>
 #include <poincare/static_hierarchy.h>
 #include <poincare/approximation_engine.h>
-
+#include <poincare/trigonometry.h>

 namespace Poincare {

@@ -26,11 +26,13 @@ private:
  /* Simplication */
  Expression * shallowReduce(Context& context, AngleUnit angleUnit) override;
  /* Evaluation */
-  template<typename T> static Complex<T> computeOnComplex(const Complex<T> c, AngleUnit angleUnit = AngleUnit::Radian);
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
+  template<typename T> static std::complex<T> computeOnComplex(const std::complex<T> c, AngleUnit angleUnit = AngleUnit::Radian) {
+    return Trigonometry::computeOnComplex(c, angleUnit, std::tan);
+  }
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<float>(this, context, angleUnit,computeOnComplex<float>);
  }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override {
    return ApproximationEngine::map<double>(this, context, angleUnit, computeOnComplex<double>);
  }
 };
--- a/poincare/include/poincare/trigonometry.h
+++ b/poincare/include/poincare/trigonometry.h
@@ -17,6 +17,8 @@ public:
  static bool ExpressionIsEquivalentToTangent(const Expression * e);
  constexpr static int k_numberOfEntries = 37;
  static Expression * table(const Expression * e, Expression::Type type, Context & context, Expression::AngleUnit angleUnit); // , Function f, bool inverse
+  template <typename T> using Approximation = std::complex<T> (*)(const std::complex<T>&);
+  template <typename T> static std::complex<T> computeOnComplex(const std::complex<T> c, Expression::AngleUnit angleUnit, Approximation<T> approximate);
 };

 }
--- a/poincare/include/poincare/undefined.h
+++ b/poincare/include/poincare/undefined.h
@@ -2,7 +2,7 @@
 #define POINCARE_UNDEFINED_H

 #include <poincare/static_hierarchy.h>
-#include <poincare/complex.h>
+#include <poincare/evaluation.h>

 namespace Poincare {

@@ -16,8 +16,8 @@ private:
  /* Layout */
  ExpressionLayout * privateCreateLayout(PrintFloat::Mode floatDisplayMode, ComplexFormat complexFormat) const override;
  /* Evaluation */
-  Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
-  Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
+  Evaluation<float> * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
+  Evaluation<double> * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
  template<typename T> Complex<T> * templatedApproximate(Context& context, AngleUnit angleUnit) const;
 };

--- a/poincare/include/poincare/variable_context.h
+++ b/poincare/include/poincare/variable_context.h
@@ -2,19 +2,29 @@
 #define POINCARE_VARIABLE_CONTEXT_H

 #include <poincare/context.h>
-#include <poincare/complex.h>
+#include <poincare/evaluation.h>

 namespace Poincare {

-template<typename T>
 class VariableContext : public Context {
 public:
  VariableContext(char name, Context * parentContext = nullptr);
+  ~VariableContext();
+  /* This function could be avoid by doing:
+   * {
+   *   Expression * e = Expression::CreateDecimal<T>(value);
+   *   variableContext.setExpressionForSymbolName(e, &xSymbol, variableContext);
+   *   delete e;
+   * }
+   * But this would imply to allocated twice the expression (as e is cloned in
+   * setExpressionForSymbolName).
+   */
+  template<typename T> void setApproximationForVariable(T value);
  void setExpressionForSymbolName(const Expression * expression, const Symbol * symbol, Context & context) override;
  const Expression * expressionForSymbol(const Symbol * symbol) override;
 private:
  char m_name;
-  Complex<T> m_value;
+  Expression * m_value;
  Context * m_parentContext;
 };

--- a/poincare/src/absolute_value.cpp
+++ b/poincare/src/absolute_value.cpp
@@ -1,5 +1,4 @@
 #include <poincare/absolute_value.h>
-#include <poincare/complex.h>
 #include <poincare/simplification_engine.h>
 #include "layout/absolute_value_layout.h"

@@ -52,8 +51,8 @@ Expression * AbsoluteValue::shallowReduce(Context& context, AngleUnit angleUnit)
 }

 template<typename T>
-Complex<T> AbsoluteValue::computeOnComplex(const Complex<T> c, AngleUnit angleUnit) {
-  return Complex<T>::Float(c.r());
+std::complex<T> AbsoluteValue::computeOnComplex(const std::complex<T> c, AngleUnit angleUnit) {
+  return std::abs(c);
 }

 }
--- a/poincare/src/addition.cpp
+++ b/poincare/src/addition.cpp
@@ -347,17 +347,17 @@ Expression * Addition::shallowBeautify(Context & context, AngleUnit angleUnit) {
 /* Evaluation */

 template<typename T>
-Complex<T> Addition::compute(const Complex<T> c, const Complex<T> d) {
-  return Complex<T>::Cartesian(c.a()+d.a(), c.b()+d.b());
+std::complex<T> Addition::compute(const std::complex<T> c, const std::complex<T> d) {
+  return c+d;
 }

-template Complex<float> Poincare::Addition::compute<float>(Poincare::Complex<float>, Poincare::Complex<float>);
-template Complex<double> Poincare::Addition::compute<double>(Poincare::Complex<double>, Poincare::Complex<double>);
+template std::complex<float> Poincare::Addition::compute<float>(std::complex<float>, std::complex<float>);
+template std::complex<double> Poincare::Addition::compute<double>(std::complex<double>, std::complex<double>);

-template Matrix* Addition::computeOnMatrices<float>(const Matrix*,const Matrix*);
-template Matrix* Addition::computeOnMatrices<double>(const Matrix*,const Matrix*);
+template MatrixComplex<float> Addition::computeOnMatrices<float>(const MatrixComplex<float>,const MatrixComplex<float>);
+template MatrixComplex<double> Addition::computeOnMatrices<double>(const MatrixComplex<double>,const MatrixComplex<double>);

-template Matrix* Addition::computeOnComplexAndMatrix<float>(Complex<float> const*, const Matrix*);
-template Matrix* Addition::computeOnComplexAndMatrix<double>(Complex<double> const*, const Matrix*);
+template MatrixComplex<float> Addition::computeOnComplexAndMatrix<float>(std::complex<float> const, const MatrixComplex<float>);
+template MatrixComplex<double> Addition::computeOnComplexAndMatrix<double>(std::complex<double> const, const MatrixComplex<double>);

 }
--- a/poincare/src/approximation_engine.cpp
+++ b/poincare/src/approximation_engine.cpp
@@ -1,4 +1,5 @@
 #include <poincare/approximation_engine.h>
+#include <poincare/evaluation.h>
 #include <poincare/matrix.h>
 #include <cmath>
 extern "C" {
@@ -7,98 +8,99 @@ extern "C" {

 namespace Poincare {

-template<typename T> Expression * ApproximationEngine::map(const Expression * expression, Context& context, Expression::AngleUnit angleUnit, ComplexCompute<T> compute) {
+template<typename T> Evaluation<T> * ApproximationEngine::map(const Expression * expression, Context& context, Expression::AngleUnit angleUnit, ComplexCompute<T> compute) {
  assert(expression->numberOfOperands() == 1);
-  Expression * input = expression->operand(0)->approximate<T>(context, angleUnit);
-  Expression * result = nullptr;
-  if (input->type() == Expression::Type::Complex) {
+  Evaluation<T> * input = expression->operand(0)->privateApproximate(T(), context, angleUnit);
+  Evaluation<T> * result = nullptr;
+  if (input->type() == Evaluation<T>::Type::Complex) {
    Complex<T> * c = static_cast<Complex<T> *>(input);
    result = new Complex<T>(compute(*c, angleUnit));
  } else {
-    assert(input->type() == Expression::Type::Matrix);
-    Expression ** operands = new Expression * [input->numberOfOperands()];
-    for (int i = 0; i < input->numberOfOperands(); i++) {
-      assert(input->operand(i)->type() == Expression::Type::Complex);
-      const Complex<T> * c = static_cast<const Complex<T> *>(input->operand(i));
-      operands[i] = new Complex<T>(compute(*c, angleUnit));
+    assert(input->type() == Evaluation<T>::Type::MatrixComplex);
+    MatrixComplex<T> * m = static_cast<MatrixComplex<T> *>(input);
+    std::complex<T> * operands = new std::complex<T> [m->numberOfComplexOperands()];
+    for (int i = 0; i < m->numberOfComplexOperands(); i++) {
+      const std::complex<T> c = m->complexOperand(i);
+      operands[i] = compute(c, angleUnit);
    }
-    result = new Matrix(operands, static_cast<Matrix *>(input)->numberOfRows(), static_cast<Matrix *>(input)->numberOfColumns(), false);
+    result = new MatrixComplex<T>(operands, m->numberOfRows(), m->numberOfColumns());
    delete[] operands;
  }
  delete input;
  return result;
 }

-template<typename T> Expression * ApproximationEngine::mapReduce(const Expression * expression, Context& context, Expression::AngleUnit angleUnit, ComplexAndComplexReduction<T> computeOnComplexes, ComplexAndMatrixReduction<T> computeOnComplexAndMatrix, MatrixAndComplexReduction<T> computeOnMatrixAndComplex, MatrixAndMatrixReduction<T> computeOnMatrices) {
-  Expression * result = expression->operand(0)->approximate<T>(context, angleUnit);
+template<typename T> Evaluation<T> * ApproximationEngine::mapReduce(const Expression * expression, Context& context, Expression::AngleUnit angleUnit, ComplexAndComplexReduction<T> computeOnComplexes, ComplexAndMatrixReduction<T> computeOnComplexAndMatrix, MatrixAndComplexReduction<T> computeOnMatrixAndComplex, MatrixAndMatrixReduction<T> computeOnMatrices) {
+  Evaluation<T> * result = expression->operand(0)->privateApproximate(T(), context, angleUnit);
  for (int i = 1; i < expression->numberOfOperands(); i++) {
-    Expression * intermediateResult = nullptr;
-    Expression * nextOperandEvaluation = expression->operand(i)->approximate<T>(context, angleUnit);
-    if (result->type() == Expression::Type::Complex && nextOperandEvaluation->type() == Expression::Type::Complex) {
+    Evaluation<T> * intermediateResult = nullptr;
+    Evaluation<T> * nextOperandEvaluation = expression->operand(i)->privateApproximate(T(), context, angleUnit);
+    if (result->type() == Evaluation<T>::Type::Complex && nextOperandEvaluation->type() == Evaluation<T>::Type::Complex) {
      const Complex<T> * c = static_cast<const Complex<T> *>(result);
      const Complex<T> * d = static_cast<const Complex<T> *>(nextOperandEvaluation);
      intermediateResult = new Complex<T>(computeOnComplexes(*c, *d));
-    } else if (result->type() == Expression::Type::Complex) {
+    } else if (result->type() == Evaluation<T>::Type::Complex) {
      const Complex<T> * c = static_cast<const Complex<T> *>(result);
-      assert(nextOperandEvaluation->type() == Expression::Type::Matrix);
-      const Matrix * n = static_cast<const Matrix *>(nextOperandEvaluation);
-      intermediateResult = computeOnComplexAndMatrix(c, n);
-    } else if (nextOperandEvaluation->type() == Expression::Type::Complex) {
-      assert(result->type() == Expression::Type::Matrix);
-      const Matrix * m = static_cast<const Matrix *>(result);
+      assert(nextOperandEvaluation->type() == Evaluation<T>::Type::MatrixComplex);
+      const MatrixComplex<T> * n = static_cast<const MatrixComplex<T> *>(nextOperandEvaluation);
+      intermediateResult = new MatrixComplex<T>(computeOnComplexAndMatrix(*c, *n));
+    } else if (nextOperandEvaluation->type() == Evaluation<T>::Type::Complex) {
+      assert(result->type() == Evaluation<T>::Type::MatrixComplex);
+      const MatrixComplex<T> * m = static_cast<const MatrixComplex<T> *>(result);
      const Complex<T> * d = static_cast<const Complex<T> *>(nextOperandEvaluation);
-      intermediateResult = computeOnMatrixAndComplex(m, d);
+      intermediateResult = new MatrixComplex<T>(computeOnMatrixAndComplex(*m, *d));
    } else {
-      assert(result->type() == Expression::Type::Matrix);
-      const Matrix * m = static_cast<const Matrix *>(result);
-      assert(nextOperandEvaluation->type() == Expression::Type::Matrix);
-      const Matrix * n = static_cast<const Matrix *>(nextOperandEvaluation);
-      intermediateResult = computeOnMatrices(m, n);
+      assert(result->type() == Evaluation<T>::Type::MatrixComplex);
+      const MatrixComplex<T> * m = static_cast<const MatrixComplex<T> *>(result);
+      assert(nextOperandEvaluation->type() == Evaluation<T>::Type::MatrixComplex);
+      const MatrixComplex<T> * n = static_cast<const MatrixComplex<T> *>(nextOperandEvaluation);
+      intermediateResult = new MatrixComplex<T>(computeOnMatrices(*m, *n));
    }
    delete result;
    delete nextOperandEvaluation;
    result = intermediateResult;
-    if (result == nullptr) {
-      return new Complex<T>(Complex<T>::Float(NAN));
+    assert(result != nullptr);
+    if (result->isUndefined()) {
+      return new Complex<T>(Complex<T>::Undefined());
    }
  }
  return result;
 }

-template<typename T> Matrix * ApproximationEngine::elementWiseOnComplexAndComplexMatrix(const Complex<T> * c, const Matrix * m, ComplexAndComplexReduction<T> computeOnComplexes) {
-  Expression ** operands = new Expression * [m->numberOfRows()*m->numberOfColumns()];
-  for (int i = 0; i < m->numberOfOperands(); i++) {
-    const Complex<T> * d = static_cast<const Complex<T> *>(m->operand(i));
-    operands[i] = new Complex<T>(computeOnComplexes(*d, *c));
+template<typename T> MatrixComplex<T> ApproximationEngine::elementWiseOnMatrixComplexAndComplex(const MatrixComplex<T> m, const std::complex<T> c, ComplexAndComplexReduction<T> computeOnComplexes) {
+  std::complex<T> * operands = new std::complex<T> [m.numberOfRows()*m.numberOfColumns()];
+  for (int i = 0; i < m.numberOfComplexOperands(); i++) {
+    const std::complex<T> d = m.complexOperand(i);
+    operands[i] = computeOnComplexes(d, c);
  }
-  Matrix * result = new Matrix(operands, m->numberOfRows(), m->numberOfColumns(), false);
+  MatrixComplex<T> result = MatrixComplex<T>(operands, m.numberOfRows(), m.numberOfColumns());
  delete[] operands;
  return result;
 }

-template<typename T> Matrix * ApproximationEngine::elementWiseOnComplexMatrices(const Matrix * m, const Matrix * n, ComplexAndComplexReduction<T> computeOnComplexes) {
-  if (m->numberOfRows() != n->numberOfRows() || m->numberOfColumns() != n->numberOfColumns()) {
-    return nullptr;
+template<typename T> MatrixComplex<T> ApproximationEngine::elementWiseOnComplexMatrices(const MatrixComplex<T> m, const MatrixComplex<T> n, ComplexAndComplexReduction<T> computeOnComplexes) {
+  if (m.numberOfRows() != n.numberOfRows() || m.numberOfColumns() != n.numberOfColumns()) {
+    return MatrixComplex<T>::Undefined();
  }
-  Expression ** operands = new Expression * [m->numberOfRows()*m->numberOfColumns()];
-  for (int i = 0; i < m->numberOfOperands(); i++) {
-    const Complex<T> * c = static_cast<const Complex<T> *>(m->operand(i));
-    const Complex<T> * d = static_cast<const Complex<T> *>(n->operand(i));
-    operands[i] = new Complex<T>(computeOnComplexes(*c, *d));
+  std::complex<T> * operands = new std::complex<T> [m.numberOfRows()*m.numberOfColumns()];
+  for (int i = 0; i < m.numberOfComplexOperands(); i++) {
+    const Complex<T> c = m.complexOperand(i);
+    const Complex<T> d = n.complexOperand(i);
+    operands[i] = computeOnComplexes(c, d);
  }
-  Matrix * result = new Matrix(operands, m->numberOfRows(), m->numberOfColumns(), false);
+  MatrixComplex<T> result = MatrixComplex<T>(operands, m.numberOfRows(), m.numberOfColumns());
  delete[] operands;
  return result;
 }

-template Poincare::Expression * Poincare::ApproximationEngine::map(const Poincare::Expression * expression, Poincare::Context& context, Poincare::Expression::AngleUnit angleUnit, Poincare::ApproximationEngine::ComplexCompute<float> compute);
-template Poincare::Expression * Poincare::ApproximationEngine::map(const Poincare::Expression * expression, Poincare::Context& context, Poincare::Expression::AngleUnit angleUnit, Poincare::ApproximationEngine::ComplexCompute<double> compute);
-template Poincare::Expression * Poincare::ApproximationEngine::mapReduce(const Poincare::Expression * expression, Poincare::Context& context, Poincare::Expression::AngleUnit angleUnit, Poincare::ApproximationEngine::ComplexAndComplexReduction<float> computeOnComplexes, Poincare::ApproximationEngine::ComplexAndMatrixReduction<float> computeOnComplexAndMatrix, Poincare::ApproximationEngine::MatrixAndComplexReduction<float> computeOnMatrixAndComplex, Poincare::ApproximationEngine::MatrixAndMatrixReduction<float> computeOnMatrices);
-template Poincare::Expression * Poincare::ApproximationEngine::mapReduce(const Poincare::Expression * expression, Poincare::Context& context, Poincare::Expression::AngleUnit angleUnit, Poincare::ApproximationEngine::ComplexAndComplexReduction<double> computeOnComplexes, Poincare::ApproximationEngine::ComplexAndMatrixReduction<double> computeOnComplexAndMatrix, Poincare::ApproximationEngine::MatrixAndComplexReduction<double> computeOnMatrixAndComplex, Poincare::ApproximationEngine::MatrixAndMatrixReduction<double> computeOnMatrices);
-template Poincare::Matrix * Poincare::ApproximationEngine::elementWiseOnComplexAndComplexMatrix<float>(Poincare::Complex<float> const*, const Poincare::Matrix *, Poincare::Complex<float> (*)(Poincare::Complex<float>, Poincare::Complex<float>));
-template Poincare::Matrix* Poincare::ApproximationEngine::elementWiseOnComplexAndComplexMatrix<double>(Poincare::Complex<double> const*, const Poincare::Matrix*, Poincare::Complex<double> (*)(Poincare::Complex<double>, Poincare::Complex<double>));
-template Poincare::Matrix* Poincare::ApproximationEngine::elementWiseOnComplexMatrices<float>(const Poincare::Matrix*, const Poincare::Matrix*, Poincare::Complex<float> (*)(Poincare::Complex<float>, Poincare::Complex<float>));
-template Poincare::Matrix* Poincare::ApproximationEngine::elementWiseOnComplexMatrices<double>(const Poincare::Matrix*, const Poincare::Matrix*, Poincare::Complex<double> (*)(Poincare::Complex<double>, Poincare::Complex<double>));
+template Poincare::Evaluation<float> * Poincare::ApproximationEngine::map(const Poincare::Expression * expression, Poincare::Context& context, Poincare::Expression::AngleUnit angleUnit, Poincare::ApproximationEngine::ComplexCompute<float> compute);
+template Poincare::Evaluation<double> * Poincare::ApproximationEngine::map(const Poincare::Expression * expression, Poincare::Context& context, Poincare::Expression::AngleUnit angleUnit, Poincare::ApproximationEngine::ComplexCompute<double> compute);
+template Poincare::Evaluation<float> * Poincare::ApproximationEngine::mapReduce(const Poincare::Expression * expression, Poincare::Context& context, Poincare::Expression::AngleUnit angleUnit, Poincare::ApproximationEngine::ComplexAndComplexReduction<float> computeOnComplexes, Poincare::ApproximationEngine::ComplexAndMatrixReduction<float> computeOnComplexAndMatrix, Poincare::ApproximationEngine::MatrixAndComplexReduction<float> computeOnMatrixAndComplex, Poincare::ApproximationEngine::MatrixAndMatrixReduction<float> computeOnMatrices);
+template Poincare::Evaluation<double> * Poincare::ApproximationEngine::mapReduce(const Poincare::Expression * expression, Poincare::Context& context, Poincare::Expression::AngleUnit angleUnit, Poincare::ApproximationEngine::ComplexAndComplexReduction<double> computeOnComplexes, Poincare::ApproximationEngine::ComplexAndMatrixReduction<double> computeOnComplexAndMatrix, Poincare::ApproximationEngine::MatrixAndComplexReduction<double> computeOnMatrixAndComplex, Poincare::ApproximationEngine::MatrixAndMatrixReduction<double> computeOnMatrices);
+template Poincare::MatrixComplex<float> Poincare::ApproximationEngine::elementWiseOnMatrixComplexAndComplex<float>(const Poincare::MatrixComplex<float>, const std::complex<float>, std::complex<float> (*)(std::complex<float>, std::complex<float>));
+template Poincare::MatrixComplex<double> Poincare::ApproximationEngine::elementWiseOnMatrixComplexAndComplex<double>(const Poincare::MatrixComplex<double>, std::complex<double> const, std::complex<double> (*)(std::complex<double>, std::complex<double>));
+template Poincare::MatrixComplex<float> Poincare::ApproximationEngine::elementWiseOnComplexMatrices<float>(const Poincare::MatrixComplex<float>, const Poincare::MatrixComplex<float>, std::complex<float> (*)(std::complex<float>, std::complex<float>));
+template Poincare::MatrixComplex<double> Poincare::ApproximationEngine::elementWiseOnComplexMatrices<double>(const Poincare::MatrixComplex<double>, const Poincare::MatrixComplex<double>, std::complex<double> (*)(std::complex<double>, std::complex<double>));


 }
--- a/poincare/src/arc_cosine.cpp
+++ b/poincare/src/arc_cosine.cpp
@@ -31,16 +31,13 @@ Expression * ArcCosine::shallowReduce(Context& context, AngleUnit angleUnit) {
 }

 template<typename T>
-Complex<T> ArcCosine::computeOnComplex(const Complex<T> c, AngleUnit angleUnit) {
+std::complex<T> ArcCosine::computeOnComplex(const std::complex<T> c, AngleUnit angleUnit) {
  assert(angleUnit != AngleUnit::Default);
-  if (c.b() != 0) {
-    return Complex<T>::Float(NAN);
+  std::complex<T> result = std::acos(c);
+  if (angleUnit == AngleUnit::Degree && result.imag() == 0.0) {
+    result *= 180/M_PI;
  }
-  T result = std::acos(c.a());
-  if (angleUnit == AngleUnit::Degree) {
-    return Complex<T>::Float(result*180/M_PI);
-  }
-  return Complex<T>::Float(result);
+  return result;
 }

 }
--- a/poincare/src/arc_sine.cpp
+++ b/poincare/src/arc_sine.cpp
@@ -31,16 +31,13 @@ Expression * ArcSine::shallowReduce(Context& context, AngleUnit angleUnit) {
 }

 template<typename T>
-Complex<T> ArcSine::computeOnComplex(const Complex<T> c, AngleUnit angleUnit) {
+std::complex<T> ArcSine::computeOnComplex(const std::complex<T> c, AngleUnit angleUnit) {
  assert(angleUnit != AngleUnit::Default);
-  if (c.b() != 0) {
-    return Complex<T>::Float(NAN);
+  std::complex<T> result = std::asin(c);
+  if (angleUnit == AngleUnit::Degree && result.imag() == 0.0) {
+    result *= 180/M_PI;
  }
-  T result = std::asin(c.a());
-  if (angleUnit == AngleUnit::Degree) {
-    return Complex<T>::Float(result*180/M_PI);
-  }
-  return Complex<T>::Float(result);
+  return result;
 }

 }
--- a/poincare/src/arc_tangent.cpp
+++ b/poincare/src/arc_tangent.cpp
@@ -31,16 +31,13 @@ Expression * ArcTangent::shallowReduce(Context& context, AngleUnit angleUnit) {
 }

 template<typename T>
-Complex<T> ArcTangent::computeOnComplex(const Complex<T> c, AngleUnit angleUnit) {
+std::complex<T> ArcTangent::computeOnComplex(const std::complex<T> c, AngleUnit angleUnit) {
  assert(angleUnit != AngleUnit::Default);
-  if (c.b() != 0) {
-    return Complex<T>::Float(NAN);
+  std::complex<T> result = std::atan(c);
+  if (angleUnit == AngleUnit::Degree && result.imag() == 0.0) {
+    result *= 180/M_PI;
  }
-  T result = std::atan(c.a());
-  if (angleUnit == AngleUnit::Degree) {
-    return Complex<T>::Float(result*180/M_PI);
-  }
-  return Complex<T>::Float(result);
+  return result;
 }

 }
--- a/poincare/src/binomial_coefficient.cpp
+++ b/poincare/src/binomial_coefficient.cpp
@@ -1,6 +1,5 @@
 #include <poincare/binomial_coefficient.h>
 #include <poincare/undefined.h>
-#include <poincare/complex.h>
 #include <poincare/rational.h>
 #include "layout/binomial_coefficient_layout.h"

@@ -82,17 +81,14 @@ ExpressionLayout * BinomialCoefficient::privateCreateLayout(PrintFloat::Mode flo
 }

 template<typename T>
-Expression * BinomialCoefficient::templatedApproximate(Context& context, AngleUnit angleUnit) const {
-  Expression * nInput = operand(0)->approximate<T>(context, angleUnit);
-  Expression * kInput = operand(1)->approximate<T>(context, angleUnit);
-  if (nInput->type() != Type::Complex || kInput->type() != Type::Complex) {
-    return new Complex<T>(Complex<T>::Float(NAN));
-  }
-  T n = static_cast<Complex<T> *>(nInput)->toScalar();
-  T k = static_cast<Complex<T> *>(kInput)->toScalar();
+Complex<T> * BinomialCoefficient::templatedApproximate(Context& context, AngleUnit angleUnit) const {
+  Evaluation<T> * nInput = operand(0)->privateApproximate(T(), context, angleUnit);
+  Evaluation<T> * kInput = operand(1)->privateApproximate(T(), context, angleUnit);
+  T n = nInput->toScalar();
+  T k = kInput->toScalar();
  delete nInput;
  delete kInput;
-  return new Complex<T>(Complex<T>::Float(compute(k, n)));
+  return new Complex<T>(compute(k, n));
 }


--- a/poincare/src/ceiling.cpp
+++ b/poincare/src/ceiling.cpp
@@ -52,11 +52,11 @@ Expression * Ceiling::shallowReduce(Context& context, AngleUnit angleUnit) {
 }

 template<typename T>
-Complex<T> Ceiling::computeOnComplex(const Complex<T> c, AngleUnit angleUnit) {
-  if (c.b() != 0) {
-    return Complex<T>::Float(NAN);
+std::complex<T> Ceiling::computeOnComplex(const std::complex<T> c, AngleUnit angleUnit) {
+  if (c.imag() != 0) {
+    return Complex<T>::Undefined();
  }
-  return Complex<T>::Float(std::ceil(c.a()));
+  return std::ceil(c.real());
 }

 ExpressionLayout * Ceiling::privateCreateLayout(PrintFloat::Mode floatDisplayMode, ComplexFormat complexFormat) const {
--- a/poincare/src/complex.cpp
+++ b/poincare/src/complex.cpp
@@ -1,362 +0,0 @@
-#include <poincare/complex.h>
-#include <poincare/undefined.h>
-#include <poincare/decimal.h>
-#include <poincare/addition.h>
-#include <poincare/multiplication.h>
-#include <poincare/symbol.h>
-#include <poincare/ieee754.h>
-#include <poincare/layout_engine.h>
-#include "layout/horizontal_layout.h"
-#include "layout/vertical_offset_layout.h"
-
-#include <cmath>
-#include <ion.h>
-
-extern "C" {
-#include <assert.h>
-#include <float.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-}
-
-namespace Poincare {
-
-template<typename T>
-Complex<T> Complex<T>::Float(T x) {
-  return Complex(x,0);
-}
-
-template<typename T>
-Complex<T> Complex<T>::Cartesian(T a, T b) {
-  return Complex(a,b);
-}
-
-template<typename T>
-Complex<T> Complex<T>::Polar(T r, T th)  {
-  // If the radius is 0, theta may be undefined but shouldn't be able to affect the result.
-  if (r == 0) {
-    return Complex(0,0);
-  }
-  T c = std::cos(th);
-  T s = std::sin(th);
-  /* Cheat: see comment on cosine.cpp.
-   * Sine and cosine openbsd immplementationd are numerical approximation.
-   * We though want to avoid evaluating e^(i*pi) to -1+1E-17i. We thus round
-   * cosine and sine results to 0 if they are negligible compared to the
-   * argument th. */
-  c = th != 0 && std::fabs(c/th) <= Expression::epsilon<T>() ? 0 : c;
-  s = th != 0 && std::fabs(s/th) <= Expression::epsilon<T>() ? 0 : s;
-  return Complex(r*c,r*s);
-}
-
-template<typename T>
-Complex<T>::Complex(const Complex<T> & other) {
-  m_a = other.m_a;
-  m_b = other.m_b;
-}
-
-template<typename T>
-Complex<T> & Complex<T>::operator=(const Complex& other) {
-  m_a = other.m_a;
-  m_b = other.m_b;
-  return *this;
-}
-
-template<typename T>
-static inline T privateFloatSetSign(T f, bool negative) {
-  if (negative) {
-    return -f;
-  }
-  return f;
-}
-
-template<typename T>
-T digitsToFloat(const char * digits, int length) {
-  if (digits == nullptr) {
-    return 0;
-  }
-  T result = 0;
-  const char * digit = digits;
-  for (int i = 0; i < length; i++) {
-    result = 10 * result;
-    result += *digit-'0';
-    digit++;
-  }
-  return result;
-}
-
-template<typename T>
-Complex<T>::Complex(const char * integralPart, int integralPartLength, bool integralNegative,
-    const char * fractionalPart, int fractionalPartLength,
-    const char * exponent, int exponentLength, bool exponentNegative) {
-  T i = digitsToFloat<T>(integralPart, integralPartLength);
-  T j = digitsToFloat<T>(fractionalPart, fractionalPartLength);
-  T l = privateFloatSetSign<T>(digitsToFloat<T>(exponent, exponentLength), exponentNegative);
-
-  m_a = privateFloatSetSign((i + j*std::pow(10, -std::ceil((T)fractionalPartLength)))* std::pow(10, l), integralNegative);
-  m_b = 0;
-}
-
-template <class T>
-T Complex<T>::a() const {
-  return m_a;
-}
-
-template <class T>
-T Complex<T>::b() const {
-  return m_b;
-}
-
-template <class T>
-T Complex<T>::r() const {
-  // We want to avoid a^2 and b^2 which could both easily overflow.
-  // min, max = minmax(abs(a), abs(b)) (*minmax returns both arguments sorted*)
-  // abs(a + bi) == sqrt(a^2 + b^2)
-  //             == sqrt(abs(a)^2 + abs(b)^2)
-  //             == sqrt(min^2 + max^2)
-  //             == sqrt((min^2 + max^2) * max^2/max^2)
-  //             == sqrt((min^2 + max^2) / max^2)*sqrt(max^2)
-  //             == sqrt(min^2/max^2 + 1) * max
-  //             == sqrt((min/max)^2 + 1) * max
-  // min >= 0 &&
-  // max >= 0 &&
-  // min <= max => min/max <= 1
-  //            => (min/max)^2 <= 1
-  //            => (min/max)^2 + 1 <= 2
-  //            => sqrt((min/max)^2 + 1) <= sqrt(2)
-  // So the calculation is guaranteed to not overflow until the final multiply.
-  // If (min/max)^2 underflows then min doesn't contribute anything significant
-  // compared to max, and the formula reduces to simply max as it should.
-  // We do need to be careful about the case where a == 0 && b == 0 which would
-  // cause a division by zero.
-  T min = std::fabs(m_a);
-  if (m_b == 0) {
-    return min;
-  }
-  T max = std::fabs(m_b);
-  if (max < min) {
-    T temp = min;
-    min = max;
-    max = temp;
-  }
-  T temp = min/max;
-  return std::sqrt(temp*temp + 1) * max;
-}
-
-template <class T>
-T Complex<T>::th() const {
-  T result = std::atan(m_b/m_a) + M_PI;
-  if (m_a >= 0) {
-    T a = m_a == 0 ? 0 : m_a;
-    result = std::atan(m_b/a);
-  }
-  if (result > M_PI + FLT_EPSILON) {
-    result = result - 2*M_PI;
-  }
-  return result;
-}
-
-template <class T>
-Complex<T> Complex<T>::conjugate() const {
-  return Cartesian(m_a, -m_b);
-}
-
-template<typename T>
-Expression::Type Complex<T>::type() const {
-  return Expression::Type::Complex;
-}
-
-template <class T>
-Complex<T> * Complex<T>::clone() const {
-  return new Complex<T>(*this);
-}
-
-template<typename T>
-T Complex<T>::toScalar() const {
-  if (m_b != 0) {
-    return NAN;
-  }
-  return m_a;
-}
-
-template <class T>
-int Complex<T>::writeTextInBuffer(char * buffer, int bufferSize, int numberOfSignificantDigits) const {
-  return convertComplexToText(buffer, bufferSize, numberOfSignificantDigits, Preferences::sharedPreferences()->displayMode(), Preferences::sharedPreferences()->complexFormat(), Ion::Charset::MultiplicationSign);
-}
-
-template <class T>
-bool Complex<T>::needParenthesisWithParent(const Expression * e) const {
-  switch (e->type()) {
-    case Type::Addition:
-      return m_a < 0.0 || (m_a == 0.0 && m_b < 0.0);
-    case Type::Subtraction:
-    case Type::Multiplication:
-    case Type::Opposite:
-      return m_a < 0.0 || m_b < 0.0 || (m_a != 0.0 && m_b != 0.0);
-    case Type::Factorial:
-    case Type::Power:
-    case Type::Division:
-      return m_a < 0.0 || m_b != 0.0;
-    default:
-      return false;
-  }
-}
-
-
-
-template <class T>
-Complex<T>::Complex(T a, T b) :
-  m_a(a),
-  m_b(b)
-{
-}
-
-template <class T>
-ExpressionLayout * Complex<T>::privateCreateLayout(PrintFloat::Mode floatDisplayMode, Expression::ComplexFormat complexFormat) const {
-  assert(floatDisplayMode != PrintFloat::Mode::Default);
-  if (complexFormat == Expression::ComplexFormat::Polar) {
-    return createPolarLayout(floatDisplayMode);
-  }
-  return createCartesianLayout(floatDisplayMode);
-}
-
-template <class T>
-Expression * Complex<T>::CreateDecimal(T f) {
-  if (std::isnan(f) || std::isinf(f)) {
-    return new Undefined();
-  }
-  int e = IEEE754<T>::exponentBase10(f);
-  int64_t mantissaf = f * std::pow((T)10, -e+PrintFloat::k_numberOfStoredSignificantDigits+1);
-  return new Decimal(Integer(mantissaf), e);
-}
-
-template <class T>
-Expression * Complex<T>::shallowReduce(Context & context, AngleUnit angleUnit) {
-  Expression * a = CreateDecimal(m_a);
-  Expression * b = CreateDecimal(m_b);
-  Multiplication * m = new Multiplication(new Symbol(Ion::Charset::IComplex), b, false);
-  Addition * add = new Addition(a, m, false);
-  a->shallowReduce(context, angleUnit);
-  b->shallowReduce(context, angleUnit);
-  m->shallowReduce(context, angleUnit);
-  return replaceWith(add, true)->shallowReduce(context, angleUnit);
-}
-
-template<typename T>
-template<typename U>
-Complex<U> * Complex<T>::templatedApproximate(Context& context, Expression::AngleUnit angleUnit) const {
-  return new Complex<U>(Complex<U>::Cartesian((U)m_a, (U)m_b));
-}
-
-template <class T>
-int Complex<T>::convertComplexToText(char * buffer, int bufferSize, int numberOfSignificantDigits, PrintFloat::Mode displayMode, Expression::ComplexFormat complexFormat, char multiplicationSpecialChar) const {
-  assert(displayMode != PrintFloat::Mode::Default);
-  int numberOfChars = 0;
-  if (std::isnan(m_a) || std::isnan(m_b) || std::isinf(m_a) || std::isinf(m_b)) {
-    return PrintFloat::convertFloatToText<T>(NAN, buffer, bufferSize, numberOfSignificantDigits, displayMode);
-  }
-  if (complexFormat == Expression::ComplexFormat::Polar) {
-    if (r() != 1 || th() == 0) {
-      numberOfChars = PrintFloat::convertFloatToText<T>(r(), buffer, bufferSize, numberOfSignificantDigits, displayMode);
-      if (r() != 0 && th() != 0 && bufferSize > numberOfChars+1) {
-        buffer[numberOfChars++] = multiplicationSpecialChar;
-        // Ensure that the string is null terminated even if buffer size is to small
-        buffer[numberOfChars] = 0;
-      }
-    }
-    if (r() != 0 && th() != 0) {
-      if (bufferSize > numberOfChars+3) {
-        buffer[numberOfChars++] = Ion::Charset::Exponential;
-        buffer[numberOfChars++] = '^';
-        buffer[numberOfChars++] = '(';
-        // Ensure that the string is null terminated even if buffer size is to small
-        buffer[numberOfChars] = 0;
-      }
-      numberOfChars += PrintFloat::convertFloatToText<T>(th(), buffer+numberOfChars, bufferSize-numberOfChars, numberOfSignificantDigits, displayMode);
-      if (bufferSize > numberOfChars+3) {
-        buffer[numberOfChars++] = multiplicationSpecialChar;
-        buffer[numberOfChars++] = Ion::Charset::IComplex;
-        buffer[numberOfChars++] = ')';
-        buffer[numberOfChars] = 0;
-      }
-    }
-    return numberOfChars;
-  }
-
-  if (m_a != 0 || m_b == 0) {
-    numberOfChars = PrintFloat::convertFloatToText<T>(m_a, buffer, bufferSize, numberOfSignificantDigits, displayMode);
-    if (m_b > 0 && !std::isnan(m_b) && bufferSize > numberOfChars+1) {
-      buffer[numberOfChars++] = '+';
-      // Ensure that the string is null terminated even if buffer size is to small
-      buffer[numberOfChars] = 0;
-    }
-  }
-  if (m_b != 1 && m_b != -1 && m_b != 0) {
-    numberOfChars += PrintFloat::convertFloatToText<T>(m_b, buffer+numberOfChars, bufferSize-numberOfChars, numberOfSignificantDigits, displayMode);
-    buffer[numberOfChars++] = multiplicationSpecialChar;
-  }
-  if (m_b == -1 && bufferSize > numberOfChars+1) {
-    buffer[numberOfChars++] = '-';
-  }
-  if (m_b != 0 && bufferSize > numberOfChars+1) {
-    buffer[numberOfChars++] = Ion::Charset::IComplex;
-    buffer[numberOfChars] = 0;
-  }
-  return numberOfChars;
-}
-
-template <class T>
-ExpressionLayout * Complex<T>::createPolarLayout(PrintFloat::Mode floatDisplayMode) const {
-  char bufferBase[PrintFloat::k_maxFloatBufferLength+2];
-  int numberOfCharInBase = 0;
-  char bufferSuperscript[PrintFloat::k_maxFloatBufferLength+2];
-  int numberOfCharInSuperscript = 0;
-
-  if (std::isnan(r()) || (std::isnan(th()) && r() != 0)) {
-    numberOfCharInBase = PrintFloat::convertFloatToText<T>(NAN, bufferBase, PrintFloat::k_maxComplexBufferLength, Preferences::sharedPreferences()->numberOfSignificantDigits(), floatDisplayMode);
-    return LayoutEngine::createStringLayout(bufferBase, numberOfCharInBase);
-  }
-  if (r() != 1 || th() == 0) {
-    numberOfCharInBase = PrintFloat::convertFloatToText<T>(r(), bufferBase, PrintFloat::k_maxFloatBufferLength, Preferences::sharedPreferences()->numberOfSignificantDigits(), floatDisplayMode);
-    if (r() != 0 && th() != 0) {
-      bufferBase[numberOfCharInBase++] = Ion::Charset::MiddleDot;
-    }
-  }
-  if (r() != 0 && th() != 0) {
-    bufferBase[numberOfCharInBase++] = Ion::Charset::Exponential;
-    bufferBase[numberOfCharInBase] = 0;
-  }
-
-  if (r() != 0 && th() != 0) {
-    numberOfCharInSuperscript = PrintFloat::convertFloatToText<T>(th(), bufferSuperscript, PrintFloat::k_maxFloatBufferLength, Preferences::sharedPreferences()->numberOfSignificantDigits(), floatDisplayMode);
-    bufferSuperscript[numberOfCharInSuperscript++] = Ion::Charset::MiddleDot;
-    bufferSuperscript[numberOfCharInSuperscript++] = Ion::Charset::IComplex;
-    bufferSuperscript[numberOfCharInSuperscript] = 0;
-  }
-  if (numberOfCharInSuperscript == 0) {
-    return LayoutEngine::createStringLayout(bufferBase, numberOfCharInBase);
-  }
-  ExpressionLayout * result = LayoutEngine::createStringLayout(bufferBase, numberOfCharInBase);
-  ExpressionLayout * exponentLayout = LayoutEngine::createStringLayout(bufferSuperscript, numberOfCharInSuperscript);
-  VerticalOffsetLayout * offsetLayout = new VerticalOffsetLayout(exponentLayout, VerticalOffsetLayout::Type::Superscript, false);
-  (static_cast<HorizontalLayout *>(result))->addChildAtIndex(offsetLayout, result->numberOfChildren());
-  return result;
-}
-
-template <class T>
-ExpressionLayout * Complex<T>::createCartesianLayout(PrintFloat::Mode floatDisplayMode) const {
-  char buffer[PrintFloat::k_maxComplexBufferLength];
-  int numberOfChars = convertComplexToText(buffer, PrintFloat::k_maxComplexBufferLength, Preferences::sharedPreferences()->numberOfSignificantDigits(), floatDisplayMode, Expression::ComplexFormat::Cartesian, Ion::Charset::MiddleDot);
-  return LayoutEngine::createStringLayout(buffer, numberOfChars);
-}
-
-template class Complex<float>;
-template class Complex<double>;
-template Complex<double>* Complex<double>::templatedApproximate<double>(Context&, Expression::AngleUnit) const;
-template Complex<float>* Complex<double>::templatedApproximate<float>(Context&, Expression::AngleUnit) const;
-template Complex<double>* Complex<float>::templatedApproximate<double>(Context&, Expression::AngleUnit) const;
-template Complex<float>* Complex<float>::templatedApproximate<float>(Context&, Expression::AngleUnit) const;
-
-}
-
--- a/poincare/src/complex_argument.cpp
+++ b/poincare/src/complex_argument.cpp
@@ -1,5 +1,4 @@
 #include <poincare/complex_argument.h>
-#include <poincare/complex.h>
 #include <poincare/simplification_engine.h>
 extern "C" {
 #include <assert.h>
@@ -32,8 +31,8 @@ Expression * ComplexArgument::shallowReduce(Context& context, AngleUnit angleUni
 }

 template<typename T>
-Complex<T> ComplexArgument::computeOnComplex(const Complex<T> c, AngleUnit angleUnit) {
-  return Complex<T>::Float(c.th());
+std::complex<T> ComplexArgument::computeOnComplex(const std::complex<T> c, AngleUnit angleUnit) {
+  return Complex<T>(std::arg(c));
 }

 }
--- a/poincare/src/confidence_interval.cpp
+++ b/poincare/src/confidence_interval.cpp
@@ -1,6 +1,5 @@
 #include <poincare/confidence_interval.h>
 #include <poincare/matrix.h>
-#include <poincare/complex.h>
 #include <poincare/addition.h>
 #include <poincare/multiplication.h>
 #include <poincare/power.h>
@@ -63,23 +62,20 @@ Expression * ConfidenceInterval::shallowReduce(Context& context, AngleUnit angle
 }

 template<typename T>
-Expression * ConfidenceInterval::templatedApproximate(Context& context, AngleUnit angleUnit) const {
-  Expression * fInput = operand(0)->approximate<T>(context, angleUnit);
-  Expression * nInput = operand(1)->approximate<T>(context, angleUnit);
-  if (fInput->type() != Type::Complex || nInput->type() != Type::Complex) {
-    return new Complex<T>(Complex<T>::Float(NAN));
-  }
+Evaluation<T> * ConfidenceInterval::templatedApproximate(Context& context, AngleUnit angleUnit) const {
+  Evaluation<T> * fInput = operand(0)->privateApproximate(T(), context, angleUnit);
+  Evaluation<T> * nInput = operand(1)->privateApproximate(T(), context, angleUnit);
  T f = static_cast<Complex<T> *>(fInput)->toScalar();
  T n = static_cast<Complex<T> *>(nInput)->toScalar();
  delete fInput;
  delete nInput;
  if (std::isnan(f) || std::isnan(n) || n != (int)n || n < 0 || f < 0 || f > 1) {
-    return new Complex<T>(Complex<T>::Float(NAN));
+    return new Complex<T>(Complex<T>::Undefined());
  }
-  Expression * operands[2];
-  operands[0] = new Complex<T>(Complex<T>::Float(f - 1/std::sqrt(n)));
-  operands[1] = new Complex<T>(Complex<T>::Float(f + 1/std::sqrt(n)));
-  return new Matrix(operands, 1, 2, false);
+  std::complex<T> operands[2];
+  operands[0] = std::complex<T>(f - 1/std::sqrt(n));
+  operands[1] = std::complex<T>(f + 1/std::sqrt(n));
+  return new MatrixComplex<T>(operands, 1, 2);
 }

 }
--- a/poincare/src/conjugate.cpp
+++ b/poincare/src/conjugate.cpp
@@ -1,5 +1,4 @@
 #include <poincare/conjugate.h>
-#include <poincare/complex.h>
 #include <poincare/simplification_engine.h>
 #include "layout/conjugate_layout.h"

@@ -43,8 +42,8 @@ Expression * Conjugate::shallowReduce(Context& context, AngleUnit angleUnit) {
 }

 template<typename T>
-Complex<T> Conjugate::computeOnComplex(const Complex<T> c, AngleUnit angleUnit) {
-  return c.conjugate();
+std::complex<T> Conjugate::computeOnComplex(const std::complex<T> c, AngleUnit angleUnit) {
+  return std::conj(c);
 }

 }
--- a/poincare/src/cosine.cpp
+++ b/poincare/src/cosine.cpp
@@ -1,6 +1,5 @@
 #include <poincare/cosine.h>
 #include <poincare/hyperbolic_cosine.h>
-#include <poincare/complex.h>
 #include <poincare/symbol.h>
 #include <poincare/rational.h>
 #include <poincare/multiplication.h>
@@ -40,29 +39,4 @@ Expression * Cosine::shallowReduce(Context& context, AngleUnit angleUnit) {
  return Trigonometry::shallowReduceDirectFunction(this, context, angleUnit);
 }

-template<typename T>
-Complex<T> Cosine::computeOnComplex(const Complex<T> c, AngleUnit angleUnit) {
-  assert(angleUnit != AngleUnit::Default);
-  if (c.b() == 0) {
-    T input = c.a();
-    if (angleUnit == AngleUnit::Degree) {
-      input *= M_PI/180.0f;
-    }
-    T result = std::cos(input);
-    /* Cheat: openbsd trigonometric functions (cos, sin & tan) are numerical
-     * implementation and thus are approximative. The error epsilon is ~1E-7
-     * on float and ~1E-15 on double. In order to avoid weird results as
-     * cos(90) = 6E-17, we neglect the result when its ratio with the argument
-     * (pi in the exemple) is smaller than epsilon.
-     * We can't do that for all evaluation as the user can operate on values as
-     * small as 1E-308 (in double) and most results still be correct. */
-    if (input !=  0 && std::fabs(result/input) <= epsilon<T>()) {
-      return Complex<T>::Float(0);
-    }
-    return Complex<T>::Float(result);
-  }
-  Complex<T> arg = Complex<T>::Cartesian(-c.b(), c.a());
-  return HyperbolicCosine::computeOnComplex(arg, angleUnit);
-}
-
 }
--- a/poincare/src/decimal.cpp
+++ b/poincare/src/decimal.cpp
@@ -1,5 +1,4 @@
 #include <poincare/decimal.h>
-#include <poincare/complex.h>
 #include <poincare/rational.h>
 #include <poincare/opposite.h>
 #include <poincare/ieee754.h>
@@ -90,10 +89,10 @@ Expression * Decimal::clone() const {
  return new Decimal(m_mantissa, m_exponent);
 }

-template<typename T> Expression * Decimal::templatedApproximate(Context& context, Expression::AngleUnit angleUnit) const {
+template<typename T> Evaluation<T> * Decimal::templatedApproximate(Context& context, Expression::AngleUnit angleUnit) const {
  T m = m_mantissa.approximate<T>();
  int numberOfDigits = Integer::numberOfDigitsWithoutSign(m_mantissa);
-  return new Complex<T>(Complex<T>::Float(m*std::pow((T)10.0, (T)(m_exponent-numberOfDigits+1))));
+  return new Complex<T>(m*std::pow((T)10.0, (T)(m_exponent-numberOfDigits+1)));
 }

 int Decimal::convertToText(char * buffer, int bufferSize, PrintFloat::Mode mode, int numberOfSignificantDigits) const {
--- a/poincare/src/derivative.cpp
+++ b/poincare/src/derivative.cpp
@@ -1,6 +1,5 @@
 #include <poincare/derivative.h>
 #include <poincare/symbol.h>
-#include <poincare/complex.h>
 #include <poincare/simplification_engine.h>
 #include <poincare/undefined.h>
 #include <cmath>
@@ -45,16 +44,16 @@ Expression * Derivative::shallowReduce(Context& context, AngleUnit angleUnit) {
 }

 template<typename T>
-Expression * Derivative::templatedApproximate(Context& context, AngleUnit angleUnit) const {
+Complex<T> * Derivative::templatedApproximate(Context& context, AngleUnit angleUnit) const {
  static T min = sizeof(T) == sizeof(double) ? DBL_MIN : FLT_MIN;
  static T epsilon = sizeof(T) == sizeof(double) ? DBL_EPSILON : FLT_EPSILON;
-  Symbol xSymbol('x');
-  T x = operand(1)->approximateToScalar<T>(context, angleUnit);
+  Evaluation<T> * xInput = operand(1)->privateApproximate(T(), context, angleUnit);
+  T x = xInput->toScalar();
+  delete xInput;
  T functionValue = operand(0)->approximateWithValueForSymbol('x', x, context);
-
  // No complex/matrix version of Derivative
  if (std::isnan(x) || std::isnan(functionValue)) {
-    return new Complex<T>(Complex<T>::Float(NAN));
+    return new Complex<T>(Complex<T>::Undefined());
  }

  T error, result;
@@ -66,13 +65,13 @@ Expression * Derivative::templatedApproximate(Context& context, AngleUnit angleU

  /* if the error is too big regarding the value, do not return the answer */
  if (std::fabs(error/result) > k_maxErrorRateOnApproximation || std::isnan(error)) {
-    return new Complex<T>(Complex<T>::Float(NAN));
+    return new Complex<T>(Complex<T>::Undefined());
  }
  if (std::fabs(error) < min) {
-    return new Complex<T>(Complex<T>::Float(result));
+    return new Complex<T>(result);
  }
  error = std::pow((T)10, std::floor(std::log10(std::fabs(error)))+2);
-  return new Complex<T>(Complex<T>::Float(std::round(result/error)*error));
+  return new Complex<T>(std::round(result/error)*error);
 }

 template<typename T>
--- a/poincare/src/determinant.cpp
+++ b/poincare/src/determinant.cpp
@@ -34,15 +34,9 @@ Expression * Determinant::shallowReduce(Context& context, AngleUnit angleUnit) {

 // TODO: handle this exactly in shallowReduce for small dimensions.
 template<typename T>
-Expression * Determinant::templatedApproximate(Context& context, AngleUnit angleUnit) const {
-  Expression * input = operand(0)->approximate<T>(context, angleUnit);
-  Expression * result = nullptr;
-  if (input->type() == Type::Complex) {
-    result = input->clone();
-  } else {
-    assert(input->type() == Type::Matrix);
-    result = static_cast<Matrix *>(input)->createDeterminant<T>(context, angleUnit);
-  }
+Evaluation<T> * Determinant::templatedApproximate(Context& context, AngleUnit angleUnit) const {
+  Evaluation<T> * input = operand(0)->privateApproximate(T(), context, angleUnit);
+  Complex<T> * result = new Complex<T>(input->createDeterminant());
  delete input;
  return result;
 }
--- a/poincare/src/division.cpp
+++ b/poincare/src/division.cpp
@@ -49,52 +49,8 @@ Expression * Division::shallowReduce(Context& context, AngleUnit angleUnit) {
 }

 template<typename T>
-Complex<T> Division::compute(const Complex<T> c, const Complex<T> d) {
-  /* We want to avoid multiplies in the middle of the calculation that could
-   * overflow.
-   * aa, ab, ba, bb, min, max = |d.a| <= |d.b| ? (c.a, c.b, -c.a, c.b, d.a, d.b)
-   *    : (c.b, c.a, c.b, -c.a, d.b, d.a)
-   * c    c.a+c.b*i   d.a-d.b*i   1/max    (c.a+c.b*i) * (d.a-d.b*i) / max
-   * - == --------- * --------- * ----- == -------------------------------
-   * d    d.a+d.b*i   d.a-d.b*i   1/max    (d.a+d.b*i) * (d.a-d.b*i) / max
-   *      (c.a*d.a - c.a*d.b*i + c.b*i*d.a - c.b*i*d.b*i) / max
-   *   == -----------------------------------------------------
-   *      (d.a*d.a - d.a*d.b*i + d.b*i*d.a - d.b*i*d.b*i) / max
-   *      (c.a*d.a - c.b*d.b*i^2 + c.b*d.b*i - c.a*d.a*i) / max
-   *   == -----------------------------------------------------
-   *                  (d.a*d.a - d.b*d.b*i^2) / max
-   *      (c.a*d.a/max + c.b*d.b/max) + (c.b*d.b/max - c.a*d.a/max)*i
-   *   == -----------------------------------------------------------
-   *                         d.a^2/max + d.b^2/max
-   *      aa*min/max + ab*max/max   bb*min/max + ba*max/max
-   *   == ----------------------- + -----------------------*i
-   *       min^2/max + max^2/max     min^2/max + max^2/max
-   *       min/max*aa + ab     min/max*bb + ba
-   *   == ----------------- + -----------------*i
-   *      min/max*min + max   min/max*min + max
-   * |min| <= |max| => |min/max| <= 1
-   *                => |min/max*x| <= |x|
-   *                => |min/max*x+y| <= |x|+|y|
-   * So the calculation is guaranteed to not overflow until the last divides as
-   * long as none of the input values have the representation's maximum exponent.
-   * Plus, the method does not propagate any error on real inputs: temp = 0,
-   * norm = d.a and then result = c.a/d.a. */
-  T aa = c.a(), ab = c.b(), ba = -aa, bb = ab;
-  T min = d.a(), max = d.b();
-  if (std::fabs(max) < std::fabs(min)) {
-    T temp = min;
-    min = max;
-    max = temp;
-    temp = aa;
-    aa = ab;
-    ab = temp;
-    temp = ba;
-    ba = bb;
-    bb = temp;
-  }
-  T temp = min/max;
-  T norm = temp*min + max;
-  return Complex<T>::Cartesian((temp*aa + ab) / norm, (temp*bb + ba) / norm);
+std::complex<T> Division::compute(const std::complex<T> c, const std::complex<T> d) {
+  return c/d;
 }

 ExpressionLayout * Division::privateCreateLayout(PrintFloat::Mode floatDisplayMode, ComplexFormat complexFormat) const {
@@ -105,25 +61,25 @@ ExpressionLayout * Division::privateCreateLayout(PrintFloat::Mode floatDisplayMo
  return new FractionLayout(numerator->createLayout(floatDisplayMode, complexFormat), denominator->createLayout(floatDisplayMode, complexFormat), false);
 }

-template<typename T> Matrix * Division::computeOnComplexAndMatrix(const Complex<T> * c, const Matrix * n) {
-  Matrix * inverse = n->createApproximateInverse<T>();
+template<typename T> MatrixComplex<T> Division::computeOnComplexAndMatrix(const std::complex<T> c, const MatrixComplex<T> n) {
+  MatrixComplex<T> * inverse = n.createInverse();
  if (inverse == nullptr) {
-    return nullptr;
+    return MatrixComplex<T>::Undefined();
  }
-  Matrix * result = Multiplication::computeOnComplexAndMatrix<T>(c, inverse);
+  MatrixComplex<T> result = Multiplication::computeOnComplexAndMatrix<T>(c, *inverse);
  delete inverse;
  return result;
 }

-template<typename T> Matrix * Division::computeOnMatrices(const Matrix * m, const Matrix * n) {
-  if (m->numberOfColumns() != n->numberOfColumns()) {
-    return nullptr;
+template<typename T> MatrixComplex<T> Division::computeOnMatrices(const MatrixComplex<T> m, const MatrixComplex<T> n) {
+  if (m.numberOfColumns() != n.numberOfColumns()) {
+    return MatrixComplex<T>::Undefined();
  }
-  Matrix * inverse = n->createApproximateInverse<T>();
+  MatrixComplex<T> * inverse = n.createInverse();
  if (inverse == nullptr) {
-    return nullptr;
+    return MatrixComplex<T>::Undefined();
  }
-  Matrix * result = Multiplication::computeOnMatrices<T>(m, inverse);
+  MatrixComplex<T> result = Multiplication::computeOnMatrices<T>(m, *inverse);
  delete inverse;
  return result;
 }
--- a/poincare/src/division_quotient.cpp
+++ b/poincare/src/division_quotient.cpp
@@ -59,16 +59,16 @@ Expression * DivisionQuotient::shallowReduce(Context& context, AngleUnit angleUn

 template<typename T>
 Complex<T> * DivisionQuotient::templatedApproximate(Context& context, AngleUnit angleUnit) const {
-  Expression * f1Input = operand(0)->approximate<T>(context, angleUnit);
-  Expression * f2Input = operand(1)->approximate<T>(context, angleUnit);
-  T f1 = f1Input->type() == Type::Complex ? static_cast<Complex<T> *>(f1Input)->toScalar() : NAN;
-  T f2 = f2Input->type() == Type::Complex ? static_cast<Complex<T> *>(f2Input)->toScalar() : NAN;
+  Evaluation<T> * f1Input = operand(0)->privateApproximate(T(), context, angleUnit);
+  Evaluation<T> * f2Input = operand(1)->privateApproximate(T(), context, angleUnit);
+  T f1 = f1Input->toScalar();
+  T f2 = f2Input->toScalar();
  delete f1Input;
  delete f2Input;
  if (std::isnan(f1) || std::isnan(f2) || f1 != (int)f1 || f2 != (int)f2) {
-    return new Complex<T>(Complex<T>::Float(NAN));
+    return new Complex<T>(Complex<T>::Undefined());
  }
-  return new Complex<T>(Complex<T>::Float(std::floor(f1/f2)));
+  return new Complex<T>(std::floor(f1/f2));
 }

 }
--- a/Show More
+++ b/Show More