diff --git a/apps/calculation/calculation.cpp b/apps/calculation/calculation.cpp
index a38f4ff36..192d7b41b 100644
--- 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;
diff --git a/apps/graph/cartesian_function.cpp b/apps/graph/cartesian_function.cpp
index 0112e2649..d3a5bc170 100644
--- 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. */
diff --git a/apps/regression/model/model.cpp b/apps/regression/model/model.cpp
index 46fb3ae88..8b9f55e17 100644
--- 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>
diff --git a/apps/regression/regression_context.cpp b/apps/regression/regression_context.cpp
index b1faf2eee..a36be05e3 100644
--- 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);
   }
diff --git a/apps/sequence/cache_context.cpp b/apps/sequence/cache_context.cpp
index f3683b5a9..29a9a1bf2 100644
--- 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*);
 
 }
diff --git a/apps/sequence/cache_context.h b/apps/sequence/cache_context.h
index c1c5536e0..f469fa20c 100644
--- 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];
 };
 
 }
diff --git a/apps/sequence/sequence.cpp b/apps/sequence/sequence.cpp
index 50f13ea2f..0a2087fa3 100644
--- 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);
diff --git a/apps/shared/store_context.h b/apps/shared/store_context.h
index 7284a84a5..4b53e0206 100644
--- 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
 };
 
 }
diff --git a/apps/statistics/statistics_context.cpp b/apps/statistics/statistics_context.cpp
index 5cb4ee356..c75f3ec3e 100644
--- 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);
   }
diff --git a/poincare/Makefile b/poincare/Makefile
index fa1b389fa..a7178cd2f 100644
--- 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\
diff --git a/poincare/include/poincare.h b/poincare/include/poincare.h
index 279682893..129a219eb 100644
--- 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>
diff --git a/poincare/include/poincare/absolute_value.h b/poincare/include/poincare/absolute_value.h
index c5b75d48d..71932af61 100644
--- 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>);
   }
 };
diff --git a/poincare/include/poincare/addition.h b/poincare/include/poincare/addition.h
index aa00d6313..76bdf357c 100644
--- 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>);
    }
 };
diff --git a/poincare/include/poincare/approximation_engine.h b/poincare/include/poincare/approximation_engine.h
index 673aa5226..33f40765d 100644
--- 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);
 
 };
 
diff --git a/poincare/include/poincare/arc_cosine.h b/poincare/include/poincare/arc_cosine.h
index 71ae77674..62fc3b053 100644
--- 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>);
   }
 };
diff --git a/poincare/include/poincare/arc_sine.h b/poincare/include/poincare/arc_sine.h
index 3aa871d05..2d2a823ad 100644
--- 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>);
   }
 };
diff --git a/poincare/include/poincare/arc_tangent.h b/poincare/include/poincare/arc_tangent.h
index e915ab332..26cbd5756 100644
--- 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>);
   }
 };
diff --git a/poincare/include/poincare/binomial_coefficient.h b/poincare/include/poincare/binomial_coefficient.h
index b5085bb67..a3c64d792 100644
--- 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;
 };
 
 }
diff --git a/poincare/include/poincare/ceiling.h b/poincare/include/poincare/ceiling.h
index e55e96f8c..cd0fc9926 100644
--- 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>);
   }
 };
diff --git a/poincare/include/poincare/complex.h b/poincare/include/poincare/complex.h
deleted file mode 100644
index 47545eaa3..000000000
--- a/poincare/include/poincare/complex.h
+++ /dev/null
@@ -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
diff --git a/poincare/include/poincare/complex_argument.h b/poincare/include/poincare/complex_argument.h
index 14d283528..0cfcc5738 100644
--- 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>);
   }
 };
diff --git a/poincare/include/poincare/confidence_interval.h b/poincare/include/poincare/confidence_interval.h
index f4e815703..ea95aa054 100644
--- 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 {
diff --git a/poincare/include/poincare/conjugate.h b/poincare/include/poincare/conjugate.h
index e43767171..10a0fa5b3 100644
--- 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>);
   }
 };
diff --git a/poincare/include/poincare/cosine.h b/poincare/include/poincare/cosine.h
index b19ec7bad..942f811f4 100644
--- 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>);
   }
 };
diff --git a/poincare/include/poincare/decimal.h b/poincare/include/poincare/decimal.h
index 33892903f..e6747ffdc 100644
--- 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
diff --git a/poincare/include/poincare/derivative.h b/poincare/include/poincare/derivative.h
index 85e10c3a4..11bd0fad1 100644
--- 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?
diff --git a/poincare/include/poincare/determinant.h b/poincare/include/poincare/determinant.h
index 7db860bc1..d59dbd111 100644
--- 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;
 };
 
 }
diff --git a/poincare/include/poincare/division.h b/poincare/include/poincare/division.h
index 0337ec310..515496540 100644
--- 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>);
   }
 };
diff --git a/poincare/include/poincare/division_quotient.h b/poincare/include/poincare/division_quotient.h
index 318d09c0e..02b96ca61 100644
--- 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;
 };
 
diff --git a/poincare/include/poincare/division_remainder.h b/poincare/include/poincare/division_remainder.h
index 05a79ea73..5f00ab180 100644
--- 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;
 
 };
diff --git a/poincare/include/poincare/empty_expression.h b/poincare/include/poincare/empty_expression.h
index 09d1cbb9d..ad3a3c926 100644
--- 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;
 };
 
diff --git a/poincare/include/poincare/equal.h b/poincare/include/poincare/equal.h
index 14d775f23..062d527d2 100644
--- 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;
 };
 
 }
diff --git a/poincare/include/poincare/evaluation.h b/poincare/include/poincare/evaluation.h
new file mode 100644
index 000000000..83e94d1c9
--- /dev/null
+++ 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
diff --git a/poincare/include/poincare/evaluation_mode.h b/poincare/include/poincare/evaluation_mode.h
new file mode 100644
index 000000000..fea357ee5
--- /dev/null
+++ 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
+};
+
+}
diff --git a/poincare/include/poincare/expression.h b/poincare/include/poincare/expression.h
index cf77d1cf0..4c2a52459 100644
--- 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;
diff --git a/poincare/include/poincare/factor.h b/poincare/include/poincare/factor.h
index 4ded7f021..bcd10db71 100644
--- 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);
   }
 };
 
diff --git a/poincare/include/poincare/factorial.h b/poincare/include/poincare/factorial.h
index 35c04b3fe..57a0ee3c4 100644
--- 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>);
   }
 
diff --git a/poincare/include/poincare/floor.h b/poincare/include/poincare/floor.h
index 8a4196f32..d75200aad 100644
--- 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>);
   }
 };
diff --git a/poincare/include/poincare/frac_part.h b/poincare/include/poincare/frac_part.h
index 2fb4f6d88..55e49c757 100644
--- 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>);
   }
 };
diff --git a/poincare/include/poincare/global_context.h b/poincare/include/poincare/global_context.h
index d6c4d480a..51de1a005 100644
--- 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;
 };
 
 }
diff --git a/poincare/include/poincare/great_common_divisor.h b/poincare/include/poincare/great_common_divisor.h
index 83d8f5d9e..cacf2b731 100644
--- 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;
 };
 
diff --git a/poincare/include/poincare/hyperbolic_arc_cosine.h b/poincare/include/poincare/hyperbolic_arc_cosine.h
index 8bc94757e..b222d871d 100644
--- 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>);
   }
 };
diff --git a/poincare/include/poincare/hyperbolic_arc_sine.h b/poincare/include/poincare/hyperbolic_arc_sine.h
index 41f7ac7bc..bd60622b7 100644
--- 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>);
   }
 };
diff --git a/poincare/include/poincare/hyperbolic_arc_tangent.h b/poincare/include/poincare/hyperbolic_arc_tangent.h
index 1beaa2947..d5cf6e2c2 100644
--- 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>);
   }
 };
diff --git a/poincare/include/poincare/hyperbolic_cosine.h b/poincare/include/poincare/hyperbolic_cosine.h
index a3325b86f..176246644 100644
--- 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>);
   }
 };
diff --git a/poincare/include/poincare/hyperbolic_sine.h b/poincare/include/poincare/hyperbolic_sine.h
index 83b4c5f2f..fadfaed5c 100644
--- 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>);
   }
 };
diff --git a/poincare/include/poincare/hyperbolic_tangent.h b/poincare/include/poincare/hyperbolic_tangent.h
index 4ec0dabed..4ff6408aa 100644
--- 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>);
   }
 };
diff --git a/poincare/include/poincare/imaginary_part.h b/poincare/include/poincare/imaginary_part.h
index 2c5757b7b..fa663e83e 100644
--- 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>);
   }
 };
diff --git a/poincare/include/poincare/integral.h b/poincare/include/poincare/integral.h
index 2c3ab3f98..49ff7955a 100644
--- 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
diff --git a/poincare/include/poincare/least_common_multiple.h b/poincare/include/poincare/least_common_multiple.h
index eddf83d93..31d700a0a 100644
--- 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;
 };
 
diff --git a/poincare/include/poincare/logarithm.h b/poincare/include/poincare/logarithm.h
index ded4cce79..0592db84b 100644
--- 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;
 };
 
 }
diff --git a/poincare/include/poincare/matrix.h b/poincare/include/poincare/matrix.h
index 90a637d4b..80853de1d 100644
--- 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;
 };
 
diff --git a/poincare/include/poincare/matrix_dimension.h b/poincare/include/poincare/matrix_dimension.h
index f6b608b4c..404cd138b 100644
--- 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;
 };
 
 }
diff --git a/poincare/include/poincare/matrix_inverse.h b/poincare/include/poincare/matrix_inverse.h
index 0acae6ef6..1a4e3ed24 100644
--- 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;
 };
 
 }
diff --git a/poincare/include/poincare/matrix_trace.h b/poincare/include/poincare/matrix_trace.h
index 14b140ad6..a5d4b8574 100644
--- 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;
 };
 
 }
diff --git a/poincare/include/poincare/matrix_transpose.h b/poincare/include/poincare/matrix_transpose.h
index 42ad6455a..e80af6ffc 100644
--- 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;
 };
 
 }
diff --git a/poincare/include/poincare/multiplication.h b/poincare/include/poincare/multiplication.h
index daf0a4090..c37318f0d 100644
--- 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>);
   }
 };
diff --git a/poincare/include/poincare/naperian_logarithm.h b/poincare/include/poincare/naperian_logarithm.h
index 89a5173d1..3e24ba86e 100644
--- 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>);
   }
 };
diff --git a/poincare/include/poincare/nth_root.h b/poincare/include/poincare/nth_root.h
index dfd8f5392..35c5d4501 100644
--- 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;
 
 };
 
diff --git a/poincare/include/poincare/opposite.h b/poincare/include/poincare/opposite.h
index 809aa4c43..a5b5ad3cd 100644
--- 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>);
   }
 };
diff --git a/poincare/include/poincare/parenthesis.h b/poincare/include/poincare/parenthesis.h
index e3589c1d1..c8360ec1a 100644
--- 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;
 
 };
 
diff --git a/poincare/include/poincare/permute_coefficient.h b/poincare/include/poincare/permute_coefficient.h
index 63ad49a78..916b96215 100644
--- 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;
 };
 
diff --git a/poincare/include/poincare/power.h b/poincare/include/poincare/power.h
index a87b668dd..dbbeb9967 100644
--- 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>);
   }
 };
diff --git a/poincare/include/poincare/prediction_interval.h b/poincare/include/poincare/prediction_interval.h
index bccc77f8c..35d0b859d 100644
--- 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;
 };
 
 }
diff --git a/poincare/include/poincare/product.h b/poincare/include/poincare/product.h
index 954a1aa20..4eaab0d83 100644
--- 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;
 };
 
 }
diff --git a/poincare/include/poincare/randint.h b/poincare/include/poincare/randint.h
index 39d99e031..7bea9a70f 100644
--- 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;
 };
 
 }
diff --git a/poincare/include/poincare/random.h b/poincare/include/poincare/random.h
index 537b89106..355484a97 100644
--- 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>());
   }
 };
 
diff --git a/poincare/include/poincare/rational.h b/poincare/include/poincare/rational.h
index 075a43060..f7cf70770 100644
--- 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);
diff --git a/poincare/include/poincare/real_part.h b/poincare/include/poincare/real_part.h
index 16b6fdfd4..00fd1b746 100644
--- 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>);
   }
 };
diff --git a/poincare/include/poincare/round.h b/poincare/include/poincare/round.h
index fa0bd81da..9e50a7659 100644
--- 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;
 };
 
diff --git a/poincare/include/poincare/sequence.h b/poincare/include/poincare/sequence.h
index 6815ce4f4..e2394d689 100644
--- 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;
 };
 
 }
diff --git a/poincare/include/poincare/simplification_root.h b/poincare/include/poincare/simplification_root.h
index f692eda8c..d9ae4ef91 100644
--- 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;
   }
diff --git a/poincare/include/poincare/sine.h b/poincare/include/poincare/sine.h
index 0c4fa57e3..5f48208bd 100644
--- 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>);
   }
 };
diff --git a/poincare/include/poincare/square_root.h b/poincare/include/poincare/square_root.h
index 52b9f33f4..540ee7ccf 100644
--- 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>);
   }
 
diff --git a/poincare/include/poincare/store.h b/poincare/include/poincare/store.h
index f47b8e342..e47f7224e 100644
--- 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); }
diff --git a/poincare/include/poincare/subtraction.h b/poincare/include/poincare/subtraction.h
index 6d611b6bb..2bef24a6a 100644
--- 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>);
   }
 };
diff --git a/poincare/include/poincare/sum.h b/poincare/include/poincare/sum.h
index 1aae88773..80caab735 100644
--- 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;
 };
 
 }
diff --git a/poincare/include/poincare/symbol.h b/poincare/include/poincare/symbol.h
index edd1d1cd1..e25850b58 100644
--- 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;
 };
 
diff --git a/poincare/include/poincare/tangent.h b/poincare/include/poincare/tangent.h
index 640a2fad3..510782bf4 100644
--- 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>);
   }
 };
diff --git a/poincare/include/poincare/trigonometry.h b/poincare/include/poincare/trigonometry.h
index d09e6ec5f..a95cfa33f 100644
--- 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);
 };
 
 }
diff --git a/poincare/include/poincare/undefined.h b/poincare/include/poincare/undefined.h
index a0b1a0b24..3c6bfae08 100644
--- 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;
 };
 
diff --git a/poincare/include/poincare/variable_context.h b/poincare/include/poincare/variable_context.h
index bbd2d9ed2..eb9c6e93f 100644
--- 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;
 };
 
diff --git a/poincare/src/absolute_value.cpp b/poincare/src/absolute_value.cpp
index 5751149b8..0023f9317 100644
--- 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);
 }
 
 }
diff --git a/poincare/src/addition.cpp b/poincare/src/addition.cpp
index 5d2dd137a..3b12e9675 100644
--- 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>);
 
 }
diff --git a/poincare/src/approximation_engine.cpp b/poincare/src/approximation_engine.cpp
index da5902d51..92321ea6b 100644
--- 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>));
 
 
 }
diff --git a/poincare/src/arc_cosine.cpp b/poincare/src/arc_cosine.cpp
index 13fc28d2e..32536a951 100644
--- 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;
 }
 
 }
