diff --git a/include/expolit/differentiation.hh b/include/expolit/differentiation.hh
index 8faed4a..ac3b0cf 100644
--- a/include/expolit/differentiation.hh
+++ b/include/expolit/differentiation.hh
@@ -1,143 +1,143 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2018 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Expolit is free software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Expolit is  distributed in the  hope that it  will be useful, but  WITHOUT
  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  * FOR A PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License for
  * more details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Expolit. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 #ifndef DIFFERENTIATION_HH
 #define DIFFERENTIATION_HH
 
 #include "exponential.hh"
 #include "litteral.hh"
 #include "operations.hh"
 #include "polynomial.hh"
 #include "types.hh"
 #include <algorithm>
 
 namespace expolit {
 
 namespace detail {
 template <std::size_t i>
 struct for_diff {
   template <typename CT, std::size_t n>
   constexpr static void diff(std::array<CT, n - 1>& res,
                              const std::array<CT, n>& a) {
-    std::get<i>(res) = a[i + 1] * static_cast<Real>(i + 1);
+    std::get<i>(res) = a[i + 1] * static_cast<CT>(i + 1);
     for_diff<i - 1>::diff(res, a);
   }
 };
 
 template <>
 struct for_diff<0> {
   template <typename CT, std::size_t n>
   constexpr static void diff(std::array<CT, n - 1>& res,
                              const std::array<CT, n>& a) {
-    std::get<0>(res) = a[1] * static_cast<Real>(1);
+    std::get<0>(res) = a[1] * static_cast<CT>(1);
   }
 };
 }  // namespace detail
 
 /// Differentiate polynomial
 template <typename CT, UInt order>
 constexpr auto differentiate(const Polynomial<CT, order>& p) {
   Polynomial<CT, order - 1> dp;
   detail::for_diff<dp.terms - 1>::diff(dp.coeffs, p.coeffs);
   return dp;
 }
 
 /// Differentiate constant
 template <typename CT>
 constexpr auto differentiate(const Constant<CT>& /*p*/) {
   return Constant<CT>({0});
 }
 
 /// Differentiate litteral
 template <typename Tag>
 constexpr auto differentiate(const Litteral<Tag>& /*p*/) {
   return 0;
 }
 
 /// Differentiate exponential
 template <typename T>
 constexpr auto differentiate(const Exponential<T>& e) {
   return differentiate(e.expression) * e;
 }
 
 /// Differentiate generic product
 template <typename T, typename U>
 constexpr auto differentiate(const Product<T, U>& p) {
   auto&& o = p.operands;
   return differentiate(o.first) * o.second + o.first * differentiate(o.second);
 }
 
 /// Differentiate constant + expression
 template <typename CT, typename U>
 constexpr auto differentiate(const Sum<Constant<CT>, U>& p) {
   return differentiate(p.operands.second);
 }
 
 /// Differentiate expression + constant
 template <typename T, typename U, typename CT>
 constexpr auto differentiate(const Sum<U, Constant<CT>>& p) {
   return differentiate(p.commute());
 }
 
 /// Differentiate litteral + expression
 template <typename Tag, typename U>
 constexpr auto differentiate(const Sum<Litteral<Tag>, U>& p) {
   return differentiate(p.operands.second);
 }
 
 /// Differentiate expression + litteral
 template <typename Tag, typename U>
 constexpr auto differentiate(const Sum<U, Litteral<Tag>>& p) {
   return differentiate(p.commute());
 }
 
 /// Differentiate constant * expression
 template <typename CT, typename U>
 constexpr auto differentiate(const Product<Constant<CT>, U>& p) {
   return p.operands.first * differentiate(p.operands.second);
 }
 
 /// Differentiate expression * constant
 template <typename CT, typename U>
 constexpr auto differentiate(const Product<U, Constant<CT>>& p) {
   return differentiate(p.commute());
 }
 
 /// Differentiate litteral * expression
 template <typename Tag, typename U>
 constexpr auto differentiate(const Product<Litteral<Tag>, U>& p) {
   return p.operands.first * differentiate(p.operands.second);
 }
 
 /// Differentiate expression * litteral
 template <typename Tag, typename U>
 constexpr auto differentiate(const Product<U, Litteral<Tag>>& p) {
   return differentiate(p.commute());
 }
 
 }  // namespace expolit
 
 #endif
