diff --git a/poincare/Makefile b/poincare/Makefile
index f5582515b..cc77414fa 100644
--- a/poincare/Makefile
+++ b/poincare/Makefile
@@ -3,6 +3,7 @@ SFLAGS += -Ipoincare/include
 #include poincare/src/simplify/Makefile
 #include poincare/src/simplification/Makefile
 objs += $(addprefix poincare/src/,\
+  addition.o\
   matrix_complex.o\
   tree_node.o\
   tree_pool.o\
diff --git a/poincare/include/poincare/addition.h b/poincare/include/poincare/addition.h
index 73c9b8eaf..c725c51d3 100644
--- a/poincare/include/poincare/addition.h
+++ b/poincare/include/poincare/addition.h
@@ -1,27 +1,41 @@
 #ifndef POINCARE_ADDITION_H
 #define POINCARE_ADDITION_H
 
-#include <poincare/dynamic_hierarchy.h>
-#include <poincare/rational.h>
-#include <poincare/layout_helper.h>
 #include <poincare/approximation_helper.h>
+#include <poincare/layout_helper.h>
+#include <poincare/n_ary_expression_node.h>
+#include <poincare/rational.h>
+#include <poincare/serialization_helper.h>
 
 namespace Poincare {
 
-class Addition : public DynamicHierarchy {
-  using DynamicHierarchy::DynamicHierarchy;
-  friend class Logarithm;
+class AdditionNode : public NAryExpressionNode {
+/* TODO friend class Logarithm;
   friend class Multiplication;
   friend class Subtraction;
   friend class Power;
   friend class Complex<float>;
-  friend class Complex<double>;
+  friend class Complex<double>;*/
 public:
-  Type type() const override;
+  using NAryExpressionNode::NAryExpressionNode;
+
+  // Tree
+  static AdditionNode * FailedAllocationStaticNode();
+  AdditionNode * failedAllocationStaticNode() override { return FailedAllocationStaticNode(); }
+  size_t size() const override { return sizeof(AdditionNode); }
+#if TREE_LOG
+  const char * description() const override { return "Addition";  }
+#endif
+
+  // ExpressionNode
+
+  // Properties
+  Type type() const override { return Type::Addition; }
   int polynomialDegree(char symbolName) const override;
-  int privateGetPolynomialCoefficients(char symbolName, Expression * coefficients[]) const override;
-  /* Evaluation */
-  template<typename T> static std::complex<T> compute(const std::complex<T> c, const std::complex<T> d);
+  int getPolynomialCoefficients(char symbolName, Expression coefficients[]) const override;
+
+  // Evaluation
+  template<typename T> static Complex<T> compute(const std::complex<T> c, const std::complex<T> d) { return Complex<T>(c+d); }
   template<typename T> static MatrixComplex<T> computeOnMatrices(const MatrixComplex<T> m, const MatrixComplex<T> n) {
     return ApproximationHelper::ElementWiseOnComplexMatrices(m, n, compute<T>);
   }
@@ -29,23 +43,23 @@ public:
     return ApproximationHelper::ElementWiseOnMatrixComplexAndComplex(m, c, compute<T>);
   }
 private:
-  /* Layout */
-  bool needsParenthesesWithParent(SerializableNode * parentNode) const override;
+  // Layout
+  bool needsParenthesesWithParent(const SerializationHelperInterface * parentNode) const override;
   LayoutRef createLayout(Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const override {
-    return LayoutHelper::Infix(this, floatDisplayMode, numberOfSignificantDigits, name());
+    return LayoutHelper::Infix(Expression(this), floatDisplayMode, numberOfSignificantDigits, name());
   }
   int serialize(char * buffer, int bufferSize, Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const override {
-    return LayoutHelper::writeInfixExpressionTextInBuffer(this, buffer, bufferSize, floatDisplayMode, numberOfSignificantDigits, name());
+    return SerializationHelper::Infix(this, buffer, bufferSize, floatDisplayMode, numberOfSignificantDigits, name());
   }
   static const char * name() { return "+"; }
 
-  /* Simplification */
-  Expression shallowReduce(Context& context, Preferences::AngleUnit angleUnit) override;
-  Expression * shallowBeautify(Context & context, Preferences::AngleUnit angleUnit) override;
-  Expression * factorizeOnCommonDenominator(Context & context, Preferences::AngleUnit angleUnit);
-  void factorizeOperands(Expression * e1, Expression * e2, Context & context, Preferences::AngleUnit angleUnit);
-  static const Rational RationalFactor(Expression * e);
-  static bool TermsHaveIdenticalNonRationalFactors(const Expression * e1, const Expression * e2);
+  // Simplification
+  Expression shallowReduce(Context& context, Preferences::AngleUnit angleUnit) const override;
+  Expression shallowBeautify(Context & context, Preferences::AngleUnit angleUnit) const override;
+  Expression factorizeOnCommonDenominator(Context & context, Preferences::AngleUnit angleUnit) const;
+  void factorizeOperands(Expression e1, Expression e2, Context & context, Preferences::AngleUnit angleUnit);
+  static const Rational RationalFactor(Expression e);
+  static bool TermsHaveIdenticalNonRationalFactors(const Expression e1, const Expression e2);
   /* Evaluation */
   template<typename T> static MatrixComplex<T> computeOnMatrixAndComplex(const MatrixComplex<T> m, const std::complex<T> c) {
     return ApproximationHelper::ElementWiseOnMatrixComplexAndComplex(m, c, compute<T>);
@@ -58,6 +72,16 @@ private:
    }
 };
 
+class Addition : public NAryExpression {
+public:
+  Addition(const AdditionNode * n) : NAryExpression(n) {}
+  Addition() : NAryExpression(TreePool::sharedPool()->createTreeNode<AdditionNode>()) {}
+
+  // Expression
+  Expression shallowReduce(Context& context, Preferences::AngleUnit angleUnit) const;
+  Expression shallowBeautify(Context& context, Preferences::AngleUnit angleUnit) const;
+};
+
 }
 
 #endif
diff --git a/poincare/include/poincare/expression.h b/poincare/include/poincare/expression.h
index 201efae79..53e6636a3 100644
--- a/poincare/include/poincare/expression.h
+++ b/poincare/include/poincare/expression.h
@@ -13,6 +13,7 @@ namespace Poincare {
 class Context;
 
 class Expression : public TreeByValue {
+  friend class AdditionNode;
   friend class ExpressionNode;
   friend class NAryExpressionNode;
   friend class SubtractionNode;
diff --git a/poincare/include/poincare/n_ary_expression_node.h b/poincare/include/poincare/n_ary_expression_node.h
index b1c69c253..cf3aef8e7 100644
--- a/poincare/include/poincare/n_ary_expression_node.h
+++ b/poincare/include/poincare/n_ary_expression_node.h
@@ -33,6 +33,7 @@ private:
 
 class NAryExpression : public Expression {
 public:
+  NAryExpression(const NAryExpressionNode * n) : Expression(n) {}
   void addChildAtIndexInPlace(TreeByValue t, int index, int currentNumberOfChildren) {
     Expression::addChildAtIndexInPlace(t, index, currentNumberOfChildren);
   }
diff --git a/poincare/src/addition.cpp b/poincare/src/addition.cpp
index d8f021b50..9c1f7d5a5 100644
--- a/poincare/src/addition.cpp
+++ b/poincare/src/addition.cpp
@@ -1,32 +1,23 @@
 #include <poincare/addition.h>
-#include <poincare/multiplication.h>
+//#include <poincare/multiplication.h>
 #include <poincare/subtraction.h>
-#include <poincare/power.h>
+//#include <poincare/power.h>
 #include <poincare/opposite.h>
 #include <poincare/undefined.h>
-#include <poincare/matrix.h>
-extern "C" {
+//#include <poincare/matrix.h>
 #include <assert.h>
-#include <stdlib.h>
-}
 
 namespace Poincare {
 
-Expression::Type Addition::type() const {
-  return Type::Addition;
+static AdditionNode * FailedAllocationStaticNode() {
+  static AllocationFailureExpressionNode<AdditionNode> failure;
+  return &failure;
 }
 
-Expression * Addition::clone() const {
-  if (numberOfChildren() == 0) {
-    return new Addition();
-  }
-  return new Addition(operands(), numberOfChildren(), true);
-}
-
-int Addition::polynomialDegree(char symbolName) const {
+int AdditionNode::polynomialDegree(char symbolName) const {
   int degree = 0;
   for (int i = 0; i < numberOfChildren(); i++) {
-    int d = operand(i)->polynomialDegree(symbolName);
+    int d = childAtIndex(i)->polynomialDegree(symbolName);
     if (d < 0) {
       return -1;
     }
@@ -35,82 +26,277 @@ int Addition::polynomialDegree(char symbolName) const {
   return degree;
 }
 
-int Addition::getPolynomialCoefficients(char symbolName, Expression coefficients[]) const {
+int AdditionNode::getPolynomialCoefficients(char symbolName, Expression coefficients[]) const {
   int deg = polynomialDegree(symbolName);
-  if (deg < 0 || deg > k_maxPolynomialDegree) {
+  if (deg < 0 || deg > Expression::k_maxPolynomialDegree) {
     return -1;
   }
   for (int k = 0; k < deg+1; k++) {
-    coefficients[k] = new Addition();
+    coefficients[k] = Addition();
   }
-  Expression * intermediateCoefficients[k_maxNumberOfPolynomialCoefficients];
+  Expression intermediateCoefficients[Expression::k_maxNumberOfPolynomialCoefficients];
   for (int i = 0; i < numberOfChildren(); i++) {
-    int d = operand(i)->privateGetPolynomialCoefficients(symbolName, intermediateCoefficients);
-    assert(d < k_maxNumberOfPolynomialCoefficients);
+    int d = childAtIndex(i)->getPolynomialCoefficients(symbolName, intermediateCoefficients);
+    assert(d < Expression::k_maxNumberOfPolynomialCoefficients);
     for (int j = 0; j < d+1; j++) {
-      static_cast<Addition *>(coefficients[j])->addOperand(intermediateCoefficients[j]);
+      static_cast<Addition *>(&coefficients[j])->addChildAtIndexInPlace(intermediateCoefficients[j], coefficients[j].numberOfChildren(), coefficients[j].numberOfChildren());
     }
   }
   return deg;
 }
 
-/* Layout */
+// Private
 
-bool Addition::needsParenthesesWithParent(const SerializationHelperInterface * e) const {
+// Layout
+bool AdditionNode::needsParenthesesWithParent(const SerializationHelperInterface * parentNode) const {
   Type types[] = {Type::Subtraction, Type::Opposite, Type::Multiplication, Type::Division, Type::Power, Type::Factorial};
-  return e->isOfType(types, 6);
+  return static_cast<const ExpressionNode *>(parentNode)->isOfType(types, 6);
 }
 
-/* Simplication */
+// Simplication
 
-Expression Addition::shallowReduce(Context& context, Preferences::AngleUnit angleUnit) {
+Expression AdditionNode::shallowReduce(Context& context, Preferences::AngleUnit angleUnit) const {
+  return Addition(this).shallowReduce(context, angleUnit);
+}
+
+Expression AdditionNode::shallowBeautify(Context & context, Preferences::AngleUnit angleUnit) const {
+  return Addition(this).shallowBeautify(context, angleUnit);
+}
+
+Expression AdditionNode::factorizeOnCommonDenominator(Context & context, Preferences::AngleUnit angleUnit) const {
+  //TODO
+  return Addition();
+#if 0
+  // We want to turn (a/b+c/d+e/b) into (a*d+b*c+e*d)/(b*d)
+
+  // Step 1: We want to compute the common denominator, b*d
+  Multiplication * commonDenominator = new Multiplication();
+  for (int i = 0; i < numberOfChildren(); i++) {
+    Expression * denominator = childAtIndex(i)->denominator(context, angleUnit);
+    if (denominator) {
+      // Make commonDenominator = LeastCommonMultiple(commonDenominator, denominator);
+      commonDenominator->addMissingFactors(denominator, context, angleUnit);
+      delete denominator;
+    }
+  }
+  if (commonDenominator->numberOfChildren() == 0) {
+    delete commonDenominator;
+    // If commonDenominator is empty this means that no child was a fraction.
+    return this;
+  }
+
+  // Step 2: Create the numerator. We start with this being a/b+c/d+e/b and we
+  // want to create numerator = a/b*b*d + c/d*b*d + e/b*b*d
+  AdditionNode * numerator = new AdditionNode();
+  for (int i=0; i < numberOfChildren(); i++) {
+    Multiplication * m = new Multiplication(childAtIndex(i), commonDenominator, true);
+    numerator->addOperand(m);
+    m->privateShallowReduce(context, angleUnit, true, false);
+  }
+  // Step 3: Add the denominator
+  Power * inverseDenominator = new Power(commonDenominator, new Rational(-1), false);
+  Multiplication * result = new Multiplication(numerator, inverseDenominator, false);
+
+  // Step 4: Simplify the numerator to a*d + c*b + e*d
+  numerator->shallowReduce(context, angleUnit);
+
+  // Step 5: Simplify the denominator (in case it's a rational number)
+  commonDenominator->deepReduce(context, angleUnit);
+  inverseDenominator->shallowReduce(context, angleUnit);
+
+  /* Step 6: We simplify the resulting multiplication forbidding any
+   * distribution of multiplication on additions (to avoid an infinite loop). */
+  return static_cast<Multiplication *>(replaceWith(result, true))->privateShallowReduce(context, angleUnit, false, true);
+#endif
+}
+
+void AdditionNode::factorizeOperands(Expression e1, Expression e2, Context & context, Preferences::AngleUnit angleUnit) {
+  // TODO
+#if 0
+  /* This function factorizes two children which only differ by a rational
+   * factor. For example, if this is AdditionNode(2*pi, 3*pi), then 2*pi and 3*pi
+   * could be merged, and this turned into AdditionNode(5*pi). */
+  assert(e1->parent() == this && e2->parent() == this);
+
+  // Step 1: Find the new rational factor
+  Rational * r = new Rational(Rational::AdditionNode(RationalFactor(e1), RationalFactor(e2)));
+
+  // Step 2: Get rid of one of the children
+  removeOperand(e2, true);
+
+  // Step 3: Use the new rational factor. Create a multiplication if needed
+  Multiplication * m = nullptr;
+  if (e1->type() == Type::Multiplication) {
+    m = static_cast<Multiplication *>(e1);
+  } else {
+    m = new Multiplication();
+    e1->replaceWith(m, false);
+    m->addOperand(e1);
+  }
+  if (m->childAtIndex(0)->type() == Type::Rational) {
+    m->replaceOperand(m->childAtIndex(0), r, true);
+  } else {
+    m->addOperand(r);
+  }
+
+  // Step 4: Reduce the multiplication (in case the new rational factor is zero)
+  m->shallowReduce(context, angleUnit);
+#endif
+}
+
+const Rational AdditionNode::RationalFactor(Expression e) {
+  if (e.type() == Type::Multiplication && e.childAtIndex(0).type() == Type::Rational) {
+    Rational result = Rational(static_cast<RationalNode *>(e.childAtIndex(0).node()));
+    return result;
+  }
+  return Rational(1);
+}
+
+static inline int NumberOfNonRationalFactors(const Expression e) {
+  // TODO should return a copy?
+  if (e.type() != ExpressionNode::Type::Multiplication) {
+    return 1; // Or (e->type() != Type::Rational);
+  }
+  int result = e.numberOfChildren();
+  if (e.childAtIndex(0).type() == ExpressionNode::Type::Rational) {
+    result--;
+  }
+  return result;
+}
+
+static inline const Expression FirstNonRationalFactor(const Expression e) {
+  // TODO should return a copy?
+  if (e.type() != ExpressionNode::Type::Multiplication) {
+    return e;
+  }
+  if (e.childAtIndex(0).type() == ExpressionNode::Type::Rational) {
+    return e.numberOfChildren() > 1 ? e.childAtIndex(1) : Expression();
+  }
+  return e.childAtIndex(0);
+}
+
+bool AdditionNode::TermsHaveIdenticalNonRationalFactors(const Expression e1, const Expression e2) {
+  return false;
+//TODO
+#if 0
+  /* Given two expressions, say wether they only differ by a rational factor.
+   * For example, 2*pi and pi do, 2*pi and 2*ln(2) don't. */
+
+  int numberOfNonRationalFactorsInE1 = NumberOfNonRationalFactors(e1);
+  int numberOfNonRationalFactorsInE2 = NumberOfNonRationalFactors(e2);
+
+  if (numberOfNonRationalFactorsInE1 != numberOfNonRationalFactorsInE2) {
+    return false;
+  }
+
+  int numberOfNonRationalFactors = numberOfNonRationalFactorsInE1;
+  if (numberOfNonRationalFactors == 1) {
+    return FirstNonRationalFactor(e1).isIdenticalTo(FirstNonRationalFactor(e2));
+  } else {
+    assert(numberOfNonRationalFactors > 1);
+    return Multiplication::HaveSameNonRationalFactors(e1, e2);
+  }
+#endif
+}
+
+Expression Addition::shallowBeautify(Context & context, Preferences::AngleUnit angleUnit) const {
+  Expression result = *this;
+  return result;
+  //TODO
+#if 0
+  /* Beautifying AdditionNode essentially consists in adding Subtractions if needed.
+   * In practice, we want to turn "a+(-1)*b" into "a-b". Or, more precisely, any
+   * "a+(-r)*b" into "a-r*b" where r is a positive Rational.
+   * Note: the process will slightly differ if the negative product occurs on
+   * the first term: we want to turn "AdditionNode(Multiplication(-1,b))" into
+   * "Opposite(b)".
+   * Last but not least, special care must be taken when iterating over children
+   * since we may remove some during the process. */
+
+  for (int i=0; i<numberOfChildren(); i++) {
+    if (childAtIndex(i)->type() != Type::Multiplication || childAtIndex(i)->childAtIndex(0)->type() != Type::Rational || childAtIndex(i)->childAtIndex(0)->sign() != Sign::Negative) {
+      // Ignore terms which are not like "(-r)*a"
+      continue;
+    }
+
+    Multiplication * m = static_cast<Multiplication *>(editableOperand(i));
+
+    if (static_cast<const Rational *>(m->childAtIndex(0))->isMinusOne()) {
+      m->removeOperand(m->childAtIndex(0), true);
+    } else {
+      m->editableOperand(0)->setSign(Sign::Positive, context, angleUnit);
+    }
+    Expression * subtractant = m->squashUnaryHierarchy();
+
+    if (i == 0) {
+      Opposite * o = new Opposite(subtractant, true);
+      replaceOperand(subtractant, o, true);
+    } else {
+      const Expression * op1 = childAtIndex(i-1);
+      removeOperand(op1, false);
+      Subtraction * s = new Subtraction(op1, subtractant, false);
+      replaceOperand(subtractant, s, false);
+      /* CAUTION: we removed a child. So we need to decrement i to make sure
+       * the next iteration is actually on the next child. */
+      i--;
+    }
+  }
+
+  return squashUnaryHierarchy();
+#endif
+}
+
+Expression Addition::shallowReduce(Context& context, Preferences::AngleUnit angleUnit) const {
+  // TODO
+  return Expression::shallowReduce(context, angleUnit);
+#if 0
   Expression * e = Expression::shallowReduce(context, angleUnit);
   if (e != this) {
     return e;
   }
-  /* Step 1: Addition is associative, so let's start by merging children which
+  /* Step 1: AdditionNode is associative, so let's start by merging children which
    * also are additions themselves. */
   int i = 0;
   int initialNumberOfOperands = numberOfChildren();
   while (i < initialNumberOfOperands) {
     Expression * o = editableOperand(i);
-    if (o->type() == Type::Addition) {
-      mergeOperands(static_cast<Addition *>(o));
+    if (o->type() == Type::AdditionNode) {
+      mergeOperands(static_cast<AdditionNode *>(o));
       continue;
     }
     i++;
   }
 
-  // Step 2: Sort the operands
+  // Step 2: Sort the children
   sortOperands(Expression::SimplificationOrder, true);
 
 #if MATRIX_EXACT_REDUCING
   /* Step 2bis: get rid of matrix */
   int n = 1;
   int m = 1;
-  /* All operands have been simplified so if any operand contains a matrix, it
-   * is at the root node of the operand. Moreover, thanks to the simplification
-   * order, all matrix operands (if any) are the last operands. */
+  /* All children have been simplified so if any child contains a matrix, it
+   * is at the root node of the child. Moreover, thanks to the simplification
+   * order, all matrix children (if any) are the last children. */
   Expression * lastOperand = editableOperand(numberOfChildren()-1);
   if (lastOperand->type() == Type::Matrix) {
-    // Create in-place the matrix of addition M (in place of the last operand)
+    // Create in-place the matrix of addition M (in place of the last child)
     Matrix * resultMatrix = static_cast<Matrix *>(lastOperand);
     n = resultMatrix->numberOfRows();
     m = resultMatrix->numberOfColumns();
     removeOperand(resultMatrix, false);
-    /* Scan (starting at the end) accross the addition operands to find any
+    /* Scan (starting at the end) accross the addition children to find any
      * other matrix */
     int i = numberOfChildren()-1;
-    while (i >= 0 && operand(i)->type() == Type::Matrix) {
+    while (i >= 0 && childAtIndex(i)->type() == Type::Matrix) {
       Matrix * currentMatrix = static_cast<Matrix *>(editableOperand(i));
       int on = currentMatrix->numberOfRows();
       int om = currentMatrix->numberOfColumns();
       if (on != n || om != m) {
         return replaceWith(new Undefined(), true);
       }
-      // Dispatch the current matrix operands in the created additions matrix
+      // Dispatch the current matrix children in the created additions matrix
       for (int j = 0; j < n*m; j++) {
-        Addition * a = new Addition();
+        AdditionNode * a = new AdditionNode();
         Expression * resultMatrixEntryJ = resultMatrix->editableOperand(j);
         resultMatrix->replaceOperand(resultMatrixEntryJ, a, false);
         a->addOperand(currentMatrix->editableOperand(j));
@@ -121,9 +307,9 @@ Expression Addition::shallowReduce(Context& context, Preferences::AngleUnit angl
       removeOperand(currentMatrix, true);
       i--;
     }
-    // Distribute the remaining addition on matrix operands
+    // Distribute the remaining addition on matrix children
     for (int i = 0; i < n*m; i++) {
-      Addition * a = static_cast<Addition *>(clone());
+      AdditionNode * a = static_cast<AdditionNode *>(clone());
       Expression * entryI = resultMatrix->editableOperand(i);
       resultMatrix->replaceOperand(entryI, a, false);
       a->addOperand(entryI);
@@ -142,7 +328,7 @@ Expression Addition::shallowReduce(Context& context, Preferences::AngleUnit angl
     if (o1->type() == Type::Rational && o2->type() == Type::Rational) {
       Rational r1 = *static_cast<Rational *>(o1);
       Rational r2 = *static_cast<Rational *>(o2);
-      Rational a = Rational::Addition(r1, r2);
+      Rational a = Rational::AdditionNode(r1, r2);
       replaceOperand(o1, new Rational(a), true);
       removeOperand(o2, true);
       continue;
@@ -156,7 +342,7 @@ Expression Addition::shallowReduce(Context& context, Preferences::AngleUnit angl
 
   /* Step 4: Let's remove zeroes if there's any. It's important to do this after
    * having factorized because factorization can lead to new zeroes. For example
-   * pi+(-1)*pi. We don't remove the last zero if it's the only operand left
+   * pi+(-1)*pi. We don't remove the last zero if it's the only child left
    * though. */
   i = 0;
   while (i < numberOfChildren()) {
@@ -168,196 +354,26 @@ Expression Addition::shallowReduce(Context& context, Preferences::AngleUnit angl
     i++;
   }
 
-  // Step 5: Let's remove the addition altogether if it has a single operand
+  // Step 5: Let's remove the addition altogether if it has a single child
   Expression * result = squashUnaryHierarchy();
 
   // Step 6: Last but not least, let's put everything under a common denominator
-  if (result == this && parent()->type() != Type::Addition) {
+  if (result == this && parent()->type() != Type::AdditionNode) {
     // squashUnaryHierarchy didn't do anything: we're not an unary hierarchy
     result = factorizeOnCommonDenominator(context, angleUnit);
   }
 
   return result;
+#endif
 }
 
-Expression * Addition::factorizeOnCommonDenominator(Context & context, Preferences::AngleUnit angleUnit) {
-  // We want to turn (a/b+c/d+e/b) into (a*d+b*c+e*d)/(b*d)
+template Complex<float> Poincare::AdditionNode::compute<float>(std::complex<float>, std::complex<float>);
+template Complex<double> Poincare::AdditionNode::compute<double>(std::complex<double>, std::complex<double>);
 
-  // Step 1: We want to compute the common denominator, b*d
-  Multiplication * commonDenominator = new Multiplication();
-  for (int i = 0; i < numberOfChildren(); i++) {
-    Expression * denominator = operand(i)->denominator(context, angleUnit);
-    if (denominator) {
-      // Make commonDenominator = LeastCommonMultiple(commonDenominator, denominator);
-      commonDenominator->addMissingFactors(denominator, context, angleUnit);
-      delete denominator;
-    }
-  }
-  if (commonDenominator->numberOfChildren() == 0) {
-    delete commonDenominator;
-    // If commonDenominator is empty this means that no operand was a fraction.
-    return this;
-  }
+template MatrixComplex<float> AdditionNode::computeOnMatrices<float>(const MatrixComplex<float>,const MatrixComplex<float>);
+template MatrixComplex<double> AdditionNode::computeOnMatrices<double>(const MatrixComplex<double>,const MatrixComplex<double>);
 
-  // Step 2: Create the numerator. We start with this being a/b+c/d+e/b and we
-  // want to create numerator = a/b*b*d + c/d*b*d + e/b*b*d
-  Addition * numerator = new Addition();
-  for (int i=0; i < numberOfChildren(); i++) {
-    Multiplication * m = new Multiplication(operand(i), commonDenominator, true);
-    numerator->addOperand(m);
-    m->privateShallowReduce(context, angleUnit, true, false);
-  }
-  // Step 3: Add the denominator
-  Power * inverseDenominator = new Power(commonDenominator, new Rational(-1), false);
-  Multiplication * result = new Multiplication(numerator, inverseDenominator, false);
-
-  // Step 4: Simplify the numerator to a*d + c*b + e*d
-  numerator->shallowReduce(context, angleUnit);
-
-  // Step 5: Simplify the denominator (in case it's a rational number)
-  commonDenominator->deepReduce(context, angleUnit);
-  inverseDenominator->shallowReduce(context, angleUnit);
-
-  /* Step 6: We simplify the resulting multiplication forbidding any
-   * distribution of multiplication on additions (to avoid an infinite loop). */
-  return static_cast<Multiplication *>(replaceWith(result, true))->privateShallowReduce(context, angleUnit, false, true);
-}
-
-void Addition::factorizeOperands(Expression * e1, Expression * e2, Context & context, Preferences::AngleUnit angleUnit) {
-  /* This function factorizes two operands which only differ by a rational
-   * factor. For example, if this is Addition(2*pi, 3*pi), then 2*pi and 3*pi
-   * could be merged, and this turned into Addition(5*pi). */
-  assert(e1->parent() == this && e2->parent() == this);
-
-  // Step 1: Find the new rational factor
-  Rational * r = new Rational(Rational::Addition(RationalFactor(e1), RationalFactor(e2)));
-
-  // Step 2: Get rid of one of the operands
-  removeOperand(e2, true);
-
-  // Step 3: Use the new rational factor. Create a multiplication if needed
-  Multiplication * m = nullptr;
-  if (e1->type() == Type::Multiplication) {
-    m = static_cast<Multiplication *>(e1);
-  } else {
-    m = new Multiplication();
-    e1->replaceWith(m, false);
-    m->addOperand(e1);
-  }
-  if (m->operand(0)->type() == Type::Rational) {
-    m->replaceOperand(m->operand(0), r, true);
-  } else {
-    m->addOperand(r);
-  }
-
-  // Step 4: Reduce the multiplication (in case the new rational factor is zero)
-  m->shallowReduce(context, angleUnit);
-}
-
-const Rational Addition::RationalFactor(Expression * e) {
-  if (e->type() == Type::Multiplication && e->operand(0)->type() == Type::Rational) {
-    return *(static_cast<const Rational *>(e->operand(0)));
-  }
-  return Rational(1);
-}
-
-static inline int NumberOfNonRationalFactors(const Expression * e) {
-  if (e->type() != Expression::Type::Multiplication) {
-    return 1; // Or (e->type() != Type::Rational);
-  }
-  int result = e->numberOfChildren();
-  if (e->operand(0)->type() == Expression::Type::Rational) {
-    result--;
-  }
-  return result;
-}
-
-static inline const Expression * FirstNonRationalFactor(const Expression * e) {
-  if (e->type() != Expression::Type::Multiplication) {
-    return e;
-  }
-  if (e->operand(0)->type() == Expression::Type::Rational) {
-    return e->numberOfChildren() > 1 ? e->operand(1) : nullptr;
-  }
-  return e->operand(0);
-}
-
-bool Addition::TermsHaveIdenticalNonRationalFactors(const Expression * e1, const Expression * e2) {
-  /* Given two expressions, say wether they only differ by a rational factor.
-   * For example, 2*pi and pi do, 2*pi and 2*ln(2) don't. */
-
-  int numberOfNonRationalFactorsInE1 = NumberOfNonRationalFactors(e1);
-  int numberOfNonRationalFactorsInE2 = NumberOfNonRationalFactors(e2);
-
-  if (numberOfNonRationalFactorsInE1 != numberOfNonRationalFactorsInE2) {
-    return false;
-  }
-
-  int numberOfNonRationalFactors = numberOfNonRationalFactorsInE1;
-  if (numberOfNonRationalFactors == 1) {
-    return FirstNonRationalFactor(e1)->isIdenticalTo(FirstNonRationalFactor(e2));
-  } else {
-    assert(numberOfNonRationalFactors > 1);
-    return Multiplication::HaveSameNonRationalFactors(e1, e2);
-  }
-}
-
-Expression * Addition::shallowBeautify(Context & context, Preferences::AngleUnit angleUnit) {
-  /* Beautifying Addition essentially consists in adding Subtractions if needed.
-   * In practice, we want to turn "a+(-1)*b" into "a-b". Or, more precisely, any
-   * "a+(-r)*b" into "a-r*b" where r is a positive Rational.
-   * Note: the process will slightly differ if the negative product occurs on
-   * the first term: we want to turn "Addition(Multiplication(-1,b))" into
-   * "Opposite(b)".
-   * Last but not least, special care must be taken when iterating over operands
-   * since we may remove some during the process. */
-
-  for (int i=0; i<numberOfChildren(); i++) {
-    if (operand(i)->type() != Type::Multiplication || operand(i)->operand(0)->type() != Type::Rational || operand(i)->operand(0)->sign() != Sign::Negative) {
-      // Ignore terms which are not like "(-r)*a"
-      continue;
-    }
-
-    Multiplication * m = static_cast<Multiplication *>(editableOperand(i));
-
-    if (static_cast<const Rational *>(m->operand(0))->isMinusOne()) {
-      m->removeOperand(m->operand(0), true);
-    } else {
-      m->editableOperand(0)->setSign(Sign::Positive, context, angleUnit);
-    }
-    Expression * subtractant = m->squashUnaryHierarchy();
-
-    if (i == 0) {
-      Opposite * o = new Opposite(subtractant, true);
-      replaceOperand(subtractant, o, true);
-    } else {
-      const Expression * op1 = operand(i-1);
-      removeOperand(op1, false);
-      Subtraction * s = new Subtraction(op1, subtractant, false);
-      replaceOperand(subtractant, s, false);
-      /* CAUTION: we removed an operand. So we need to decrement i to make sure
-       * the next iteration is actually on the next operand. */
-      i--;
-    }
-  }
-
-  return squashUnaryHierarchy();
-}
-
-/* Evaluation */
-
-template<typename T>
-std::complex<T> Addition::compute(const std::complex<T> c, const std::complex<T> d) {
-  return c+d;
-}
-
-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 MatrixComplex<float> Addition::computeOnMatrices<float>(const MatrixComplex<float>,const MatrixComplex<float>);
-template MatrixComplex<double> Addition::computeOnMatrices<double>(const MatrixComplex<double>,const MatrixComplex<double>);
-
-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>);
+template MatrixComplex<float> AdditionNode::computeOnComplexAndMatrix<float>(std::complex<float> const, const MatrixComplex<float>);
+template MatrixComplex<double> AdditionNode::computeOnComplexAndMatrix<double>(std::complex<double> const, const MatrixComplex<double>);
 
 }