diff --git a/poincare/src/arc_sine.cpp b/poincare/src/arc_sine.cpp
index a593e92f1..f66d53629 100644
--- 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;
 }
 
 }
diff --git a/poincare/src/arc_tangent.cpp b/poincare/src/arc_tangent.cpp
index 7d5fd6d96..8def53bf1 100644
--- 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;
 }
 
 }
diff --git a/poincare/src/binomial_coefficient.cpp b/poincare/src/binomial_coefficient.cpp
index 51c0c62ba..05e2fe5fb 100644
--- 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));
 }
 
 
diff --git a/poincare/src/ceiling.cpp b/poincare/src/ceiling.cpp
index 4d8af70bb..ee097015a 100644
--- 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 {
diff --git a/poincare/src/complex.cpp b/poincare/src/complex.cpp
deleted file mode 100644
index 7ca360614..000000000
--- a/poincare/src/complex.cpp
+++ /dev/null
@@ -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;
-
-}
-
diff --git a/poincare/src/complex_argument.cpp b/poincare/src/complex_argument.cpp
index 99d0c4d4f..3c3eec2f7 100644
--- 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));
 }
 
 }
diff --git a/poincare/src/confidence_interval.cpp b/poincare/src/confidence_interval.cpp
index eed71f1e8..d97afff41 100644
--- 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);
 }
 
 }
diff --git a/poincare/src/conjugate.cpp b/poincare/src/conjugate.cpp
index 57461ae00..b98eee588 100644
--- 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);
 }
 
 }
diff --git a/poincare/src/cosine.cpp b/poincare/src/cosine.cpp
index b45c006e7..67d4cb099 100644
--- 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);
-}
-
 }
diff --git a/poincare/src/decimal.cpp b/poincare/src/decimal.cpp
index 96506dad3..05c4f7558 100644
--- 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 {
diff --git a/poincare/src/derivative.cpp b/poincare/src/derivative.cpp
index b7cacecdc..390db4f64 100644
--- 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>
diff --git a/poincare/src/determinant.cpp b/poincare/src/determinant.cpp
index 74cbeb04d..622565cb5 100644
--- 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;
 }
diff --git a/poincare/src/division.cpp b/poincare/src/division.cpp
index de9b576a8..64a1ea120 100644
--- 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;
 }
diff --git a/poincare/src/division_quotient.cpp b/poincare/src/division_quotient.cpp
index 23032d0d4..b671c2bf1 100644
--- 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));
 }
 
 }
diff --git a/poincare/src/division_remainder.cpp b/poincare/src/division_remainder.cpp
index bf6157b8a..1aec5b1a8 100644
--- a/poincare/src/division_remainder.cpp
+++ b/poincare/src/division_remainder.cpp
@@ -59,16 +59,16 @@ Expression * DivisionRemainder::shallowReduce(Context& context, AngleUnit angleU
 
 template<typename T>
 Complex<T> * DivisionRemainder::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::round(f1-f2*std::floor(f1/f2))));
+  return new Complex<T>(std::round(f1-f2*std::floor(f1/f2)));
 }
 
 }
diff --git a/poincare/src/empty_expression.cpp b/poincare/src/empty_expression.cpp
index 9e416269c..4b32ed986 100644
--- a/poincare/src/empty_expression.cpp
+++ b/poincare/src/empty_expression.cpp
@@ -22,7 +22,7 @@ ExpressionLayout * EmptyExpression::privateCreateLayout(PrintFloat::Mode floatDi
 }
 
 template<typename T> Complex<T> * EmptyExpression::templatedApproximate(Context& context, AngleUnit angleUnit) const {
-  return new Complex<T>(Complex<T>::Float(NAN));
+  return new Complex<T>(Complex<T>::Undefined());
 }
 
 }
diff --git a/poincare/src/equal.cpp b/poincare/src/equal.cpp
index 58db799c9..84c183450 100644
--- a/poincare/src/equal.cpp
+++ b/poincare/src/equal.cpp
@@ -6,7 +6,6 @@ extern "C" {
 }
 #include <poincare/equal.h>
 #include <ion.h>
-#include <poincare/complex.h>
 #include <poincare/rational.h>
 #include <poincare/addition.h>
 #include <poincare/division.h>
@@ -61,8 +60,8 @@ ExpressionLayout * Equal::privateCreateLayout(PrintFloat::Mode floatDisplayMode,
 }
 
 template<typename T>
-Expression * Equal::templatedApproximate(Context& context, AngleUnit angleUnit) const {
-  return new Complex<T>(Complex<T>::Float(NAN));
+Evaluation<T> * Equal::templatedApproximate(Context& context, AngleUnit angleUnit) const {
+  return new Complex<T>(Complex<T>::Undefined());
 }
 
 }
diff --git a/poincare/src/evaluation.cpp b/poincare/src/evaluation.cpp
new file mode 100644
index 000000000..73e5ee1d9
--- /dev/null
+++ b/poincare/src/evaluation.cpp
@@ -0,0 +1,117 @@
+extern "C" {
+#include <assert.h>
+#include <float.h>
+#include <stdlib.h>
+}
+#include <poincare/evaluation.h>
+#include <poincare/division.h>
+#include <poincare/matrix.h>
+#include <cmath>
+
+namespace Poincare {
+
+template<typename T>
+T Complex<T>::toScalar() const {
+  if (this->imag() == 0.0) {
+    return this->real();
+  }
+  return NAN;
+}
+
+template<typename T>
+Complex<T> * Complex<T>::createInverse() const {
+  return new Complex<T>(Division::compute(std::complex<T>(1.0), *this));
+}
+
+template<typename T>
+MatrixComplex<T>::MatrixComplex(std::complex<T> * operands, int numberOfRows, int numberOfColumns) :
+  m_numberOfRows(numberOfRows),
+  m_numberOfColumns(numberOfColumns)
+{
+  m_operands = new std::complex<T> [numberOfRows*numberOfColumns];
+  for (int i=0; i<numberOfRows*numberOfColumns; i++) {
+    m_operands[i] = operands[i];
+  }
+}
+
+template<typename T>
+std::complex<T> MatrixComplex<T>::createTrace() const {
+  if (numberOfRows() != numberOfColumns()) {
+    return std::complex<T>(NAN, NAN);
+  }
+  int dim = numberOfRows();
+  std::complex<T> c = std::complex<T>(0);
+  for (int i = 0; i < dim; i++) {
+    c += complexOperand(i*dim+i);
+  }
+  return c;
+}
+
+template<typename T>
+std::complex<T> MatrixComplex<T>::createDeterminant() const {
+  if (numberOfRows() != numberOfColumns()) {
+    return std::complex<T>(NAN, NAN);
+  }
+  std::complex<T> * operandsCopy = new std::complex<T> [m_numberOfRows*m_numberOfColumns];
+  for (int i=0; i<m_numberOfRows*m_numberOfColumns; i++) {
+    operandsCopy[i] = m_operands[i];
+  }
+  std::complex<T> determinant = std::complex<T>(1);
+  Matrix::ArrayRowCanonize(operandsCopy, m_numberOfRows, m_numberOfColumns, &determinant);
+  delete[] operandsCopy;
+  return determinant;
+}
+
+template<typename T>
+MatrixComplex<T> * MatrixComplex<T>::createInverse() const {
+  std::complex<T> * operandsCopy = new std::complex<T> [m_numberOfRows*m_numberOfColumns];
+  for (int i=0; i<m_numberOfRows*m_numberOfColumns; i++) {
+    operandsCopy[i] = m_operands[i];
+  }
+  int result = Matrix::ArrayInverse(operandsCopy, m_numberOfRows, m_numberOfColumns);
+  MatrixComplex<T> * inverse = nullptr;
+  if (result == 0) {
+    // Intentionally swapping dimensions for inverse, although it doesn't make a difference because it is square
+    inverse = new MatrixComplex<T>(operandsCopy, m_numberOfColumns, m_numberOfRows);
+  }
+  delete [] operandsCopy;
+  return inverse;
+}
+
+template<typename T>
+MatrixComplex<T> * MatrixComplex<T>::createTranspose() const {
+  std::complex<T> * operands = new std::complex<T> [numberOfComplexOperands()];
+  for (int i = 0; i < numberOfRows(); i++) {
+    for (int j = 0; j < numberOfColumns(); j++) {
+      operands[j*numberOfRows()+i] = complexOperand(i*numberOfColumns()+j);
+    }
+  }
+  // Intentionally swapping dimensions for transpose
+  MatrixComplex<T> * matrix = new MatrixComplex<T>(operands, numberOfColumns(), numberOfRows());
+  delete[] operands;
+  return matrix;
+}
+
+template<typename T>
+MatrixComplex<T> MatrixComplex<T>::createIdentity(int dim) {
+  std::complex<T> * operands = new std::complex<T> [dim*dim];
+  for (int i = 0; i < dim; i++) {
+    for (int j = 0; j < dim; j++) {
+      if (i == j) {
+        operands[i*dim+j] = std::complex<T>(1);
+      } else {
+        operands[i*dim+j] = std::complex<T>(0);
+      }
+    }
+  }
+  MatrixComplex<T> matrix = MatrixComplex(operands, dim, dim);
+  delete [] operands;
+  return matrix;
+}
+
+template class Complex<float>;
+template class Complex<double>;
+template class MatrixComplex<float>;
+template class MatrixComplex<double>;
+
+}
diff --git a/poincare/src/expression.cpp b/poincare/src/expression.cpp
index 22e30cfcd..34367e64a 100644
--- a/poincare/src/expression.cpp
+++ b/poincare/src/expression.cpp
@@ -9,9 +9,14 @@
 #include <poincare/simplification_root.h>
 #include <poincare/rational.h>
 #include <poincare/matrix.h>
-#include <poincare/complex.h>
 #include <poincare/opposite.h>
 #include <poincare/variable_context.h>
+#include <poincare/decimal.h>
+#include <poincare/ieee754.h>
+#include <poincare/addition.h>
+#include <poincare/multiplication.h>
+#include <poincare/power.h>
+#include <ion.h>
 #include <cmath>
 #include <float.h>
 #include "expression_parser.hpp"
@@ -178,7 +183,7 @@ bool Expression::recursivelyMatches(ExpressionTest test, Context & context) cons
 
 bool Expression::isApproximate(Context & context) const {
   return recursivelyMatches([](const Expression * e, Context & context) {
-        return e->type() == Expression::Type::Decimal || e->type() == Expression::Type::Complex || Expression::IsMatrix(e, context) || (e->type() == Expression::Type::Symbol && static_cast<const Symbol *>(e)->isApproximate(context));
+        return e->type() == Expression::Type::Decimal || Expression::IsMatrix(e, context) || (e->type() == Expression::Type::Symbol && static_cast<const Symbol *>(e)->isApproximate(context));
     }, context);
 }
 