diff --git a/include/expolit/integration.hh b/include/expolit/integration.hh
index b002233..77ecb12 100644
--- a/include/expolit/integration.hh
+++ b/include/expolit/integration.hh
@@ -1,155 +1,155 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2018 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Expolit is free software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Expolit is  distributed in the  hope that it  will be useful, but  WITHOUT
  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  * FOR A PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License for
  * more details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Expolit. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 #ifndef INTEGRATION_HH
 #define INTEGRATION_HH
 
 #include "differentiation.hh"
 #include "exponential.hh"
 #include "litteral.hh"
 #include "polynomial.hh"
 #include "types.hh"
 
 namespace expolit {
 
 namespace detail {
 template <std::size_t i>
 struct for_integ {
   template <typename CT, std::size_t n>
   static constexpr void integ(std::array<CT, n+1>& res,
 			      const std::array<CT, n>& a) {
-    std::get<i>(res) = a[i-1] / static_cast<Real>(i);
+    std::get<i>(res) = a[i-1] / static_cast<CT>(i);
     for_integ<i-1>::integ(res, a);
   }
 };
 
 template <>
 struct for_integ<1> {
   template <typename CT, std::size_t n>
   static constexpr void integ(std::array<CT, n+1>& res,
 			      const std::array<CT, n>& a) {
     std::get<1>(res) = a[0];
   }
 };
 }  // namespace detail
 
 /// Integrate polynomial
 template <typename CT, UInt order>
 constexpr auto integrate(const Polynomial<CT, order>& p) {
   Polynomial<CT, order + 1> ip;
   detail::for_integ<ip.terms-1>::integ(ip.coeffs, p.coeffs);
   return ip;
 }
 
 /// Integrate litteral
 template <typename Tag>
 constexpr auto integrate(const Litteral<Tag>& q) {
   constexpr Polynomial<Int, 1> id({0, 1});
   return q * id;
 }
 
 /// Integrate constant * expression
 template <typename CT, typename Derived>
 constexpr auto integrate(const Product<Constant<CT>, Derived>& p) {
   return p.operands.first * integrate(p.operands.second);
 }
 
 /// Integrate expression * constant
 template <typename CT, typename Derived>
 constexpr auto integrate(const Product<Derived, Constant<CT>>& p) {
   return integrate(p.commute());
 }
 
 /// Integrate litteral * expression
 template <typename Tag, typename Derived>
 constexpr auto integrate(const Product<Litteral<Tag>, Derived>& p) {
   return p.operands.first * integrate(p.operands.second);
 }
 
 /// Integrate expression * litteral
 template <typename Tag, typename Derived>
 constexpr auto integrate(const Product<Derived, Litteral<Tag>>& p) {
   return integrate(p.commute());
 }
 
 /// Integrate litteral / expression
 template <typename Tag, typename Derived>
 constexpr auto integrate(const Division<Derived, Litteral<Tag>>& p) {
   return integrate(p.operands.first) / p.operands.second;
 }
 
 
 /// Integrate exponential of order 1 polynomial
 template <typename CT>
 constexpr auto integrate(const Exponential<Polynomial<CT, 1>>& e) {
   return Constant<CT>({1 / e.expression.coeffs.back()}) * e;
 }
 
 /// Integrate exponential of litteral * order 1 polynomial
 template <typename Tag, typename CT>
 constexpr auto
 integrate(const Exponential<Product<Litteral<Tag>, Polynomial<CT, 1>>>& e) {
   return e * Constant<CT>({1 / e.expression.operands.second.coeffs.back()}) /
          e.expression.operands.first;
 }
 
 /// Integrate exponential of order 1 polynomial * litteral
 template <typename Tag, typename CT>
 constexpr auto
 integrate(const Exponential<Product<Polynomial<CT, 1>, Litteral<Tag>>>& e) {
   return integrate(e.op(e.expression.commute()));
 }
 
 /// Integration by parts of expression * polynomial
 template <typename Derived, typename CT, UInt order>
 constexpr auto integrate(const Product<Derived, Polynomial<CT, order>>& p) {
   auto u = integrate(p.operands.first);
   return u * p.operands.second -
          integrate(u * differentiate(p.operands.second));
 }
 
 /// Integration by parts of polynomial * expression
 template <typename Derived, typename CT, UInt order>
 constexpr auto integrate(const Product<Polynomial<CT, order>, Derived>& p) {
   return integrate(p.commute());
 }
 
 /// Integration of sum
 template <typename Der1, typename Der2>
 constexpr auto integrate(const Sum<Der1, Der2>& s) {
   return integrate(s.operands.first) + integrate(s.operands.second);
 }
 
 /// Compute definite integral
 template <typename T, typename Derived>
 constexpr auto definite_integral(const std::pair<T, T>& bounds,
                                  const Expression<Derived>& exp) {
   auto integ = integrate(exp.downcast());
   return integ(bounds.second) - integ(bounds.first);
 }
 
 }  // namespace expolit
 #endif
diff --git a/include/expolit/types.hh b/include/expolit/types.hh
index fd25d69..d295c7e 100644
--- a/include/expolit/types.hh
+++ b/include/expolit/types.hh
@@ -1,53 +1,52 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2018 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Expolit is free software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Expolit is  distributed in the  hope that it  will be useful, but  WITHOUT
  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  * FOR A PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License for
  * more details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Expolit. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 #ifndef TYPES_HH
 #define TYPES_HH
 
 #include <cstddef>
 #include <ostream>
 #include <type_traits>
 #include <utility>
 
 namespace expolit {
-using Real = double;
 using Int = int;
 using UInt = unsigned int;
 
 /// Expression class for expression templates
 template <typename Derived>
 struct Expression {
   constexpr const Derived& downcast() const { return static_cast<const Derived&>(*this); }
   constexpr Derived& downcast() { return static_cast<Derived&>(*this); }
 };
 
 template <typename Derived>
 std::ostream& operator<<(std::ostream& os, const Expression<Derived>& e) {
   os << e.downcast();
   return os;
 }
 }  // namespace expolit
 #endif
diff --git a/tests/tests_gtest.cc b/tests/tests_gtest.cc
index bd47dca..258f1b9 100644
--- a/tests/tests_gtest.cc
+++ b/tests/tests_gtest.cc
@@ -1,174 +1,175 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2018 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Expolit is free software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Expolit is  distributed in the  hope that it  will be useful, but  WITHOUT
  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  * FOR A PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License for
  * more details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Expolit. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 #include <expolit/expolit>
 #include <complex>
 #include <gtest/gtest.h>
 
 using namespace expolit;