@@ -439,21 +444,35 @@ Expression * Expression::deepBeautify(Context & context, AngleUnit angleUnit) {
 
 /* Evaluation */
 
-template<typename T> Expression * Expression::approximate(Context& context, AngleUnit angleUnit) const {
-  switch (angleUnit) {
-    case AngleUnit::Default:
-      return privateApproximate(T(), context, Preferences::sharedPreferences()->angleUnit());
-    default:
-      return privateApproximate(T(), context, angleUnit);
+template<typename T> Expression * Expression::approximate(Context& context, AngleUnit angleUnit, ComplexFormat complexFormat) const {
+  if (angleUnit == AngleUnit::Default) {
+    angleUnit = Preferences::sharedPreferences()->angleUnit();
   }
+  if (complexFormat == ComplexFormat::Default) {
+    complexFormat = Preferences::sharedPreferences()->complexFormat();
+  }
+  Expression * result = nullptr;
+  Evaluation<T> * e = privateApproximate(T(), context, angleUnit);
+  if (e->type() == Evaluation<T>::Type::Complex) {
+    result = complexToExpression(*(static_cast<Complex<T> *>(e)), complexFormat);
+  } else {
+    MatrixComplex<T> * matrix = static_cast<MatrixComplex<T> *>(e);
+    Expression ** operands = new Expression * [matrix->numberOfComplexOperands()];
+    for (int i = 0; i < matrix->numberOfComplexOperands(); i++) {
+      operands[i] = complexToExpression(matrix->complexOperand(i), complexFormat);
+    }
+    result = new Matrix(operands, matrix->numberOfRows(), matrix->numberOfColumns(), false);
+    delete[] operands;
+  }
+  delete e;
+  return result;
 }
 
-template<typename T> T Expression::approximateToScalar(Context& context, AngleUnit angleUnit) const {
-  Expression * evaluation = approximate<T>(context, angleUnit);
-  assert(evaluation->type() == Type::Complex || evaluation->type() == Type::Matrix);
+template<typename T> T Expression::approximateToScalar(Context& context, AngleUnit angleUnit, ComplexFormat complexFormat) const {
+  Expression * evaluation = approximate<T>(context, angleUnit, complexFormat);
   T result = NAN;
-  if (evaluation->type() == Type::Complex) {
-    result = static_cast<const Complex<T> *>(evaluation)->toScalar();
+  if (evaluation->type() == Type::Decimal) {
+    //result = static_cast<const Decimal *>(evaluation)->toScalar(); //TODO
   }
   /*if (evaluation->type() == Type::Matrix) {
     if (numberOfOperands() == 1) {
@@ -464,9 +483,9 @@ template<typename T> T Expression::approximateToScalar(Context& context, AngleUn
   return result;
 }
 
-template<typename T> T Expression::approximateToScalar(const char * text, Context& context, AngleUnit angleUnit) {
+template<typename T> T Expression::approximateToScalar(const char * text, Context& context, AngleUnit angleUnit, ComplexFormat complexFormat) {
   Expression * exp = ParseAndSimplify(text, context, angleUnit);
-  T result = exp->approximateToScalar<T>(context, angleUnit);
+  T result = exp->approximateToScalar<T>(context, angleUnit, complexFormat);
   delete exp;
   return result;
 }
@@ -799,21 +818,41 @@ double Expression::brentRoot(char symbol, double ax, double bx, double precision
 
 template<typename T>
 T Expression::approximateWithValueForSymbol(char symbol, T x, Context & context) const {
-  VariableContext<T> variableContext = VariableContext<T>(symbol, &context);
-  Symbol xSymbol(symbol);
-  Complex<T> e = Complex<T>::Float(x);
-  variableContext.setExpressionForSymbolName(&e, &xSymbol, variableContext);
-  return approximateToScalar<T>(variableContext);
+  VariableContext variableContext = VariableContext(symbol, &context);
+  variableContext.setApproximationForVariable(x);
+  T value = approximateToScalar<T>(variableContext);
+  return value;
+}
+
+template <typename T>
+Expression * Expression::CreateDecimal(T f) {
+  if (std::isnan(f) || std::isinf(f)) {
+    return new Undefined();
+  }
+  return new Decimal(f);
+}
+
+template<typename T> Expression * Expression::complexToExpression(std::complex<T> c, ComplexFormat complexFormat) {
+  if (complexFormat == ComplexFormat::Default) {
+    complexFormat = Preferences::sharedPreferences()->complexFormat();
+  }
+  switch (complexFormat) {
+    case ComplexFormat::Cartesian:
+      return new Addition(CreateDecimal(c.real()), new Multiplication(new Symbol(Ion::Charset::IComplex), CreateDecimal(c.imag()), false), false);
+    default:
+      assert(complexFormat == ComplexFormat::Polar);
+      return new Multiplication(CreateDecimal(std::abs(c)), new Power(new Symbol(Ion::Charset::Exponential), new Multiplication(new Symbol(Ion::Charset::IComplex), CreateDecimal(std::arg(c)), false), false), false);
+  }
 }
 
 }
 
-template Poincare::Expression * Poincare::Expression::approximate<double>(Context& context, AngleUnit angleUnit) const;
-template Poincare::Expression * Poincare::Expression::approximate<float>(Context& context, AngleUnit angleUnit) const;
-template double Poincare::Expression::approximateToScalar<double>(char const*, Poincare::Context&, Poincare::Expression::AngleUnit);
-template float Poincare::Expression::approximateToScalar<float>(char const*, Poincare::Context&, Poincare::Expression::AngleUnit);
-template double Poincare::Expression::approximateToScalar<double>(Poincare::Context&, Poincare::Expression::AngleUnit) const;
-template float Poincare::Expression::approximateToScalar<float>(Poincare::Context&, Poincare::Expression::AngleUnit) const;
+template Poincare::Expression * Poincare::Expression::approximate<double>(Context& context, AngleUnit angleUnit, ComplexFormat complexFormat) const;
+template Poincare::Expression * Poincare::Expression::approximate<float>(Context& context, AngleUnit angleUnit, ComplexFormat complexFormat) const;
+template double Poincare::Expression::approximateToScalar<double>(char const*, Poincare::Context&, Poincare::Expression::AngleUnit, ComplexFormat complexFormat);
+template float Poincare::Expression::approximateToScalar<float>(char const*, Poincare::Context&, Poincare::Expression::AngleUnit, ComplexFormat complexFormat);
+template double Poincare::Expression::approximateToScalar<double>(Poincare::Context&, Poincare::Expression::AngleUnit, ComplexFormat complexFormat) const;
+template float Poincare::Expression::approximateToScalar<float>(Poincare::Context&, Poincare::Expression::AngleUnit, ComplexFormat complexFormat) const;
 template double Poincare::Expression::epsilon<double>();
 template float Poincare::Expression::epsilon<float>();
 template double Poincare::Expression::approximateWithValueForSymbol(char, double, Poincare::Context&) const;
diff --git a/poincare/src/expression_debug.cpp b/poincare/src/expression_debug.cpp
index 0e96306af..ccc6c0160 100644
--- a/poincare/src/expression_debug.cpp
+++ b/poincare/src/expression_debug.cpp
@@ -3,7 +3,6 @@
 #include <poincare/symbol.h>
 #include <poincare/integer.h>
 #include <poincare/rational.h>
-#include <poincare/complex.h>
 #include <ion.h>
 #include <iostream>
 
@@ -39,13 +38,6 @@ void print_expression(const Expression * e, int indentationLevel) {
     case Expression::Type::Ceiling:
       std::cout << "Ceiling";
       break;
-    case Expression::Type::Complex:
-      std::cout << "Complex(";
-      std::cout << static_cast<const Complex<double> *>(e)->a();
-      std::cout << ", ";
-      std::cout << static_cast<const Complex<double> *>(e)->b();
-      std::cout << ")";
-      break;
     case Expression::Type::ComplexArgument:
       std::cout << "ComplexArgument";
       break;
diff --git a/poincare/src/factorial.cpp b/poincare/src/factorial.cpp
index bb1c048b7..a39bce388 100644
--- a/poincare/src/factorial.cpp
+++ b/poincare/src/factorial.cpp
@@ -77,19 +77,19 @@ Expression * Factorial::shallowBeautify(Context& context, AngleUnit angleUnit) {
 }
 
 template<typename T>
-Complex<T> Factorial::computeOnComplex(const Complex<T> c, AngleUnit angleUnit) {
-  T n = c.a();
-  if (c.b() != 0 || std::isnan(n) || n != (int)n || n < 0) {
-    return Complex<T>::Float(NAN);
+std::complex<T> Factorial::computeOnComplex(const std::complex<T> c, AngleUnit angleUnit) {
+  T n = c.real();
+  if (c.imag() != 0 || std::isnan(n) || n != (int)n || n < 0) {
+    return Complex<T>::Undefined();
   }
   T result = 1;
   for (int i = 1; i <= (int)n; i++) {
     result *= (T)i;
     if (std::isinf(result)) {
-      return Complex<T>::Float(result);
+      return Complex<T>(result);
     }
   }
-  return Complex<T>::Float(std::round(result));
+  return Complex<T>(std::round(result));
 }
 
 ExpressionLayout * Factorial::privateCreateLayout(PrintFloat::Mode floatDisplayMode, ComplexFormat complexFormat) const {
diff --git a/poincare/src/floor.cpp b/poincare/src/floor.cpp
index fa9b02daa..08e5300cb 100644
--- a/poincare/src/floor.cpp
+++ b/poincare/src/floor.cpp
@@ -49,11 +49,11 @@ Expression * Floor::shallowReduce(Context& context, AngleUnit angleUnit) {
 }
 
 template<typename T>
-Complex<T> Floor::computeOnComplex(const Complex<T> c, AngleUnit angleUnit) {
-  if (c.b() != 0) {
-    return Complex<T>::Float(NAN);
+std::complex<T> Floor::computeOnComplex(const std::complex<T> c, AngleUnit angleUnit) {
+  if (c.imag() != 0) {
+    return Complex<T>::Undefined();
   }
-  return Complex<T>::Float(std::floor(c.a()));
+  return Complex<T>(std::floor(c.real()));
 }
 
 ExpressionLayout * Floor::privateCreateLayout(PrintFloat::Mode floatDisplayMode, ComplexFormat complexFormat) const {
diff --git a/poincare/src/frac_part.cpp b/poincare/src/frac_part.cpp
index 0231b4e75..c6313e195 100644
--- a/poincare/src/frac_part.cpp
+++ b/poincare/src/frac_part.cpp
@@ -37,11 +37,11 @@ Expression * FracPart::shallowReduce(Context& context, AngleUnit angleUnit) {
 }
 
 template<typename T>
-Complex<T> FracPart::computeOnComplex(const Complex<T> c, AngleUnit angleUnit) {
-  if (c.b() != 0) {
-    return Complex<T>::Float(NAN);
+std::complex<T> FracPart::computeOnComplex(const std::complex<T> c, AngleUnit angleUnit) {
+  if (c.imag() != 0) {
+    return Complex<T>::Undefined();
   }
-  return Complex<T>::Float(c.a()-std::floor(c.a()));
+  return Complex<T>(c.real()-std::floor(c.real()));
 }
 
 }
diff --git a/poincare/src/global_context.cpp b/poincare/src/global_context.cpp
index 479e61cfe..93d4cc099 100644
--- a/poincare/src/global_context.cpp
+++ b/poincare/src/global_context.cpp
@@ -1,6 +1,6 @@
 #include <poincare/global_context.h>
 #include <poincare/matrix.h>
-#include <poincare/matrix.h>
+#include <poincare/undefined.h>
 #include <assert.h>
 #include <cmath>
 #include <ion.h>
@@ -8,9 +8,8 @@
 namespace Poincare {
 
 GlobalContext::GlobalContext() :
-  m_pi(Complex<double>::Float(M_PI)),
-  m_e(Complex<double>::Float(M_E)),
-  m_i(Complex<double>::Cartesian(0.0, 1.0))
+  m_pi(M_PI),
+  m_e(M_E)
 {
   for (int i = 0; i < k_maxNumberOfScalarExpressions; i++) {
     m_expressions[i] = nullptr;
@@ -40,8 +39,8 @@ GlobalContext::~GlobalContext() {
   }
 }
 
-Complex<double> * GlobalContext::defaultExpression() {
-  static Complex<double> defaultExpression(Complex<double>::Float(0.0));
+Decimal * GlobalContext::defaultExpression() {
+  static Decimal defaultExpression(0.0f);
   return &defaultExpression;
 }
 
@@ -62,9 +61,6 @@ const Expression * GlobalContext::expressionForSymbol(const Symbol * symbol) {
   if (symbol->name() == Ion::Charset::Exponential) {
     return &m_e;
   }
-  if (symbol->name() == Ion::Charset::IComplex) {
-    return &m_i;
-  }
   int index = symbolIndex(symbol);
   if (symbol->isMatrixSymbol()) {
     return m_matrixExpressions[index];
@@ -104,7 +100,7 @@ void GlobalContext::setExpressionForSymbolName(const Expression * expression, co
       m_matrixLayout[indexMatrix] = nullptr;
     }
     if (evaluation != nullptr) {
-      if (evaluation->type() == Expression::Type::Complex) {
+      if (evaluation->type() != Expression::Type::Matrix) {
         m_matrixExpressions[indexMatrix] = new Matrix(&evaluation, 1, 1, false);
       } else {
         m_matrixExpressions[indexMatrix] = static_cast<Matrix *>(evaluation);
@@ -123,11 +119,11 @@ void GlobalContext::setExpressionForSymbolName(const Expression * expression, co
   if (evaluation == nullptr) {
     return;
   }
-  if (evaluation->type() == Expression::Type::Complex) {
-    m_expressions[index] = static_cast<Complex<double> *>(evaluation);
-  } else {
-    m_expressions[index] = new Complex<double>(Complex<double>::Float(NAN));
+  if (evaluation->type() == Expression::Type::Matrix) {
+    m_expressions[index] = new Undefined();
     delete evaluation;
+  } else {
+    m_expressions[index] = evaluation;
   }
 }
 
diff --git a/poincare/src/great_common_divisor.cpp b/poincare/src/great_common_divisor.cpp
index 6cd10604c..99865c8b4 100644
--- a/poincare/src/great_common_divisor.cpp
+++ b/poincare/src/great_common_divisor.cpp
@@ -1,5 +1,4 @@
 #include <poincare/great_common_divisor.h>
-#include <poincare/complex.h>
 #include <poincare/undefined.h>
 #include <poincare/rational.h>
 #include <poincare/arithmetic.h>
@@ -58,14 +57,14 @@ Expression * GreatCommonDivisor::shallowReduce(Context& context, AngleUnit angle
 
 template<typename T>
 Complex<T> * GreatCommonDivisor::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());
   }
   int a = (int)f2;
   int b = (int)f1;
@@ -79,7 +78,7 @@ Complex<T> * GreatCommonDivisor::templatedApproximate(Context& context, AngleUni
     a = b;
     b = r;
   }
-  return new Complex<T>(Complex<T>::Float(std::round((T)a)));
+  return new Complex<T>(std::round((T)a));
 }
 
 }
diff --git a/poincare/src/hyperbolic_arc_cosine.cpp b/poincare/src/hyperbolic_arc_cosine.cpp
index 0d5fbb7c3..90af51886 100644
--- a/poincare/src/hyperbolic_arc_cosine.cpp
+++ b/poincare/src/hyperbolic_arc_cosine.cpp
@@ -31,11 +31,8 @@ Expression * HyperbolicArcCosine::shallowReduce(Context& context, AngleUnit angl
 }
 
 template<typename T>
-Complex<T> HyperbolicArcCosine::computeOnComplex(const Complex<T> c, AngleUnit angleUnit) {
-  if (c.b() != 0) {
-    return Complex<T>::Float(NAN);
-  }
-  return Complex<T>::Float(std::acosh(c.a()));
+std::complex<T> HyperbolicArcCosine::computeOnComplex(const std::complex<T> c, AngleUnit angleUnit) {
+  return std::acosh(c);
 }
 
 }
diff --git a/poincare/src/hyperbolic_arc_sine.cpp b/poincare/src/hyperbolic_arc_sine.cpp
index 6f68f1f21..476374da2 100644
--- a/poincare/src/hyperbolic_arc_sine.cpp
+++ b/poincare/src/hyperbolic_arc_sine.cpp
@@ -31,11 +31,8 @@ Expression * HyperbolicArcSine::shallowReduce(Context& context, AngleUnit angleU
 }
 
 template<typename T>
-Complex<T> HyperbolicArcSine::computeOnComplex(const Complex<T> c, AngleUnit angleUnit) {
-  if (c.b() != 0) {
-    return Complex<T>::Float(NAN);
-  }
-  return Complex<T>::Float(std::asinh(c.a()));
+std::complex<T> HyperbolicArcSine::computeOnComplex(const std::complex<T> c, AngleUnit angleUnit) {
+  return std::asinh(c);
 }
 
 }
diff --git a/poincare/src/hyperbolic_arc_tangent.cpp b/poincare/src/hyperbolic_arc_tangent.cpp
index 6717cedb0..b7e6bb74b 100644
--- a/poincare/src/hyperbolic_arc_tangent.cpp
+++ b/poincare/src/hyperbolic_arc_tangent.cpp
@@ -31,11 +31,8 @@ Expression * HyperbolicArcTangent::shallowReduce(Context& context, AngleUnit ang
 }
 
 template<typename T>
-Complex<T> HyperbolicArcTangent::computeOnComplex(const Complex<T> c, AngleUnit angleUnit) {
-  if (c.b() != 0) {
-    return Complex<T>::Float(NAN);
-  }
-  return Complex<T>::Float(std::atanh(c.a()));
+std::complex<T> HyperbolicArcTangent::computeOnComplex(const std::complex<T> c, AngleUnit angleUnit) {
+  return std::atanh(c);
 }
 
 }
diff --git a/poincare/src/hyperbolic_cosine.cpp b/poincare/src/hyperbolic_cosine.cpp
index a0e7d7c5e..dbccbd7c7 100644
--- a/poincare/src/hyperbolic_cosine.cpp
+++ b/poincare/src/hyperbolic_cosine.cpp
@@ -1,5 +1,4 @@
 #include <poincare/hyperbolic_cosine.h>
-#include <poincare/complex.h>
 #include <poincare/addition.h>
 #include <poincare/power.h>
 #include <poincare/division.h>
@@ -36,15 +35,8 @@ Expression * HyperbolicCosine::shallowReduce(Context& context, AngleUnit angleUn
 }
 
 template<typename T>
-Complex<T> HyperbolicCosine::computeOnComplex(const Complex<T> c, AngleUnit angleUnit) {
-  if (c.b() == 0) {
-    return Complex<T>::Float(std::cosh(c.a()));
-  }
-  Complex<T> e = Complex<T>::Float(M_E);
-  Complex<T> exp1 = Power::compute(e, c);
-  Complex<T> exp2 = Power::compute(e, Complex<T>::Cartesian(-c.a(), -c.b()));
-  Complex<T> sum = Addition::compute(exp1, exp2);
-  return Division::compute(sum, Complex<T>::Float(2));
+std::complex<T> HyperbolicCosine::computeOnComplex(const std::complex<T> c, AngleUnit angleUnit) {
+  return std::cosh(c);
 }
 
 }
diff --git a/poincare/src/hyperbolic_sine.cpp b/poincare/src/hyperbolic_sine.cpp
index 1b5ce48c9..588b084eb 100644
--- a/poincare/src/hyperbolic_sine.cpp
+++ b/poincare/src/hyperbolic_sine.cpp
@@ -1,5 +1,4 @@
 #include <poincare/hyperbolic_sine.h>
-#include <poincare/complex.h>
 #include <poincare/subtraction.h>
 #include <poincare/power.h>
 #include <poincare/division.h>
@@ -36,15 +35,8 @@ Expression * HyperbolicSine::shallowReduce(Context& context, AngleUnit angleUnit
 }
 
 template<typename T>
-Complex<T> HyperbolicSine::computeOnComplex(const Complex<T> c, AngleUnit angleUnit) {
-  if (c.b() == 0) {
-    return Complex<T>::Float(std::sinh(c.a()));
-  }
-  Complex<T> e = Complex<T>::Float(M_E);
-  Complex<T> exp1 = Power::compute(e, c);
-  Complex<T> exp2 = Power::compute(e, Complex<T>::Cartesian(-c.a(), -c.b()));
-  Complex<T> sub = Subtraction::compute(exp1, exp2);
-  return Division::compute(sub, Complex<T>::Float(2));
+std::complex<T> HyperbolicSine::computeOnComplex(const std::complex<T> c, AngleUnit angleUnit) {
+  return std::sinh(c);
 }
 
 }
diff --git a/poincare/src/hyperbolic_tangent.cpp b/poincare/src/hyperbolic_tangent.cpp
index 96ccbb6ae..1df271e57 100644
--- a/poincare/src/hyperbolic_tangent.cpp
+++ b/poincare/src/hyperbolic_tangent.cpp
@@ -1,7 +1,6 @@
 #include <poincare/hyperbolic_tangent.h>
 #include <poincare/hyperbolic_cosine.h>
 #include <poincare/hyperbolic_sine.h>
-#include <poincare/complex.h>
 #include <poincare/division.h>
 #include <poincare/simplification_engine.h>
 extern "C" {
@@ -35,13 +34,8 @@ Expression * HyperbolicTangent::shallowReduce(Context& context, AngleUnit angleU
 }
 
 template<typename T>
-Complex<T> HyperbolicTangent::computeOnComplex(const Complex<T> c, AngleUnit angleUnit) {
-  if (c.b() == 0) {
-    return Complex<T>::Float(std::tanh(c.a()));
-  }
-  Complex<T> arg1 = HyperbolicSine::computeOnComplex(c, angleUnit);
-  Complex<T> arg2 = HyperbolicCosine::computeOnComplex(c, angleUnit);
-  return Division::compute(arg1, arg2);
+std::complex<T> HyperbolicTangent::computeOnComplex(const std::complex<T> c, AngleUnit angleUnit) {
+  return std::tanh(c);
 }
 
 }
diff --git a/poincare/src/imaginary_part.cpp b/poincare/src/imaginary_part.cpp
index 506141ef3..216fd20ac 100644
--- a/poincare/src/imaginary_part.cpp
+++ b/poincare/src/imaginary_part.cpp
@@ -1,5 +1,4 @@
 #include <poincare/imaginary_part.h>
-#include <poincare/complex.h>
 #include <poincare/simplification_engine.h>
 #include <poincare/rational.h>
 #include <cmath>
@@ -37,8 +36,8 @@ Expression * ImaginaryPart::shallowReduce(Context& context, AngleUnit angleUnit)
 }
 
 template<typename T>
-Complex<T> ImaginaryPart::computeOnComplex(const Complex<T> c, AngleUnit angleUnit) {
-  return Complex<T>::Float(c.b());
+std::complex<T> ImaginaryPart::computeOnComplex(const std::complex<T> c, AngleUnit angleUnit) {
+  return std::imag(c);
 }
 
 }
diff --git a/poincare/src/integer.cpp b/poincare/src/integer.cpp
index a55a19dac..3e959cda2 100644
--- a/poincare/src/integer.cpp
+++ b/poincare/src/integer.cpp
@@ -1,5 +1,4 @@
 #include <poincare/integer.h>
-#include <poincare/complex.h>
 #include <poincare/ieee754.h>
 #include <poincare/layout_engine.h>
 #include <cmath>
@@ -9,6 +8,9 @@ extern "C" {
 #include <string.h>
 #include <assert.h>
 }
+#include <cmath>
+#include <poincare/ieee754.h>
+#include <utility>
 
 namespace Poincare {
 
diff --git a/poincare/src/integral.cpp b/poincare/src/integral.cpp
index 3246ee5fd..f9269dcde 100644
--- a/poincare/src/integral.cpp
+++ b/poincare/src/integral.cpp
@@ -49,21 +49,21 @@ Expression * Integral::shallowReduce(Context& context, AngleUnit angleUnit) {
 
 template<typename T>
 Complex<T> * Integral::templatedApproximate(Context & context, AngleUnit angleUnit) const {
-  Expression * aInput = operand(1)->approximate<T>(context, angleUnit);
-  T a = aInput->type() == Type::Complex ? static_cast<Complex<T> *>(aInput)->toScalar() : NAN;
+  Evaluation<T> * aInput = operand(1)->privateApproximate(T(), context, angleUnit);
+  T a = aInput->toScalar();
   delete aInput;
-  Expression * bInput = operand(2)->approximate<T>(context, angleUnit);
-  T b = bInput->type() == Type::Complex ? static_cast<Complex<T> *>(bInput)->toScalar() : NAN;
+  Evaluation<T> * bInput = operand(2)->privateApproximate(T(), context, angleUnit);
+  T b = bInput->toScalar();
   delete bInput;
   if (std::isnan(a) || std::isnan(b)) {
-    return new Complex<T>(Complex<T>::Float(NAN));
+    return new Complex<T>(Complex<T>::Undefined());
   }
 #ifdef LAGRANGE_METHOD
   T result = lagrangeGaussQuadrature<T>(a, b, context, angleUnit);
 #else
   T result = adaptiveQuadrature<T>(a, b, 0.1, k_maxNumberOfIterations, context, angleUnit);
 #endif
-  return new Complex<T>(Complex<T>::Float(result));
+  return new Complex<T>(result);
 }
 
 ExpressionLayout * Integral::privateCreateLayout(PrintFloat::Mode floatDisplayMode, ComplexFormat complexFormat) const {
diff --git a/poincare/src/least_common_multiple.cpp b/poincare/src/least_common_multiple.cpp
index 43bb02e92..fe2d82702 100644
--- a/poincare/src/least_common_multiple.cpp
+++ b/poincare/src/least_common_multiple.cpp
@@ -1,5 +1,4 @@
 #include <poincare/least_common_multiple.h>
-#include <poincare/complex.h>
 #include <poincare/rational.h>
 #include <poincare/undefined.h>
 #include <poincare/arithmetic.h>
@@ -58,17 +57,17 @@ Expression * LeastCommonMultiple::shallowReduce(Context& context, AngleUnit angl
 
 template<typename T>
 Complex<T> * LeastCommonMultiple::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());
   }
   if (f1 == 0.0f || f2 == 0.0f) {
-    return new Complex<T>(Complex<T>::Float(0));
+    return new Complex<T>(0);
   }
   int a = (int)f2;
   int b = (int)f1;
@@ -83,7 +82,7 @@ Complex<T> * LeastCommonMultiple::templatedApproximate(Context& context, AngleUn
     a = b;
     b = r;
   }
-  return new Complex<T>(Complex<T>::Float(product/a));
+  return new Complex<T>(product/a);
 }
 
 }
diff --git a/poincare/src/logarithm.cpp b/poincare/src/logarithm.cpp
index a8e979dc0..40fb8f9c2 100644
--- a/poincare/src/logarithm.cpp
+++ b/poincare/src/logarithm.cpp
@@ -30,11 +30,8 @@ Expression * Logarithm::clone() const {
 }
 
 template<typename T>
-Complex<T> Logarithm::computeOnComplex(const Complex<T> c, AngleUnit angleUnit) {
-  if (c.b() != 0) {
-    return Complex<T>::Float(NAN);
-  }
-  return Complex<T>::Float(std::log10(c.a()));
+std::complex<T> Logarithm::computeOnComplex(const std::complex<T> c, AngleUnit angleUnit) {
+  return std::log10(c);
 }
 
 Expression * Logarithm::simpleShallowReduce(Context & context, AngleUnit angleUnit) {
@@ -211,14 +208,14 @@ Expression * Logarithm::shallowBeautify(Context & context, AngleUnit angleUnit)
 }
 
 template<typename T>
-Expression * Logarithm::templatedApproximate(Context& context, AngleUnit angleUnit) const {
+Evaluation<T> * Logarithm::templatedApproximate(Context& context, AngleUnit angleUnit) const {
   if (numberOfOperands() == 1) {
     return ApproximationEngine::map(this, context, angleUnit, computeOnComplex<T>);
   }
-  Expression * x = operand(0)->approximate<T>(context, angleUnit);
-  Expression * n = operand(1)->approximate<T>(context, angleUnit);
-  Complex<T> result = Complex<T>::Float(NAN);
-  if (x->type() == Type::Complex && n->type() == Type::Complex) {
+  Evaluation<T> * x = operand(0)->privateApproximate(T(), context, angleUnit);
+  Evaluation<T> * n = operand(1)->privateApproximate(T(), context, angleUnit);
+  std::complex<T> result = std::complex<T>(NAN, NAN);
+  if (x->type() == Evaluation<T>::Type::Complex && n->type() == Evaluation<T>::Type::Complex) {
     Complex<T> * xc = static_cast<Complex<T> *>(x);
     Complex<T> * nc = static_cast<Complex<T> *>(n);
     result = Division::compute<T>(computeOnComplex(*xc, angleUnit), computeOnComplex(*nc, angleUnit));
diff --git a/poincare/src/matrix.cpp b/poincare/src/matrix.cpp
index d559d0c83..e8ad9e4f8 100644
--- a/poincare/src/matrix.cpp
+++ b/poincare/src/matrix.cpp
@@ -4,7 +4,6 @@ extern "C" {
 }
 #include <poincare/global_context.h>
 #include <poincare/matrix.h>
-#include <poincare/complex.h>
 #include <poincare/addition.h>
 #include <poincare/decimal.h>
 #include <poincare/undefined.h>
@@ -135,7 +134,7 @@ void Matrix::rowCanonize(Context & context, AngleUnit angleUnit, Multiplication
       }
       /* Set to 1 M[h][k] by linear combination */
       Expression * divisor = matrixOperand(h, k);
-      // Update determinant: det *= 1/divisor
+      // Update determinant: det *= divisor
       if (determinant) { determinant->addOperand(divisor->clone()); }
       for (int j = k+1; j < n; j++) {
         Expression * opHJ = matrixOperand(h, j);
@@ -165,19 +164,21 @@ void Matrix::rowCanonize(Context & context, AngleUnit angleUnit, Multiplication
 }
 
 template<typename T>
-void Matrix::ArrayRowCanonize(T * array, int numberOfRows, int numberOfColumns) {
+void Matrix::ArrayRowCanonize(T * array, int numberOfRows, int numberOfColumns, T * determinant) {
   int h = 0; // row pivot
   int k = 0; // column pivot
 
   while (h < numberOfRows && k < numberOfColumns) {
     // Find the first non-null pivot
     int iPivot = h;
-    while (iPivot < numberOfRows && std::fabs(array[iPivot*numberOfColumns+k]) < Expression::epsilon<T>()) {
+    while (iPivot < numberOfRows && std::abs(array[iPivot*numberOfColumns+k]) < Expression::epsilon<double>()) {
       iPivot++;
     }
     if (iPivot == numberOfRows) {
       // No non-null coefficient in this column, skip
       k++;
+      // Update determinant: det *= 0
+      if (determinant) { *determinant *= 0.0; }
     } else {
       // Swap row h and iPivot
       if (iPivot != h) {
@@ -187,9 +188,13 @@ void Matrix::ArrayRowCanonize(T * array, int numberOfRows, int numberOfColumns)
           array[iPivot*numberOfColumns+col] = array[h*numberOfColumns+col];
           array[h*numberOfColumns+col] = temp;
         }
+        // Update determinant: det *= -1
+        if (determinant) { *determinant *= -1.0; }
       }
       /* Set to 1 array[h][k] by linear combination */
       T divisor = array[h*numberOfColumns+k];
+      // Update determinant: det *= divisor
+      if (determinant) { *determinant *= divisor; }
       for (int j = k+1; j < numberOfColumns; j++) {
         array[h*numberOfColumns+j] /= divisor;
       }
@@ -246,36 +251,66 @@ int Matrix::rank(Context & context, AngleUnit angleUnit, bool inPlace) {
 }
 
 template<typename T>
-Complex<T> * Matrix::createTrace() const {
-  if (numberOfRows() != numberOfColumns()) {
-    return new Complex<T>(Complex<T>::Float(NAN));
+int Matrix::ArrayInverse(T * array, int numberOfRows, int numberOfColumns) {
+  if (numberOfRows != numberOfColumns) {
+    return -1;
   }
-  int dim = numberOfRows();
-  Complex<T> c = Complex<T>::Float(0);
+  int dim = numberOfRows;
+  /* Create the matrix inv = (A|I) with A the input matrix and I the dim identity matrix */
+  T * operands = new T[dim*2*dim];
   for (int i = 0; i < dim; i++) {
-    assert(operand(i*dim+i)->type() == Type::Complex);
-    c = Addition::compute(c, *(static_cast<const Complex<T> *>(operand(i*dim+i))));
+    for (int j = 0; j < dim; j++) {
+      operands[i*2*dim+j] = array[i*numberOfColumns+j];
+    }
+    for (int j = dim; j < 2*dim; j++) {
+      operands[i*2*dim+j] = j-dim == i ? 1 : 0;
+    }
   }
-  return new Complex<T>(c);
+  ArrayRowCanonize(operands, dim, 2*dim);
+  // Check inversibility
+  for (int i = 0; i < dim; i++) {
+    T one = 1.0;
+    if (std::abs(operands[i*2*dim+i] - one) > Expression::epsilon<float>()) {
+      return -2;
+    }
+  }
+  for (int i = 0; i < dim; i++) {
+    for (int j = 0; j < dim; j++) {
+      array[i*numberOfColumns+j] = operands[i*2*dim+j+dim];
+    }
+  }
+  delete [] operands;
+  return 0;
 }
 
-// TODO: 2. implement determinant/inverse for any expression (do not evaluate first)
-template<typename T>
-Complex<T> * Matrix::createDeterminant(Context & context, AngleUnit angleUnit) const {
-  if (numberOfRows() != numberOfColumns()) {
-    return new Complex<T>(Complex<T>::Float(NAN));
+#if MATRIX_EXACT_REDUCING
+Matrix * Matrix::createTranspose() const {
+  const Expression ** operands = new const Expression * [numberOfOperands()];
+  for (int i = 0; i < numberOfRows(); i++) {
+    for (int j = 0; j < numberOfColumns(); j++) {
+      operands[j*numberOfRows()+i] = operand(i*numberOfColumns()+j);
+    }
   }
-  Matrix * copy = static_cast<Matrix *>(clone());
-  Multiplication * m = new Multiplication();
-  copy->rowCanonize(context, angleUnit, m);
-  delete copy;
-  Expression * determinant = m->approximate<T>(context, angleUnit);
-  delete m;
-  if (determinant->type() != Type::Complex) {
-    delete determinant;
-    return new Complex<T>(Complex<T>::Float(NAN));
+  // Intentionally swapping dimensions for transpose
+  Matrix * matrix = new Matrix(operands, numberOfColumns(), numberOfRows(), true);
+  delete[] operands;
+  return matrix;
+}
+
+Matrix * Matrix::createIdentity(int dim) {
+  Expression ** operands = new Expression * [dim*dim];
+  for (int i = 0; i < dim; i++) {
+    for (int j = 0; j < dim; j++) {
+      if (i == j) {
+        operands[i*dim+j] = new Rational(1);
+      } else {
+        operands[i*dim+j] = new Rational(0);
+      }
+    }
   }
-  return static_cast<Complex<T> *>(determinant);
+  Matrix * matrix = new Matrix(operands, dim, dim, false);
+  delete [] operands;
+  return matrix;
 }
 
 Expression * Matrix::createInverse(Context & context, AngleUnit angleUnit) const {
@@ -316,130 +351,31 @@ Expression * Matrix::createInverse(Context & context, AngleUnit angleUnit) const
   return inverse;
 }
 
+#endif
+
 template<typename T>
-Matrix * Matrix::createApproximateInverse() const {
+Evaluation<T> * Matrix::templatedApproximate(Context& context, AngleUnit angleUnit) const {
+  std::complex<T> * operands = new std::complex<T> [numberOfOperands()];
   for (int i = 0; i < numberOfOperands(); i++) {
-    assert(operand(i)->type() == Type::Complex);
-  }
-  /* This method is called on already approximate matrices. The context and
-   * angle unit parameters are therefore useless in createInverse and
-   * approximate (as we work on Complex expressions only - no symbols or
-   * trigonometric functions would required a context and angle unit). We feed
-   * them with default parameters.*/
-  GlobalContext context;
-  Expression * inv = createInverse(context, AngleUnit::Degree);
-  Expression * invApproximation = inv->approximate<T>(context, AngleUnit::Degree);
-  delete inv;
-  if (invApproximation->type() == Type::Matrix) {
-    return static_cast<Matrix *>(invApproximation);
-  }
-  delete invApproximation;
-  return nullptr;
-}
-
-template<typename T>
-int Matrix::ArrayInverse(T * array, int numberOfRows, int numberOfColumns) {
-  if (numberOfRows != numberOfColumns) {
-    return -1;
-  }
-  int dim = numberOfRows;
-  /* Create the matrix inv = (A|I) with A the input matrix and I the dim identity matrix */
-  T * operands = new T[dim*2*dim];
-  for (int i = 0; i < dim; i++) {
-    for (int j = 0; j < dim; j++) {
-      operands[i*2*dim+j] = array[i*numberOfColumns+j];
-    }
-    for (int j = dim; j < 2*dim; j++) {
-      operands[i*2*dim+j] = j-dim == i ? 1 : 0;
-    }
-  }
-  ArrayRowCanonize(operands, dim, 2*dim);
-  // Check inversibility
-  for (int i = 0; i < dim; i++) {
-    if (operands[i*2*dim+i] != 1) {
-      return -2;
-    }
-  }
-  for (int i = 0; i < dim; i++) {
-    for (int j = 0; j < dim; j++) {
-      array[i*numberOfColumns+j] = operands[i*2*dim+j+dim];
-    }
-  }
-  delete [] operands;
-  return 0;
-}
-
-Matrix * Matrix::createTranspose() const {
-  const Expression ** operands = new const Expression * [numberOfOperands()];
-  for (int i = 0; i < numberOfRows(); i++) {
-    for (int j = 0; j < numberOfColumns(); j++) {
-      operands[j*numberOfRows()+i] = operand(i*numberOfColumns()+j);
-    }
-  }
-  // Intentionally swapping dimensions for transpose
-  Matrix * matrix = new Matrix(operands, numberOfColumns(), numberOfRows(), true);
-  delete[] operands;
-  return matrix;
-}
-
-Matrix * Matrix::createIdentity(int dim) {
-  Expression ** operands = new Expression * [dim*dim];
-  for (int i = 0; i < dim; i++) {
-    for (int j = 0; j < dim; j++) {
-      if (i == j) {
-        operands[i*dim+j] = new Rational(1);
-      } else {
-        operands[i*dim+j] = new Rational(0);
-      }
-    }
-  }
-  Matrix * matrix = new Matrix(operands, dim, dim, false);
-  delete [] operands;
-  return matrix;
-}
-
-template<typename T>
-Matrix * Matrix::createApproximateIdentity(int dim) {
-  Expression ** operands = new Expression * [dim*dim];
-  for (int i = 0; i < dim; i++) {
-    for (int j = 0; j < dim; j++) {
-      if (i == j) {
-        operands[i*dim+j] = new Complex<T>(Complex<T>::Float(1));
-      } else {
-        operands[i*dim+j] = new Complex<T>(Complex<T>::Float(0));
-      }
-    }
-  }
-  Matrix * matrix = new Matrix(operands, dim, dim, false);
-  delete [] operands;
-  return matrix;
-}
-
-template<typename T>
-Expression * Matrix::templatedApproximate(Context& context, AngleUnit angleUnit) const {
-  Expression ** operands = new Expression * [numberOfOperands()];
-  for (int i = 0; i < numberOfOperands(); i++) {
-    Expression * operandEvaluation = operand(i)->approximate<T>(context, angleUnit);
-    if (operandEvaluation->type() != Type::Complex) {
-      operands[i] = new Complex<T>(Complex<T>::Float(NAN));
-      delete operandEvaluation;
+    Evaluation<T> * operandEvaluation = operand(i)->privateApproximate(T(), context, angleUnit);
+    if (operandEvaluation->type() != Evaluation<T>::Type::Complex) {
+      operands[i] = Complex<T>::Undefined();
     } else {
-      operands[i] = operandEvaluation;
+      std::complex<T> * c = static_cast<Complex<T> *>(operandEvaluation);
+      operands[i] = *c;
     }
+    delete operandEvaluation;
   }
-  Expression * matrix = new Matrix(operands, numberOfRows(), numberOfColumns(), false);
+  MatrixComplex<T> * matrix = new MatrixComplex<T>(operands, numberOfRows(), numberOfColumns());
   delete[] operands;
   return matrix;
 }
 
-template Complex<float>* Matrix::createTrace<float>() const;
-template Complex<double>* Matrix::createTrace<double>() const;
-template Complex<float>* Matrix::createDeterminant<float>(Context &, AngleUnit) const;
-template Complex<double>* Matrix::createDeterminant<double>(Context &, AngleUnit) const;
-template Matrix* Matrix::createApproximateInverse<float>() const;
-template Matrix* Matrix::createApproximateInverse<double>() const;
-template Matrix* Matrix::createApproximateIdentity<float>(int);
-template Matrix* Matrix::createApproximateIdentity<double>(int);
 template int Matrix::ArrayInverse<float>(float *, int, int);
 template int Matrix::ArrayInverse<double>(double *, int, int);
+template int Matrix::ArrayInverse<std::complex<float>>(std::complex<float> *, int, int);
+template int Matrix::ArrayInverse<std::complex<double>>(std::complex<double> *, int, int);
+template void Matrix::ArrayRowCanonize<std::complex<float> >(std::complex<float>*, int, int, std::complex<float>*);
+template void Matrix::ArrayRowCanonize<std::complex<double> >(std::complex<double>*, int, int, std::complex<double>*);
+
 }
diff --git a/poincare/src/matrix_data.cpp b/poincare/src/matrix_data.cpp
index 4ef93f708..af46cf356 100644
--- a/poincare/src/matrix_data.cpp
+++ b/poincare/src/matrix_data.cpp
@@ -3,7 +3,6 @@ extern "C" {
 #include <stdlib.h>
 }
 #include <poincare/matrix_data.h>
-#include <poincare/complex.h>
 
 namespace Poincare {
 
diff --git a/poincare/src/matrix_dimension.cpp b/poincare/src/matrix_dimension.cpp
index f582f93cc..f431f4b9a 100644
--- a/poincare/src/matrix_dimension.cpp
+++ b/poincare/src/matrix_dimension.cpp
@@ -40,18 +40,18 @@ Expression * MatrixDimension::shallowReduce(Context& context, AngleUnit angleUni
 }
 
 template<typename T>
-Expression * MatrixDimension::templatedApproximate(Context& context, AngleUnit angleUnit) const {
-  Expression * input = operand(0)->approximate<T>(context, angleUnit);
-  Expression * operands[2];
-  if (input->type() == Type::Matrix) {
-    operands[0] = new Complex<T>(Complex<T>::Float((T)static_cast<Matrix *>(input)->numberOfRows()));
-    operands[1] = new Complex<T>(Complex<T>::Float((T)static_cast<Matrix *>(input)->numberOfColumns()));
+Evaluation<T> * MatrixDimension::templatedApproximate(Context& context, AngleUnit angleUnit) const {
+  Evaluation<T> * input = operand(0)->privateApproximate(T(), context, angleUnit);
+  std::complex<T> operands[2];
+  if (input->type() == Evaluation<T>::Type::MatrixComplex) {
+    operands[0] = std::complex<T>(static_cast<MatrixComplex<T> *>(input)->numberOfRows());
+    operands[1] = std::complex<T>(static_cast<MatrixComplex<T> *>(input)->numberOfColumns());
   } else {
-    operands[0] = new Complex<T>(Complex<T>::Float(1.0));
-    operands[1] = new Complex<T>(Complex<T>::Float(1.0));
+    operands[0] = std::complex<T>(1.0);
+    operands[1] = std::complex<T>(1.0);
   }
   delete input;
-  return new Matrix(operands, 1, 2, false);
+  return new MatrixComplex<T>(operands, 1, 2);
 }
 
 
diff --git a/poincare/src/matrix_inverse.cpp b/poincare/src/matrix_inverse.cpp
index 1f60cfbce..3dc4468f3 100644
--- a/poincare/src/matrix_inverse.cpp
+++ b/poincare/src/matrix_inverse.cpp
@@ -1,6 +1,5 @@
 #include <poincare/matrix_inverse.h>
 #include <poincare/matrix.h>
-#include <poincare/complex.h>
 #include <poincare/division.h>
 #include <poincare/undefined.h>
 #include <poincare/power.h>
@@ -46,21 +45,14 @@ Expression * MatrixInverse::shallowReduce(Context& context, AngleUnit angleUnit)
 
 // TODO: handle this exactly in shallowReduce for small dimensions.
 template<typename T>
-Expression * MatrixInverse::templatedApproximate(Context& context, AngleUnit angleUnit) const {
-  Expression * input = operand(0)->approximate<T>(context, angleUnit);
-  Expression * result = nullptr;
-  if (input->type() == Type::Complex) {
-    Complex<T> * c = static_cast<Complex<T> *>(input);
-    result = new Complex<T>(Division::compute(Complex<T>::Cartesian(1, 0), *c));
-  } else {
-    assert(input->type() == Type::Matrix);
-    result = static_cast<Matrix *>(input)->createApproximateInverse<T>();
-  }
-  if (result == nullptr) {
-    result = new Complex<T>(Complex<T>::Float(NAN));
+Evaluation<T> * MatrixInverse::templatedApproximate(Context& context, AngleUnit angleUnit) const {
+  Evaluation<T> * input = operand(0)->privateApproximate(T(), context, angleUnit);
+  Evaluation<T> * inverse = input->createInverse();
+  if (inverse == nullptr) {
+    inverse = new Complex<T>(Complex<T>::Undefined());
   }
   delete input;
-  return result;
+  return inverse;
 }
 
 }
diff --git a/poincare/src/matrix_trace.cpp b/poincare/src/matrix_trace.cpp
index 0525b1590..8af363820 100644
--- a/poincare/src/matrix_trace.cpp
+++ b/poincare/src/matrix_trace.cpp
@@ -49,15 +49,9 @@ Expression * MatrixTrace::shallowReduce(Context& context, AngleUnit angleUnit) {
 }
 
 template<typename T>
-Expression * MatrixTrace::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)->createTrace<T>();
-  }
+Complex<T> * MatrixTrace::templatedApproximate(Context& context, AngleUnit angleUnit) const {
+  Evaluation<T> * input = operand(0)->privateApproximate(T(), context, angleUnit);
+  Complex<T> * result = new Complex<T>(input->createTrace());
   delete input;
   return result;
 }
diff --git a/poincare/src/matrix_transpose.cpp b/poincare/src/matrix_transpose.cpp
index c6deb029d..3e548333c 100644
--- a/poincare/src/matrix_transpose.cpp
+++ b/poincare/src/matrix_transpose.cpp
@@ -1,6 +1,5 @@
 #include <poincare/matrix_transpose.h>
 #include <poincare/matrix.h>
-#include <poincare/complex.h>
 #include <poincare/division.h>
 extern "C" {
 #include <assert.h>
@@ -39,17 +38,12 @@ Expression * MatrixTranspose::shallowReduce(Context& context, AngleUnit angleUni
 }
 
 template<typename T>
-Expression * MatrixTranspose::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)->createTranspose();
-  }
+Evaluation<T> * MatrixTranspose::templatedApproximate(Context& context, AngleUnit angleUnit) const {
+  Evaluation<T> * input = operand(0)->privateApproximate(T(), context, angleUnit);
+  Evaluation<T> * transpose = input->createTranspose();
+  assert(transpose != nullptr);
   delete input;
-  return result;
+  return transpose;
 }
 
 }
diff --git a/poincare/src/multiplication.cpp b/poincare/src/multiplication.cpp
index b0cf74ffd..ec3edfb49 100644
--- a/poincare/src/multiplication.cpp
+++ b/poincare/src/multiplication.cpp
@@ -114,30 +114,26 @@ Expression * Multiplication::setSign(Sign s, Context & context, AngleUnit angleU
 }
 
 template<typename T>
-Complex<T> Multiplication::compute(const Complex<T> c, const Complex<T> d) {
-  return Complex<T>::Cartesian(c.a()*d.a()-c.b()*d.b(), c.b()*d.a() + c.a()*d.b());
+std::complex<T> Multiplication::compute(const std::complex<T> c, const std::complex<T> d) {
+  return c*d;
 }
 
 template<typename T>
-Matrix * Multiplication::computeOnMatrices(const Matrix * m, const Matrix * n) {
-  if (m->numberOfColumns() != n->numberOfRows()) {
-    return nullptr;
+MatrixComplex<T> Multiplication::computeOnMatrices(const MatrixComplex<T> m, const MatrixComplex<T> n) {
+  if (m.numberOfColumns() != n.numberOfRows()) {
+    return MatrixComplex<T>::Undefined();
   }
-  Expression ** operands = new Expression * [m->numberOfRows()*n->numberOfColumns()];
-  for (int i = 0; i < m->numberOfRows(); i++) {
-    for (int j = 0; j < n->numberOfColumns(); j++) {
-      T a = 0.0f;
-      T b = 0.0f;
-      for (int k = 0; k < m->numberOfColumns(); k++) {
-        const Complex<T> * mEntry = static_cast<const Complex<T> *>(m->operand(i*m->numberOfColumns()+k));
-        const Complex<T> * nEntry = static_cast<const Complex<T> *>(n->operand(k*n->numberOfColumns()+j));
-        a += mEntry->a()*nEntry->a() - mEntry->b()*nEntry->b();
-        b += mEntry->b()*nEntry->a() + mEntry->a()*nEntry->b();
+  std::complex<T> * operands = new std::complex<T> [m.numberOfRows()*n.numberOfColumns()];
+  for (int i = 0; i < m.numberOfRows(); i++) {
+    for (int j = 0; j < n.numberOfColumns(); j++) {
+      std::complex<T> c(0.0);
+      for (int k = 0; k < m.numberOfColumns(); k++) {
+        c += m.complexOperand(i*m.numberOfColumns()+k)*n.complexOperand(k*n.numberOfColumns()+j);
       }
-      operands[i*n->numberOfColumns()+j] = new Complex<T>(Complex<T>::Cartesian(a, b));
+      operands[i*n.numberOfColumns()+j] = c;
     }
   }
-  Matrix * result = new Matrix(operands, m->numberOfRows(), n->numberOfColumns(), false);
+  MatrixComplex<T> result = MatrixComplex<T>(operands, m.numberOfRows(), n.numberOfColumns());
   delete[] operands;
   return result;
 }
@@ -728,10 +724,10 @@ void Multiplication::addMissingFactors(Expression * factor, Context & context, A
   sortOperands(SimplificationOrder, false);
 }
 
-template Matrix * Multiplication::computeOnComplexAndMatrix<float>(Complex<float> const*, const Matrix*);
-template Matrix * Multiplication::computeOnComplexAndMatrix<double>(Complex<double> const*, const Matrix*);
-template Complex<float> Multiplication::compute<float>(Complex<float>, Complex<float>);
-template Complex<double> Multiplication::compute<double>(Complex<double>, Complex<double>);
+template MatrixComplex<float> Multiplication::computeOnComplexAndMatrix<float>(std::complex<float> const, const MatrixComplex<float>);
+template MatrixComplex<double> Multiplication::computeOnComplexAndMatrix<double>(std::complex<double> const, const MatrixComplex<double>);
+template std::complex<float> Multiplication::compute<float>(const std::complex<float>, const std::complex<float>);
+template std::complex<double> Multiplication::compute<double>(const std::complex<double>, const std::complex<double>);
 template void Multiplication::computeOnArrays<double>(double * m, double * n, double * result, int mNumberOfColumns, int mNumberOfRows, int nNumberOfColumns);
 
 }
diff --git a/poincare/src/naperian_logarithm.cpp b/poincare/src/naperian_logarithm.cpp
index 9a4ff62b0..4278c6465 100644
--- a/poincare/src/naperian_logarithm.cpp
+++ b/poincare/src/naperian_logarithm.cpp
@@ -37,11 +37,8 @@ Expression * NaperianLogarithm::shallowReduce(Context& context, AngleUnit angleU
 }
 
 template<typename T>
-Complex<T> NaperianLogarithm::computeOnComplex(const Complex<T> c, AngleUnit angleUnit) {
-  if (c.b() != 0) {
-    return Complex<T>::Float(NAN);
-  }
-  return Complex<T>::Float(std::log(c.a()));
+std::complex<T> NaperianLogarithm::computeOnComplex(const std::complex<T> c, AngleUnit angleUnit) {
+  return std::log(c);
 }
 
 }
diff --git a/poincare/src/nth_root.cpp b/poincare/src/nth_root.cpp
index 8f8005aec..22e738dba 100644
--- a/poincare/src/nth_root.cpp
+++ b/poincare/src/nth_root.cpp
@@ -1,5 +1,4 @@
 #include <poincare/nth_root.h>
-#include <poincare/complex.h>
 #include <poincare/division.h>
 #include <poincare/power.h>
 #include <poincare/undefined.h>
@@ -45,23 +44,14 @@ ExpressionLayout * NthRoot::privateCreateLayout(PrintFloat::Mode floatDisplayMod
 }
 
 template<typename T>
-Complex<T> NthRoot::compute(const Complex<T> c, const Complex<T> d) {
-  if (c.a() >= 0 && c.b() == 0 && d.b() == 0) {
-    return Complex<T>::Float(std::pow(c.a(), 1/d.a()));
-  }
-  Complex<T> invIndex = Division::compute(Complex<T>::Float(1), d);
-  return Power::compute(c, invIndex);
-}
-
-template<typename T>
-Expression * NthRoot::templatedApproximate(Context& context, AngleUnit angleUnit) const {
-  Expression * base = operand(0)->approximate<T>(context, angleUnit);
-  Expression * index = operand(1)->approximate<T>(context, angleUnit);
-  Complex<T> result = Complex<T>::Float(NAN);
-  if (base->type() == Type::Complex && index->type() == Type::Complex) {
+Evaluation<T> * NthRoot::templatedApproximate(Context& context, AngleUnit angleUnit) const {
+  Evaluation<T> * base = operand(0)->privateApproximate(T(), context, angleUnit);
+  Evaluation<T> * index = operand(1)->privateApproximate(T(), context, angleUnit);
+  Complex<T> result = Complex<T>::Undefined();
+  if (base->type() == Evaluation<T>::Type::Complex && index->type() == Evaluation<T>::Type::Complex) {
     Complex<T> * basec = static_cast<Complex<T> *>(base);
     Complex<T> * indexc = static_cast<Complex<T> *>(index);
-    result = compute(*basec, *indexc);
+    result = std::pow(*basec, std::complex<T>(1)/(*indexc));
   }
   delete base;
   delete index;
diff --git a/poincare/src/opposite.cpp b/poincare/src/opposite.cpp
index 249ad84ce..bd8e75ea8 100644
--- a/poincare/src/opposite.cpp
+++ b/poincare/src/opposite.cpp
@@ -2,7 +2,6 @@
 #include "layout/char_layout.h"
 #include "layout/horizontal_layout.h"
 #include <cmath>
-#include <poincare/complex.h>
 #include <poincare/layout_engine.h>
 #include <poincare/multiplication.h>
 #include <poincare/rational.h>
@@ -45,8 +44,8 @@ bool Opposite::needParenthesisWithParent(const Expression * e) const {
 }
 
 template<typename T>
-Complex<T> Opposite::compute(const Complex<T> c, AngleUnit angleUnit) {
-  return Complex<T>::Cartesian(-c.a(), -c.b());
+std::complex<T> Opposite::compute(const std::complex<T> c, AngleUnit angleUnit) {
+  return -c;
 }
 
 Expression * Opposite::shallowReduce(Context& context, AngleUnit angleUnit) {
@@ -94,5 +93,5 @@ int Opposite::writeTextInBuffer(char * buffer, int bufferSize, int numberOfSigni
 
 }
 
-template Poincare::Complex<float> Poincare::Opposite::compute<float>(Poincare::Complex<float>, AngleUnit angleUnit);
-template Poincare::Complex<double> Poincare::Opposite::compute<double>(Poincare::Complex<double>, AngleUnit angleUnit);
+template std::complex<float> Poincare::Opposite::compute<float>(const std::complex<float>, AngleUnit angleUnit);
+template std::complex<double> Poincare::Opposite::compute<double>(const std::complex<double>, AngleUnit angleUnit);
diff --git a/poincare/src/parenthesis.cpp b/poincare/src/parenthesis.cpp
index 09240a108..cada4ce7c 100644
--- a/poincare/src/parenthesis.cpp
+++ b/poincare/src/parenthesis.cpp
@@ -36,8 +36,8 @@ Expression * Parenthesis::shallowReduce(Context& context, AngleUnit angleUnit) {
 }
 
 template<typename T>
-Expression * Parenthesis::templatedApproximate(Context& context, AngleUnit angleUnit) const {
-  return operand(0)->approximate<T>(context, angleUnit);
+Evaluation<T> * Parenthesis::templatedApproximate(Context& context, AngleUnit angleUnit) const {
+  return operand(0)->privateApproximate(T(), context, angleUnit);
 }
 
 }
diff --git a/poincare/src/permute_coefficient.cpp b/poincare/src/permute_coefficient.cpp
index fd4536e97..bba8cd97b 100644
--- a/poincare/src/permute_coefficient.cpp
+++ b/poincare/src/permute_coefficient.cpp
@@ -68,30 +68,26 @@ Expression * PermuteCoefficient::shallowReduce(Context& context, AngleUnit angle
 }
 
 template<typename T>
-Complex<T> * PermuteCoefficient::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> * PermuteCoefficient::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;
   if (std::isnan(n) || std::isnan(k) || n != std::round(n) || k != std::round(k) || n < 0.0f || k < 0.0f) {
-    return new Complex<T>(Complex<T>::Float(NAN));
+    return new Complex<T>(Complex<T>::Undefined());
   }
   if (k > n) {
-    return new Complex<T>(Complex<T>::Float(0));
+    return new Complex<T>(0);
   }
   T result = 1;
   for (int i = (int)n-(int)k+1; i <= (int)n; i++) {
     result *= i;
     if (std::isinf(result) || std::isnan(result)) {
-      return new Complex<T>(Complex<T>::Float(result));
+      return new Complex<T>(result);
     }
   }
-  return new Complex<T>(Complex<T>::Float(std::round(result)));
+  return new Complex<T>(std::round(result));
 }
 
 }
diff --git a/poincare/src/power.cpp b/poincare/src/power.cpp
index b20c500c4..80df1944c 100644
--- a/poincare/src/power.cpp
+++ b/poincare/src/power.cpp
@@ -118,58 +118,47 @@ Expression * Power::setSign(Sign s, Context & context, AngleUnit angleUnit) {
 }
 
 template<typename T>
-Complex<T> Power::compute(const Complex<T> c, const Complex<T> d) {
-  // c == c.r * e^(c.th*i)
-  // d == d.a + d.b*i
-  // c^d == e^(ln(c^d))
-  //     == e^(ln(c) * d)
-  //     == e^(ln(c.r * e^(c.th*i))  *  (d.a + d.b*i))
-  //     == e^((ln(c.r) + ln(e^(c.th*i)))  *  (d.a + d.b*i))
-  //     == e^((ln(c.r) + c.th*i)  *  (d.a + d.b*i))
-  //     == e^(ln(c.r)*d.a + ln(c.r)*d.b*i + c.th*i*d.a + c.th*i*d.b*i)
-  //     == e^((ln(c.r^d.a) + ln(e^(c.th*d.b*i^2)))  +  (ln(c.r)*d.b + c.th*d.a)*i)
-  //     == e^(ln(c.r^d.a * e^(-c.th*d.b)))  *  e^((ln(c.r)*d.b + c.th*d.a)*i)
-  //     == c.r^d.a*e^(-c.th*d.b) * e^((ln(c.r)*d.b + c.th*d.a)*i)
-  if (c.a() == 0 && c.b() == 0) { // ln(c.r) and c.th are undefined
-    return Complex<T>::Float(d.a() > 0 ? 0 : NAN);
-  }
-  T radius = std::pow(c.r(), d.a()) * std::exp(-c.th() * d.b());
-  T theta = std::log(c.r())*d.b() + c.th()*d.a();
-  return Complex<T>::Polar(radius, theta);
+std::complex<T> Power::compute(const std::complex<T> c, const std::complex<T> d) {
+  return std::pow(c, d);
 }
 
-template<typename T> Matrix * Power::computeOnMatrixAndComplex(const Matrix * m, const Complex<T> * d) {
- if (m->numberOfRows() != m->numberOfColumns()) {
-    return nullptr;
+template<typename T> MatrixComplex<T> Power::computeOnComplexAndMatrix(const std::complex<T> c, const MatrixComplex<T> n) {
+  return MatrixComplex<T>::Undefined();
+}
+
+template<typename T> MatrixComplex<T> Power::computeOnMatrixAndComplex(const MatrixComplex<T> m, const std::complex<T> d) {
+ if (m.numberOfRows() != m.numberOfColumns()) {
+    return MatrixComplex<T>::Undefined();
   }
-  T power = d->toScalar();
+  T power = Complex<T>(d).toScalar();
   if (std::isnan(power) || std::isinf(power) || power != (int)power || std::fabs(power) > k_maxApproximatePowerMatrix) {
-    return nullptr;
+    return MatrixComplex<T>::Undefined();
   }
   if (power < 0) {
-    Matrix * inverse = m->createApproximateInverse<T>();
+    MatrixComplex<T> * inverse = m.createInverse();
     if (inverse == nullptr) {
-      return nullptr;
+      return MatrixComplex<T>::Undefined();
     }
-    Complex<T> minusC = Opposite::compute(*d, AngleUnit::Default);
-    Matrix * result = Power::computeOnMatrixAndComplex(inverse, &minusC);
+    Complex<T> minusC = Complex<T>(Opposite::compute(d, AngleUnit::Default));
+    MatrixComplex<T> result = Power::computeOnMatrixAndComplex(*inverse, minusC);
     delete inverse;
     return result;
   }
-  Matrix * result = Matrix::createApproximateIdentity<T>(m->numberOfRows());
+  MatrixComplex<T> result = MatrixComplex<T>::createIdentity(m.numberOfRows());
   // TODO: implement a quick exponentiation
   for (int k = 0; k < (int)power; k++) {
     if (shouldStopProcessing()) {
-      delete result;
-      return nullptr;
+      return MatrixComplex<T>::Undefined();
     }
-    Matrix * mult = Multiplication::computeOnMatrices<T>(result, m);
-    delete result;
-    result = mult;
+    result = Multiplication::computeOnMatrices<T>(result, m);
   }
   return result;
 }
 
+template<typename T> MatrixComplex<T> Power::computeOnMatrices(const MatrixComplex<T> m, const MatrixComplex<T> n) {
+  return MatrixComplex<T>::Undefined();
+}
+
 bool Power::needParenthesisWithParent(const Expression * e) const {
   Type types[] = {Type::Power, Type::Factorial};
   return e->isOfType(types, 2);
@@ -257,14 +246,14 @@ Expression * Power::shallowReduce(Context& context, AngleUnit angleUnit) {
 
   /* Step 0: if both operands are true complexes, the result is undefined.
    * We can assert that evaluations is a complex as matrix are not simplified */
-  Complex<float> * op0 = static_cast<Complex<float> *>(operand(0)->approximate<float>(context, angleUnit));
-  Complex<float> * op1 = static_cast<Complex<float> *>(operand(1)->approximate<float>(context, angleUnit));
-  bool bothOperandsComplexes = op0->b() != 0 && op1->b() != 0;
-  bool nonComplexNegativeOperand0 = op0->b() == 0 && op0->a() < 0;
+  Complex<float> * op0 = static_cast<Complex<float> *>(operand(0)->privateApproximate(float(), context, angleUnit));
+  Complex<float> * op1 = static_cast<Complex<float> *>(operand(1)->privateApproximate(float(), context, angleUnit));
+  bool bothOperandsComplexes = op0->imag() != 0 && op1->imag() != 0;
+  bool nonComplexNegativeOperand0 = op0->imag() == 0 && op0->real() < 0;
   delete op0;
   delete op1;
   if (bothOperandsComplexes) {
-    return replaceWith(new Undefined(), true);
+    return this;
   }
 
   /* Step 1: We handle simple cases as x^0, x^1, 0^x and 1^x first for 2 reasons:
@@ -893,6 +882,6 @@ bool Power::RationalExponentShouldNotBeReduced(const Rational * b, const Rationa
   return false;
 }
 
-template Complex<float> Power::compute<float>(Complex<float>, Complex<float>);
-template Complex<double> Power::compute<double>(Complex<double>, Complex<double>);
+template std::complex<float> Power::compute<float>(std::complex<float>, std::complex<float>);
+template std::complex<double> Power::compute<double>(std::complex<double>, std::complex<double>);
 }
diff --git a/poincare/src/prediction_interval.cpp b/poincare/src/prediction_interval.cpp
index 45c53ed61..88116abfb 100644
--- a/poincare/src/prediction_interval.cpp
+++ b/poincare/src/prediction_interval.cpp
@@ -5,7 +5,6 @@
 #include <poincare/power.h>
 #include <poincare/undefined.h>
 #include <poincare/division.h>
-#include <poincare/complex.h>
 extern "C" {
 #include <assert.h>
 }
@@ -72,23 +71,20 @@ Expression * PredictionInterval::shallowReduce(Context& context, AngleUnit angle
 }
 
 template<typename T>
-Expression * PredictionInterval::templatedApproximate(Context& context, AngleUnit angleUnit) const {
-  Expression * pInput = operand(0)->approximate<T>(context, angleUnit);
-  Expression * nInput = operand(1)->approximate<T>(context, angleUnit);
-  if (pInput->type() != Type::Complex || nInput->type() != Type::Complex) {
-    return new Complex<T>(Complex<T>::Float(NAN));
-  }
+Evaluation<T> * PredictionInterval::templatedApproximate(Context& context, AngleUnit angleUnit) const {
+  Evaluation<T> * pInput = operand(0)->privateApproximate(T(), context, angleUnit);
+  Evaluation<T> * nInput = operand(1)->privateApproximate(T(), context, angleUnit);
   T p = static_cast<Complex<T> *>(pInput)->toScalar();
   T n = static_cast<Complex<T> *>(nInput)->toScalar();
   delete pInput;
   delete nInput;
   if (std::isnan(p) || std::isnan(n) || n != (int)n || n < 0 || p < 0 || p > 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(p - 1.96*std::sqrt(p*(1.0-p))/std::sqrt(n)));
-  operands[1] = new Complex<T>(Complex<T>::Float(p + 1.96*std::sqrt(p*(1.0-p))/std::sqrt(n)));
-  return new Matrix(operands, 1, 2, false);
+  std::complex<T> operands[2];
+  operands[0] = std::complex<T>(p - 1.96*std::sqrt(p*(1.0-p))/std::sqrt(n));
+  operands[1] = std::complex<T>(p + 1.96*std::sqrt(p*(1.0-p))/std::sqrt(n));
+  return new MatrixComplex<T>(operands, 1, 2);
 }
 
 }
diff --git a/poincare/src/preferences.cpp b/poincare/src/preferences.cpp
index 2dbc12bce..c43eae780 100644
--- a/poincare/src/preferences.cpp
+++ b/poincare/src/preferences.cpp
@@ -1,5 +1,4 @@
 #include <poincare/preferences.h>
-#include <poincare/complex.h>
 
 namespace Poincare {
 
diff --git a/poincare/src/product.cpp b/poincare/src/product.cpp
index ab07fbc0d..89a8869d4 100644
--- a/poincare/src/product.cpp
+++ b/poincare/src/product.cpp
@@ -31,23 +31,23 @@ ExpressionLayout * Product::createSequenceLayoutWithArgumentLayouts(ExpressionLa
 }
 
 template<typename T>
-Expression * Product::templatedApproximateWithNextTerm(Expression * a, Expression * b) const {
-  if (a->type() == Type::Complex && b->type() == Type::Complex) {
+Evaluation<T> * Product::templatedApproximateWithNextTerm(Evaluation<T> * a, Evaluation<T> * b) const {
+  if (a->type() == Evaluation<T>::Type::Complex && b->type() == Evaluation<T>::Type::Complex) {
     Complex<T> * c = static_cast<Complex<T> *>(a);
     Complex<T> * d = static_cast<Complex<T> *>(b);
-    return new Complex<T>(Multiplication::compute(*c, *d));
+    return new Complex<T>((*c)*(*d));
   }
-  if (a->type() == Type::Complex) {
+  if (a->type() == Evaluation<T>::Type::Complex) {
     Complex<T> * c = static_cast<Complex<T> *>(a);
-    assert(b->type() == Type::Matrix);
-    Matrix * m = static_cast<Matrix *>(b);
-    return Multiplication::computeOnComplexAndMatrix(c, m);
+    assert(b->type() == Evaluation<T>::Type::MatrixComplex);
+    MatrixComplex<T> * m = static_cast<MatrixComplex<T> *>(b);
+    return new MatrixComplex<T>(Multiplication::computeOnComplexAndMatrix(*c, *m));
   }
-  assert(a->type() == Type::Matrix);
-  assert(b->type() == Type::Matrix);
-  Matrix * m = static_cast<Matrix *>(a);
-  Matrix * n = static_cast<Matrix *>(b);
-  return Multiplication::computeOnMatrices<T>(m, n);
+  assert(a->type() == Evaluation<T>::Type::MatrixComplex);
+  assert(b->type() == Evaluation<T>::Type::MatrixComplex);
+  MatrixComplex<T> * m = static_cast<MatrixComplex<T> *>(a);
+  MatrixComplex<T> * n = static_cast<MatrixComplex<T> *>(b);
+  return new MatrixComplex<T>(Multiplication::computeOnMatrices<T>(*m, *n));
 }
 
 }
diff --git a/poincare/src/randint.cpp b/poincare/src/randint.cpp
index 2cdec3665..b94bb0636 100644
--- a/poincare/src/randint.cpp
+++ b/poincare/src/randint.cpp
@@ -18,22 +18,19 @@ Expression * Randint::clone() const {
   return a;
 }
 
-template <typename T> Expression * Randint::templateApproximate(Context & context, AngleUnit angleUnit) const {
-  Expression * aInput = operand(0)->approximate<T>(context, angleUnit);
-  Expression * bInput = operand(1)->approximate<T>(context, angleUnit);
-  if (aInput->type() != Type::Complex || bInput->type() != Type::Complex) {
-    return new Complex<T>(Complex<T>::Float(NAN));
-  }
-  T a = static_cast<Complex<T> *>(aInput)->toScalar();
-  T b = static_cast<Complex<T> *>(bInput)->toScalar();
+template <typename T> Evaluation<T> * Randint::templateApproximate(Context & context, AngleUnit angleUnit) const {
+  Evaluation<T> * aInput = operand(0)->privateApproximate(T(), context, angleUnit);
+  Evaluation<T> * bInput = operand(1)->privateApproximate(T(), context, angleUnit);
+  T a = aInput->toScalar();
+  T b = bInput->toScalar();
   delete aInput;
   delete bInput;
   if (std::isnan(a) || std::isnan(b) || a != std::round(a) || b != std::round(b) || a > b) {
-    return new Complex<T>(Complex<T>::Float(NAN));
+    return new Complex<T>(Complex<T>::Undefined());
 
   }
   T result = std::floor(Random::random<T>()*(b+1.0-a)+a);
-  return new Complex<T>(Complex<T>::Float(result));
+  return new Complex<T>(result);
 }
 
 }
diff --git a/poincare/src/rational.cpp b/poincare/src/rational.cpp
index 03e9c1954..c5ae52a9c 100644
--- a/poincare/src/rational.cpp
+++ b/poincare/src/rational.cpp
@@ -140,7 +140,7 @@ int Rational::simplificationOrderSameType(const Expression * e, bool canBeInterr
 template<typename T> Complex<T> * Rational::templatedApproximate(Context& context, Expression::AngleUnit angleUnit) const {
   T n = m_numerator.approximate<T>();
   T d = m_denominator.approximate<T>();
-  return new Complex<T>(Complex<T>::Float(n/d));
+  return new Complex<T>(n/d);
 }
 
 bool Rational::needParenthesisWithParent(const Expression * e) const {
diff --git a/poincare/src/real_part.cpp b/poincare/src/real_part.cpp
index aa8d4904f..a600cd737 100644
--- a/poincare/src/real_part.cpp
+++ b/poincare/src/real_part.cpp
@@ -1,5 +1,4 @@
 #include <poincare/real_part.h>
-#include <poincare/complex.h>
 #include <poincare/simplification_engine.h>
 extern "C" {
 #include <assert.h>
@@ -35,8 +34,8 @@ Expression * RealPart::shallowReduce(Context& context, AngleUnit angleUnit) {
 }
 
 template<typename T>
-Complex<T> RealPart::computeOnComplex(const Complex<T> c, AngleUnit angleUnit) {
-  return Complex<T>::Float(c.a());
+std::complex<T> RealPart::computeOnComplex(const std::complex<T> c, AngleUnit angleUnit) {
+  return std::real(c);
 }
 
 }
diff --git a/poincare/src/round.cpp b/poincare/src/round.cpp
index ee38ad9a8..7a04f667e 100644
--- a/poincare/src/round.cpp
+++ b/poincare/src/round.cpp
@@ -54,17 +54,17 @@ Expression * Round::shallowReduce(Context& context, AngleUnit angleUnit) {
 
 template<typename T>
 Complex<T> * Round::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(f2) || f2 != std::round(f2)) {
-    return new Complex<T>(Complex<T>::Float(NAN));
+    return new Complex<T>(Complex<T>::Undefined());
   }
   T err = std::pow(10, std::floor(f2));
-  return new Complex<T>(Complex<T>::Float(std::round(f1*err)/err));
+  return new Complex<T>(std::round(f1*err)/err);
 }
 
 }
diff --git a/poincare/src/sequence.cpp b/poincare/src/sequence.cpp
index 07aeb5159..ee67259e4 100644
--- a/poincare/src/sequence.cpp
+++ b/poincare/src/sequence.cpp
@@ -1,5 +1,5 @@
 #include <poincare/sequence.h>
-#include <poincare/complex.h>
+#include <poincare/decimal.h>
 #include <poincare/symbol.h>
 #include <poincare/undefined.h>
 #include <poincare/variable_context.h>
@@ -18,33 +18,31 @@ ExpressionLayout * Sequence::privateCreateLayout(PrintFloat::Mode floatDisplayMo
 }
 
 template<typename T>
-Expression * Sequence::templatedApproximate(Context& context, AngleUnit angleUnit) const {
-  Expression * aInput = operand(1)->approximate<T>(context, angleUnit);
-  Expression * bInput = operand(2)->approximate<T>(context, angleUnit);
-  T start = aInput->type() == Type::Complex ? static_cast<Complex<T> *>(aInput)->toScalar() : NAN;
-  T end = bInput->type() == Type::Complex ? static_cast<Complex<T> *>(bInput)->toScalar() : NAN;
+Evaluation<T> * Sequence::templatedApproximate(Context& context, AngleUnit angleUnit) const {
+  Evaluation<T> * aInput = operand(1)->privateApproximate(T(), context, angleUnit);
+  Evaluation<T> * bInput = operand(2)->privateApproximate(T(), context, angleUnit);
+  T start = aInput->toScalar();
+  T end = bInput->toScalar();
   delete aInput;
   delete bInput;
   if (std::isnan(start) || std::isnan(end) || start != (int)start || end != (int)end || end - start > k_maxNumberOfSteps) {
-    return new Complex<T>(Complex<T>::Float(NAN));
+    return new Complex<T>(Complex<T>::Undefined());
   }
-  VariableContext<T> nContext = VariableContext<T>('n', &context);
-  Symbol nSymbol('n');
-  Expression * result = new Complex<T>(Complex<T>::Float(emptySequenceValue()));
+  VariableContext nContext = VariableContext('n', &context);
+  Evaluation<T> * result = new Complex<T>(emptySequenceValue());
   for (int i = (int)start; i <= (int)end; i++) {
     if (shouldStopProcessing()) {
       delete result;
-      return new Complex<T>(Complex<T>::Float(NAN));
+      return new Complex<T>(Complex<T>::Undefined());
     }
-    Complex<T> iExpression = Complex<T>::Float(i);
-    nContext.setExpressionForSymbolName(&iExpression, &nSymbol, nContext);
-    Expression * expression = operand(0)->approximate<T>(nContext, angleUnit);
-    Expression * newResult = evaluateWithNextTerm(T(), result, expression);
+    nContext.setApproximationForVariable((T)i);
+    Evaluation<T> * expression = operand(0)->privateApproximate(T(), nContext, angleUnit);
+    Evaluation<T> * newResult = evaluateWithNextTerm(T(), result, expression);
     delete result;
     delete expression;
     result = newResult;
     if (result == nullptr) {
-      return new Complex<T>(Complex<T>::Float(NAN));
+      return new Complex<T>(Complex<T>::Undefined());
     }
   }
   return result;
diff --git a/poincare/src/sine.cpp b/poincare/src/sine.cpp
index 324dc02a0..ac325cf54 100644
--- a/poincare/src/sine.cpp
+++ b/poincare/src/sine.cpp
@@ -1,7 +1,6 @@
 #include <poincare/sine.h>
 #include <poincare/trigonometry.h>
 #include <poincare/hyperbolic_sine.h>
-#include <poincare/complex.h>
 #include <poincare/multiplication.h>
 #include <poincare/symbol.h>
 #include <poincare/simplification_engine.h>
@@ -40,24 +39,4 @@ Expression * Sine::shallowReduce(Context& context, AngleUnit angleUnit) {
   return Trigonometry::shallowReduceDirectFunction(this, context, angleUnit);
 }
 
-template<typename T>
-Complex<T> Sine::computeOnComplex(const Complex<T> c, AngleUnit angleUnit) {
-  if (c.b() == 0) {
-    T input = c.a();
-    if (angleUnit == AngleUnit::Degree) {
-      input *= M_PI/180;
-    }
-    T result = std::sin(input);
-    /* Cheat: see comment in cosine.cpp
-     * We cheat to avoid returning sin(Pi) = epsilon */
-    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());
-  Complex<T> sinh = HyperbolicSine::computeOnComplex(arg, angleUnit);
-  return Multiplication::compute(Complex<T>::Cartesian(0, -1), sinh);
-}
-
 }
diff --git a/poincare/src/square_root.cpp b/poincare/src/square_root.cpp
index 7b8105b7a..f74e53c03 100644
--- a/poincare/src/square_root.cpp
+++ b/poincare/src/square_root.cpp
@@ -1,5 +1,4 @@
 #include <poincare/square_root.h>
-#include <poincare/complex.h>
 #include <poincare/power.h>
 #include <poincare/simplification_engine.h>
 #include "layout/nth_root_layout.h"
@@ -26,11 +25,8 @@ int SquareRoot::writeTextInBuffer(char * buffer, int bufferSize, int numberOfSig
 }
 
 template<typename T>
-Complex<T> SquareRoot::computeOnComplex(const Complex<T> c, AngleUnit angleUnit) {
-  if (c.b() == 0 && c.a() >= 0) {
-    return Complex<T>::Float(std::sqrt(c.a()));
-  }
-  return Power::compute(c, Complex<T>::Float(0.5));
+std::complex<T> SquareRoot::computeOnComplex(const std::complex<T> c, AngleUnit angleUnit) {
+  return std::sqrt(c);
 }
 
 Expression * SquareRoot::shallowReduce(Context& context, AngleUnit angleUnit) {
diff --git a/poincare/src/store.cpp b/poincare/src/store.cpp
index 31164a008..6a16e3e12 100644
--- a/poincare/src/store.cpp
+++ b/poincare/src/store.cpp
@@ -5,7 +5,6 @@ extern "C" {
 }
 #include <poincare/store.h>
 #include <ion.h>
-#include <poincare/complex.h>
 #include <poincare/context.h>
 #include "layout/char_layout.h"
 #include "layout/horizontal_layout.h"
@@ -45,12 +44,12 @@ ExpressionLayout * Store::privateCreateLayout(PrintFloat::Mode floatDisplayMode,
 }
 
 template<typename T>
-Expression * Store::templatedApproximate(Context& context, AngleUnit angleUnit) const {
+Evaluation<T> * Store::templatedApproximate(Context& context, AngleUnit angleUnit) const {
   context.setExpressionForSymbolName(value(), symbol(), context);
   if (context.expressionForSymbol(symbol()) != nullptr) {
-    return context.expressionForSymbol(symbol())->approximate<T>(context, angleUnit);
+    return context.expressionForSymbol(symbol())->privateApproximate(T(), context, angleUnit);
   }
-  return new Complex<T>(Complex<T>::Float(NAN));
+  return new Complex<T>(Complex<T>::Undefined());
 }
 
 }
diff --git a/poincare/src/subtraction.cpp b/poincare/src/subtraction.cpp
index e517138f4..077fb1fff 100644
--- a/poincare/src/subtraction.cpp
+++ b/poincare/src/subtraction.cpp
@@ -39,23 +39,18 @@ bool Subtraction::needParenthesisWithParent(const Expression * e) const {
 }
 
 template<typename T>
-Complex<T> Subtraction::compute(const Complex<T> c, const Complex<T> d) {
-  return Complex<T>::Cartesian(c.a()-d.a(), c.b() - d.b());
+std::complex<T> Subtraction::compute(const std::complex<T> c, const std::complex<T> d) {
+  return c - d;
 }
 
-template<typename T> Matrix * Subtraction::computeOnComplexAndMatrix(const Complex<T> * c, const Matrix * m) {
-  Matrix * opposite = computeOnMatrixAndComplex(m, c);
-  if (opposite == nullptr) {
-    return nullptr;
+template<typename T> MatrixComplex<T> Subtraction::computeOnComplexAndMatrix(const std::complex<T> c, const MatrixComplex<T> m) {
+  MatrixComplex<T> opposite = computeOnMatrixAndComplex(m, c);
+  std::complex<T> * operands = new std::complex<T> [opposite.numberOfComplexOperands()];
+  for (int i = 0; i < opposite.numberOfComplexOperands(); i++) {
+    operands[i] = -opposite.complexOperand(i);
   }
-  Expression ** operands = new Expression * [opposite->numberOfRows() * opposite->numberOfColumns()];
-  for (int i = 0; i < opposite->numberOfOperands(); i++) {
-    const Complex<T> * entry = static_cast<const Complex<T> *>(opposite->operand(i));
-    operands[i] = new Complex<T>(Complex<T>::Cartesian(-entry->a(), -entry->b()));
-  }
-  Matrix * result = new Matrix(operands, m->numberOfRows(), m->numberOfColumns(), false);
+  MatrixComplex<T> result = MatrixComplex<T>(operands, opposite.numberOfRows(), opposite.numberOfColumns());
   delete[] operands;
-  delete opposite;
   return result;
 }
 
diff --git a/poincare/src/sum.cpp b/poincare/src/sum.cpp
index fe973047a..61369e352 100644
--- a/poincare/src/sum.cpp
+++ b/poincare/src/sum.cpp
@@ -31,23 +31,23 @@ ExpressionLayout * Sum::createSequenceLayoutWithArgumentLayouts(ExpressionLayout
 }
 
 template<typename T>
-Expression * Sum::templatedApproximateWithNextTerm(Expression * a, Expression * b) const {
-  if (a->type() == Type::Complex && b->type() == Type::Complex) {
+Evaluation<T> * Sum::templatedApproximateWithNextTerm(Evaluation<T> * a, Evaluation<T> * b) const {
+  if (a->type() == Evaluation<T>::Type::Complex && b->type() == Evaluation<T>::Type::Complex) {
     Complex<T> * c = static_cast<Complex<T> *>(a);
     Complex<T> * d = static_cast<Complex<T> *>(b);
-    return new Complex<T>(Addition::compute(*c, *d));
+    return new Complex<T>((*c)+(*d));
   }
-  if (a->type() == Type::Complex) {
+  if (a->type() == Evaluation<T>::Type::Complex) {
     Complex<T> * c = static_cast<Complex<T> *>(a);
-    assert(b->type() == Type::Matrix);
-    Matrix * m = static_cast<Matrix *>(b);
-    return Addition::computeOnComplexAndMatrix(c, m);
+    assert(b->type() == Evaluation<T>::Type::MatrixComplex);
+    MatrixComplex<T> * m = static_cast<MatrixComplex<T> *>(b);
+    return new MatrixComplex<T>(Addition::computeOnComplexAndMatrix(*c, *m));
   }
-  assert(a->type() == Type::Matrix);
-  assert(b->type() == Type::Matrix);
-  Matrix * m = static_cast<Matrix *>(a);
-  Matrix * n = static_cast<Matrix *>(b);
-  return Addition::computeOnMatrices<T>(m, n);
+  assert(a->type() == Evaluation<T>::Type::MatrixComplex);
+  assert(b->type() == Evaluation<T>::Type::MatrixComplex);
+  MatrixComplex<T> * m = static_cast<MatrixComplex<T> *>(a);
+  MatrixComplex<T> * n = static_cast<MatrixComplex<T> *>(b);
+  return new MatrixComplex<T>(Addition::computeOnMatrices<T>(*m, *n));
 }
 
 }
diff --git a/poincare/src/symbol.cpp b/poincare/src/symbol.cpp
index 21a349b6a..928ea5f7c 100644
--- a/poincare/src/symbol.cpp
+++ b/poincare/src/symbol.cpp
@@ -1,7 +1,6 @@
 #include <poincare/symbol.h>
 
 #include <poincare/context.h>
-#include <poincare/complex.h>
 #include <poincare/division.h>
 #include <poincare/multiplication.h>
 #include <poincare/parenthesis.h>
@@ -229,11 +228,14 @@ Expression * Symbol::shallowReduce(Context& context, AngleUnit angleUnit) {
 }
 
 template<typename T>
-Expression * Symbol::templatedApproximate(Context& context, AngleUnit angleUnit) const {
-  if (context.expressionForSymbol(this) != nullptr) {
-    return context.expressionForSymbol(this)->approximate<T>(context, angleUnit);
+Evaluation<T> * Symbol::templatedApproximate(Context& context, AngleUnit angleUnit) const {
+  if (m_name == Ion::Charset::IComplex) {
+    return new Complex<T>(0.0, 1.0);
   }
-  return new Complex<T>(Complex<T>::Float(NAN));
+  if (context.expressionForSymbol(this) != nullptr) {
+    return context.expressionForSymbol(this)->privateApproximate(T(), context, angleUnit);
+  }
+  return new Complex<T>(Complex<T>::Undefined());
 }
 
 Expression::Type Symbol::type() const {
diff --git a/poincare/src/tangent.cpp b/poincare/src/tangent.cpp
index ef97c8741..99c627c4b 100644
--- a/poincare/src/tangent.cpp
+++ b/poincare/src/tangent.cpp
@@ -1,5 +1,4 @@
 #include <poincare/tangent.h>
-#include <poincare/complex.h>
 #include <poincare/sine.h>
 #include <poincare/cosine.h>
 #include <poincare/division.h>
@@ -50,14 +49,4 @@ Expression * Tangent::shallowReduce(Context& context, AngleUnit angleUnit) {
   return newExpression;
 }
 
-template<typename T>
-Complex<T> Tangent::computeOnComplex(const Complex<T> c, AngleUnit angleUnit) {
-  Complex<T> result = Division::compute(Sine::computeOnComplex(c, angleUnit), Cosine::computeOnComplex(c, angleUnit));
-  if (!std::isnan(result.a()) && !std::isnan(result.b())) {
-    return result;
-  }
-  Complex<T> tanh = HyperbolicTangent::computeOnComplex(Multiplication::compute(Complex<T>::Cartesian(0, -1), c), angleUnit);
-  return Multiplication::compute(Complex<T>::Cartesian(0, 1), tanh);
-}
-
 }
diff --git a/poincare/src/trigonometry.cpp b/poincare/src/trigonometry.cpp
index f6e8fccc0..ff40d8e2c 100644
--- a/poincare/src/trigonometry.cpp
+++ b/poincare/src/trigonometry.cpp
@@ -1,12 +1,13 @@
 #include <poincare/trigonometry.h>
 #include <poincare/hyperbolic_cosine.h>
-#include <poincare/complex.h>
 #include <poincare/symbol.h>
+#include <poincare/preferences.h>
 #include <poincare/undefined.h>
 #include <poincare/rational.h>
 #include <poincare/multiplication.h>
 #include <poincare/subtraction.h>
 #include <poincare/derivative.h>
+#include <poincare/decimal.h>
 #include <ion.h>
 extern "C" {
 #include <assert.h>
@@ -35,7 +36,7 @@ float Trigonometry::characteristicXRange(const Expression * e, Context & context
   assert(d == 1);
   /* To compute a, the slope of the expression op, we compute the derivative of
    * op for any x value. */
-  Poincare::Complex<float> x = Poincare::Complex<float>::Float(1.0f);
+  Poincare::Decimal x(1.0f);
   const Poincare::Expression * args[2] = {op, &x};
   Poincare::Derivative derivative(args, true);
   float a = derivative.approximateToScalar<float>(context);
@@ -240,4 +241,31 @@ Expression * Trigonometry::table(const Expression * e, Expression::Type type, Co
   return nullptr;
 }
 
+template <typename T>
+std::complex<T> Trigonometry::computeOnComplex(const std::complex<T> c, Expression::AngleUnit angleUnit, Approximation<T> approximate) {
+  assert(angleUnit != Expression::AngleUnit::Default);
+  std::complex<T> input(c);
+  if (angleUnit == Expression::AngleUnit::Degree && input.imag() == 0.0) {
+    input = input*std::complex<T>(M_PI/180.0);
+  }
+  std::complex<T> result = approximate(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 (std::abs(input) !=  0 && std::fabs(result.real()/std::abs(input)) <= Expression::epsilon<T>()) {
+    result.real(0);
+  }
+  if (std::abs(input) !=  0 && std::fabs(result.imag()/std::abs(input)) <= Expression::epsilon<T>()) {
+    result.imag(0);
+  }
+  return result;
+}
+
+template std::complex<double> Trigonometry::computeOnComplex<double>(const std::complex<double>, Expression::AngleUnit, Approximation<double>);
+template std::complex<float> Trigonometry::computeOnComplex<float>(const std::complex<float>, Expression::AngleUnit, Approximation<float>);
+
 }
diff --git a/poincare/src/undefined.cpp b/poincare/src/undefined.cpp
index c5600b15c..aefe6c373 100644
--- a/poincare/src/undefined.cpp
+++ b/poincare/src/undefined.cpp
@@ -21,7 +21,7 @@ int Undefined::polynomialDegree(char symbolName) const {
 }
 
 template<typename T> Complex<T> * Undefined::templatedApproximate(Context& context, AngleUnit angleUnit) const {
-  return new Complex<T>(Complex<T>::Float(NAN));
+  return new Complex<T>(Complex<T>::Undefined());
 }
 
 ExpressionLayout * Undefined::privateCreateLayout(PrintFloat::Mode floatDisplayMode, ComplexFormat complexFormat) const {
diff --git a/poincare/src/variable_context.cpp b/poincare/src/variable_context.cpp
index a85ff9072..f89686d66 100644
--- a/poincare/src/variable_context.cpp
+++ b/poincare/src/variable_context.cpp
@@ -1,44 +1,58 @@
 #include <poincare/variable_context.h>
 #include <poincare/preferences.h>
+#include <poincare/undefined.h>
+#include <poincare/expression.h>
 #include <assert.h>
 #include <cmath>
 
 namespace Poincare {
 
-template<typename T>
-VariableContext<T>::VariableContext(char name, Context * parentContext) :
+VariableContext::VariableContext(char name, Context * parentContext) :
   m_name(name),
-  m_value(Complex<T>::Float(NAN)),
+  m_value(nullptr),
   m_parentContext(parentContext)
 {
 }
 
-template<typename T>
-void VariableContext<T>::setExpressionForSymbolName(const Expression * expression, const Symbol * symbol, Context & context) {
+VariableContext::~VariableContext() {
+  if (m_value) {
+    delete m_value;
+    m_value = nullptr;
+  }
+}
+
+template<typename T> void VariableContext::setApproximationForVariable(T value) {
+  if (m_value) {
+    delete m_value;
+    m_value = nullptr;
+  }
+  m_value = Expression::CreateDecimal<T>(value);
+}
+
+void VariableContext::setExpressionForSymbolName(const Expression * expression, const Symbol * symbol, Context & context) {
   if (symbol->name() == m_name) {
     if (expression == nullptr) {
       return;
     }
-    if (expression->type() == Expression::Type::Complex) {
-      m_value = Complex<T>::Float(static_cast<const Complex<T> *>(expression)->toScalar());
-    } else {
-      m_value = Complex<T>::Float(NAN);
+    if (m_value) {
+      delete m_value;
+      m_value = nullptr;
     }
+    assert(false); // TODO: remove
+    m_value = expression->clone();
   } else {
     m_parentContext->setExpressionForSymbolName(expression, symbol, context);
   }
 }
 
-template<typename T>
-const Expression * VariableContext<T>::expressionForSymbol(const Symbol * symbol) {
+const Expression * VariableContext::expressionForSymbol(const Symbol * symbol) {
   if (symbol->name() == m_name) {
-    return &m_value;
+    return m_value;
   } else {
     return m_parentContext->expressionForSymbol(symbol);
   }
 }
 
-template class Poincare::VariableContext<float>;
-template class Poincare::VariableContext<double>;
-
+template void VariableContext::setApproximationForVariable<double>(double);
+template void VariableContext::setApproximationForVariable<float>(float);
 }
diff --git a/poincare/test/helper.cpp b/poincare/test/helper.cpp
index 7100b2461..45913138c 100644
--- a/poincare/test/helper.cpp
+++ b/poincare/test/helper.cpp
@@ -70,6 +70,7 @@ template<typename T>
 void assert_parsed_expression_evaluates_to(const char * expression, Complex<T> * results, int numberOfRows, int numberOfColumns, Expression::AngleUnit angleUnit) {
   GlobalContext globalContext;
   Expression * a = parse_expression(expression);
+  Expression::Simplify(&a, globalContext, angleUnit);
   Expression * m = a->approximate<T>(globalContext, angleUnit);
   assert(m);
   assert(m->numberOfOperands() == 0 || m->numberOfOperands() == numberOfRows*numberOfColumns);