+using Real = double;
 
 static constexpr Polynomial<Int, 1> Zi({0, 1});
 static constexpr Polynomial<std::complex<Real>, 1> Zc({0, 1});
 static constexpr Polynomial<Real, 1> Zr({0, 1});
 
 using lin = std::remove_cv_t<decltype(Zr)>;
 
 template <typename CT, UInt order>
 constexpr void compare(const Polynomial<CT, order>& p1,
                        const Polynomial<CT, order>& p2) {
   for (UInt i = 0; i < p1.terms; ++i)
     ASSERT_NEAR(std::abs(p1.coeffs[i] - p2.coeffs[i]), 0, 1e-14)
         << "Polynomials are not equal (i = " << i << ")";
 }
 
 template <typename T, UInt order>
 constexpr void compare(const Polynomial<Int, order>& p1,
                        const Polynomial<Int, order>& p2) {
   for (UInt i = 0; i < p1.terms; ++i)
     ASSERT_EQ(p1.coeffs[i], p2.coeffs[i])
         << "Polynomials are not equal (i = " << i << ")";
 }
 
 TEST(Polynomial, initialization) {
   constexpr auto p = 1 + 2 * pow<2>(Zi) + pow<3>(Zi);
   constexpr Int z = 3;
   constexpr Int sol = 1 + 2 * z * z + z * z * z;
   constexpr Int res = p(z);
   static_assert(sol == res, "Bad coefficients");
 
   // This guy cannot be consexpr because std::complex is stupid
   auto p2 = 1. + Zc * Zc;
   std::complex<Real> z2(0, 1);
   ASSERT_NEAR(std::abs(p2(z2)), 0, 1e-14) << "Bad coefficents (complex)";
 }
 
 TEST(Polynomial, sum) {
   constexpr auto p = 1 + 2 * pow<2>(Zi) + pow<3>(Zi);
   constexpr auto p2 = Zi + 2 * Zi * Zi;
   constexpr auto res = 1 + Zi + 4 * Zi * Zi + Zi * Zi * Zi;
   compare(p + p2, res);
 }
 
 TEST(Polynomial, product) {
   constexpr auto p = 1 + 2 * pow<2>(Zi) + pow<3>(Zi);
   constexpr auto p2 = Zi + 2 * pow<2>(Zi);
   constexpr auto res =
       Zi + 2 * pow<2>(Zi) + 2 * pow<3>(Zi) + 5 * pow<4>(Zi) + 2 * pow<5>(Zi);
   compare(p * p2, res);
 }
 
 TEST(Polynomial, derivative) {
   constexpr auto p = 1 + 2 * pow<2>(Zi) + pow<3>(Zi);
   constexpr auto res = 4 * Zi + 3 * Zi * Zi;
   compare(differentiate(p), res);
 }
 
 TEST(Polynomial, integration) {
   constexpr auto p = 1. + 2. * pow<2>(Zr) + pow<3>(Zr);
   constexpr auto res = Zr + 2 / 3. * pow<3>(Zr) + 1 / 4. * pow<4>(Zr);
   compare(integrate(p), res);
 }
 
 TEST(Exponential, evaluation) {
   constexpr auto e = exp(2. * Zr);
   constexpr Real z = 3.4;
   const Real res = e(z);  // cannot be constexpr because of exponential
   ASSERT_NEAR(std::exp(2 * z), res, 1e-14) << "Exponential evaluation fail";
 }
 
 TEST(Exponential, derivative) {
   constexpr auto e = exp(2. * Zr);
   constexpr auto diff = differentiate(e);
   constexpr auto res = Constant<Real>({2}) * e;
   compare(res.operands.first, diff.operands.first);
   compare(e.expression, diff.operands.second.expression);
 }
 
 TEST(Exponential, integration) {
   constexpr auto e = exp(2. * Zr);
   constexpr auto integ = integrate(e);
   constexpr auto res = Constant<Real>({1 / 2.}) * e;
   compare(res.operands.first, integ.operands.first);
   compare(e.expression, integ.operands.second.expression);
 }
 
 TEST(Integration, parts) {
   constexpr auto e = exp(2. * Zr);
   constexpr auto expr = e * (2. * Zr);
   constexpr auto res = Zr * e - Constant<Real>({1 / 2.}) * e;
   constexpr auto integ = integrate(expr);
   compare(integ.operands.first.operands.first,
           res.operands.first.operands.first);
   compare(integ.operands.second.operands.first,
           res.operands.second.operands.first);
 }
 
 TEST(Integration, definite_integral) {
   constexpr auto p = Zr;
   constexpr std::pair<Real, Real> bounds(0, 1);
   constexpr auto res = definite_integral(bounds, p);
   ASSERT_NEAR(0.5, res, 1e-14) << "Definite integral fail";
 
   constexpr auto pi = 2 * Zi;
   static_assert(definite_integral(std::make_pair(0, 1), pi) == 1,
                 "Definite integral fail");
 }
 
 using Q = Litteral<struct Q_>;
 
 TEST(Litteral, substitution) {
   constexpr Q q;
   constexpr auto p = q * Zi;
   constexpr auto c = Constant<Int>({2});
   constexpr auto p_sub = substitute<Q::tag>(c, p);
   compare(2 * Zi, p_sub);
 
   constexpr auto expo = exp(q * Zr);
   constexpr auto x_sub = substitute<Q::tag>(c, expo);
   compare(2. * Zr, x_sub.expression);
 }
 
 TEST(Litteral, differentiation) {
   constexpr Q q;
   constexpr auto c = Constant<Real>({2});
   constexpr auto expo = exp(q * Zr);
   constexpr auto dexpo = differentiate(expo);
   constexpr auto res = substitute<Q::tag>(c, dexpo);
   compare(c, res.operands.first);
   compare(2. * Zr, res.operands.second.expression);
 }
 
 TEST(Litteral, integration) {
   constexpr Q q;
   constexpr auto c = Constant<Real>({2});
   constexpr auto res = substitute<Q::tag>(c, integrate(2 * Zr * exp(q * Zr)));
   constexpr auto sol =
       Zr * exp(2. * Zr) - Constant<Real>({1 / 2.}) * exp(2. * Zr);
   compare(sol.operands.first.operands.first,
           res.operands.first.operands.first);
   compare(sol.operands.second.operands.first,
           res.operands.second.operands.first);
 }