diff --git a/include/expolit/differentiation.hh b/include/expolit/differentiation.hh
index c00e2a0..8faed4a 100644
--- a/include/expolit/differentiation.hh
+++ b/include/expolit/differentiation.hh
@@ -1,123 +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);
+    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);
+  }
+};
+}  // namespace detail
+
 /// Differentiate polynomial
 template <typename CT, UInt order>
 constexpr auto differentiate(const Polynomial<CT, order>& p) {
   Polynomial<CT, order - 1> dp;
-  for (UInt i = 0; i < dp.terms; ++i)
-    dp.coeffs[i] = p.coeffs[i + 1] * static_cast<Real>(i + 1);
+  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 e818965..b002233 100644
--- a/include/expolit/integration.hh
+++ b/include/expolit/integration.hh
@@ -1,134 +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);
+    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;
-  for (UInt i = 1; i < ip.terms; ++i)
-    ip.coeffs[i] = p.coeffs[i - 1] / static_cast<Real>(i);
+  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/litteral.hh b/include/expolit/litteral.hh
index 40b40c6..7023274 100644
--- a/include/expolit/litteral.hh
+++ b/include/expolit/litteral.hh
@@ -1,79 +1,80 @@
 /**
  * @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 LITTERAL_HH
 #define LITTERAL_HH
 
 #include "operations.hh"
 #include "polynomial.hh"
 #include "types.hh"
 
 namespace expolit {
 template <typename Tag>
-struct Litteral : public Expression<Litteral<Tag>> {
+struct Litteral : Expression<Litteral<Tag>> {
   using tag = Tag;
+  Litteral() = default;
 };
 
 /// Substitue litteral
 template <typename Tag, typename New>
 constexpr auto substitute(const Expression<New>& val,
 			  const Litteral<Tag>& /*t*/) {
   return val.downcast();
 }
 
 /// Substitute litteral in binary operation
 template <typename Tag, typename New, typename Der1, typename Der2,
           typename Op>
 constexpr auto substitute(const Expression<New>& val,
                           const BinaryOperation<Der1, Der2, Op>& op) {
   return op.op(substitute<Tag>(val, op.operands.first),
                substitute<Tag>(val, op.operands.second));
 }
 
 /// Substitute litteral in polynomial (aka do nothing)
 template <typename Tag, typename New, typename CT, UInt order>
 constexpr auto substitute(const Expression<New>& /*val*/,
 			  const Polynomial<CT, order>& p) {
   return p;
 }
 
 /// Substitute litteral in polynomial (aka do nothing)
 template <typename Tag, typename New, typename T, typename Op>
 constexpr auto substitute(const Expression<New>& val,
 			  const UnaryOperation<T, Op>& p) {
   return p.op(substitute<Tag>(val, p.expression));
 }
 
 /// Substitute generic expression by downcast
 // template <typename Tag, typename New, typename Derived>
 // constexpr auto substitute(const Expression<New>& val,
 // 			  const Expression<Derived>& e) {
 //   return substitute<Tag>(val, e.downcast());
 // }
 
 }  // namespace expolit
 
 #endif
diff --git a/include/expolit/operations.hh b/include/expolit/operations.hh
index f3d3546..32407a9 100644
--- a/include/expolit/operations.hh
+++ b/include/expolit/operations.hh
@@ -1,180 +1,189 @@
 /**
  * @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 OPERATIONS_HH
 #define OPERATIONS_HH
 
 #include "types.hh"
 #include <functional>
 #include <ostream>
 #include <utility>
 
 namespace expolit {
 /**
  * @brief Type for unary operation
  */
 template <typename T, typename Op>
-struct UnaryOperation : public Expression<UnaryOperation<T, Op>> {
+struct UnaryOperation : Expression<UnaryOperation<T, Op>> {
   /// Construct from an expression
   constexpr UnaryOperation(const Expression<T>& t) : expression(t.downcast()) {}
   /// Evaluate operation
   template <typename Val>
   constexpr auto operator()(Val&& v) const {
     return op(expression(v));
   }
 
   const T expression;
   const Op op = Op();
 };
 
 /**
  * @brief Type for binary operation
  *
  * Handles binary operations between two arbitrary expressions
  */
 template <typename T, typename U, typename Op>
-struct BinaryOperation : public Expression<BinaryOperation<T, U, Op>> {
+struct BinaryOperation : Expression<BinaryOperation<T, U, Op>> {
   /// Construct from two expressions
   constexpr BinaryOperation(const Expression<T>& t, const Expression<U>& u)
       : operands(t.downcast(), u.downcast()) {}
 
   /// Evaluate operation
   template <typename Val>
   constexpr auto operator()(Val&& v) const {
     return op(operands.first(v), operands.second(v));
   }
 
   constexpr auto commute() const {
     return BinaryOperation<U, T, Op>(operands.second, operands.first);
   }
 
   /// Product operarands
   const std::pair<T, U> operands;
   const Op op = Op();
 };
 
 /// Symmetric of T, U
 template <typename T, typename U, typename Op>
 using SBinaryOperation = BinaryOperation<U, T, Op>;
 
 /// Product expression
 template <typename T, typename U>
 using Product = BinaryOperation<T, U, std::multiplies<>>;
 
 /// Sum expression
 template <typename T, typename U>
 using Sum = BinaryOperation<T, U, std::plus<>>;
 
 /// Difference expression
 template <typename T, typename U>
 using Difference = BinaryOperation<T, U, std::minus<>>;
 
 /// Division expression
 template <typename T, typename U>
 using Division = BinaryOperation<T, U, std::divides<>>;
 
 /* -------------------------------------------------------------------------- */
 /* Algebraic rules */
 /* -------------------------------------------------------------------------- */
 /// Product operator
 template <typename Der1, typename Der2>
 constexpr auto operator*(const Expression<Der1>& e1,
                          const Expression<Der2>& e2) {
   return Product<Der1, Der2>(e1, e2);
 }
 
 /// Sum operator
 template <typename Der1, typename Der2>
 constexpr auto operator+(const Expression<Der1>& e1,
                          const Expression<Der2>& e2) {
   return Sum<Der1, Der2>(e1, e2);
 }
 
 /// Difference operator
 template <typename Der1, typename Der2>
 constexpr auto operator-(const Expression<Der1>& e1,
                          const Expression<Der2>& e2) {
   return Difference<Der1, Der2>(e1, e2);
 }
 
 /// Division operator
 template <typename Der1, typename Der2>
 constexpr auto operator/(const Expression<Der1>& e1,
-			 const Expression<Der2>& e2) {
+                         const Expression<Der2>& e2) {
   return Division<Der1, Der2>(e1, e2);
 }
 
 namespace detail {
-template <UInt n, typename DerAcc, typename Derived>
-constexpr auto power_acc(const Expression<DerAcc>& acc,
-                         const Expression<Derived>& e) {
-  if constexpr (n > 1)
-    return power_acc<n - 1>(acc.downcast() * e.downcast(), e);
-  else
+template <UInt n>
+struct power_acc {
+  template <typename DerAcc, typename Derived>
+  static constexpr auto accumulate(const Expression<DerAcc>& acc,
+                                   const Expression<Derived>& e) {
+    return power_acc<n - 1>::accumulate(acc.downcast() * e.downcast(), e);
+  }
+};
+
+template <>
+struct power_acc<1> {
+  template <typename DerAcc, typename Derived>
+  static constexpr auto accumulate(const Expression<DerAcc>& acc,
+                                   const Expression<Derived>&) {
     return acc.downcast();
-}
+  }
+};
+
 }  // namespace detail
 
 template <UInt n, typename Derived>
 constexpr auto pow(const Expression<Derived>& e) {
-  return detail::power_acc<n>(e, e);
+  return detail::power_acc<n>::accumulate(e, e);
 }
 
-
 /* -------------------------------------------------------------------------- */
 /* Output */
 /* -------------------------------------------------------------------------- */
 
 template <typename T, typename U>
 std::ostream& operator<<(std::ostream& os, const Product<T, U>& p) {
   os << "(" << p.operands.first << ")"
      << " * "
      << "(" << p.operands.second << ")";
   return os;
 }
 
 template <typename T, typename U>
 std::ostream& operator<<(std::ostream& os, const Division<T, U>& p) {
   os << "(" << p.operands.first << ")"
      << " / "
      << "(" << p.operands.second << ")";
   return os;
 }
 
 template <typename T, typename U>
 std::ostream& operator<<(std::ostream& os, const Difference<T, U>& p) {
   os << p.operands.first << " - " << p.operands.second;
   return os;
 }
 
 template <typename T, typename U>
 std::ostream& operator<<(std::ostream& os, const Sum<T, U>& p) {
   os << p.operands.first << " + " << p.operands.second;
   return os;
 }
 
 }  // namespace expolit
 #endif
diff --git a/include/expolit/polynomial.hh b/include/expolit/polynomial.hh
index 960e65b..2743f73 100644
--- a/include/expolit/polynomial.hh
+++ b/include/expolit/polynomial.hh
@@ -1,232 +1,318 @@
 /**
  * @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 POLYNOMIALS_HH
 #define POLYNOMIALS_HH
 
 #include "operations.hh"
 #include "types.hh"
 #include <algorithm>
 #include <array>
 #include <numeric>
 #include <ostream>
 
 namespace expolit {
 
 template <typename From, typename To>
-using compatible = std::enable_if_t<std::is_convertible_v<From, To>>;
+using compatible = std::enable_if_t<std::is_convertible<From, To>::value>;
 
 /**
  * @brief Polynomial type
  *
  * Type for polynomial expressions of arbitrary order, variable type and
  * coefficient type
  */
 template <typename CoeffType, UInt order>
-struct Polynomial : public Expression<Polynomial<CoeffType, order>> {
+struct Polynomial : Expression<Polynomial<CoeffType, order>> {
   /// Number of terms
   constexpr static UInt terms = order + 1;
 
   constexpr Polynomial() = default;
   constexpr explicit Polynomial(const std::array<CoeffType, terms>& coeffs)
       : coeffs(coeffs) {}
 
   /// Evaluating polynomial
   template <typename Val, typename = compatible<Val, CoeffType>>
   constexpr auto operator()(Val z) const {
     using Res = decltype(std::declval<CoeffType>() * z);
     Res res(coeffs.front());
     Val power(z);
     for (UInt i = 1; i < terms; ++i) {
       res += coeffs[i] * power;
       power *= z;
     }
     return res;
   }
 
   /// Storing coefficients
   std::array<CoeffType, terms> coeffs = {0};
 };
 
 template <typename CoeffType>
 using Constant = Polynomial<CoeffType, 0>;
 
 /* -------------------------------------------------------------------------- */
 /* Algebraic rules */
 /* -------------------------------------------------------------------------- */
 
 /// Sum with scalar
 template <typename CT, UInt order, typename DT, typename = compatible<DT, CT>>
 constexpr auto operator+(DT x, const Polynomial<CT, order>& p) {
   Polynomial<decltype(std::declval<CT>() + x), order> res(p);
-  res.coeffs.front() += x;
+  std::get<0>(res.coeffs) += x;
   return res;
 }
 
 /// Sum with scalar
 template <typename CT, UInt order, typename DT, typename = compatible<DT, CT>>
 constexpr auto operator+(const Polynomial<CT, order>& p, DT x) {
   return x + p;
 }
 
+namespace detail {
+template <std::size_t i>
+struct for_scalar {
+  template <typename CT, std::size_t n, typename DT>
+  constexpr static void prod(std::array<CT, n>& res, DT x) {
+    std::get<i>(res) *= x;
+    for_scalar<i-1>::prod(res, x);
+  }
+
+  template <typename CT, std::size_t n, typename DT>
+  constexpr static void div(std::array<CT, n>& res, DT x) {
+    std::get<i>(res) /= x;
+    for_scalar<i-1>::prod(res, x);
+  }
+};
+
+template <>
+struct for_scalar<0> {
+  template <typename CT, std::size_t n, typename DT>
+  constexpr static void prod(std::array<CT, n>& res, DT x) {
+    std::get<0>(res) *= x;
+  }
+
+  template <typename CT, std::size_t n, typename DT>
+  constexpr static void div(std::array<CT, n>& res, DT x) {
+    std::get<0>(res) /= x;
+  }
+};
+}  // namespace detail
+
 /// Product with scalar
 template <typename CT, UInt order, typename DT, typename = compatible<DT, CT>>
 constexpr auto operator*(DT x, const Polynomial<CT, order>& p) {
   Polynomial<decltype(std::declval<CT>() * x), order> res(p);
-  for (auto&& v : res.coeffs)
-    v *= x;
+  detail::for_scalar<res.terms-1>::prod(res.coeffs, x);
   return res;
 }
 
 /// Product with scalar
 template <typename CT, UInt order, typename DT, typename = compatible<DT, CT>>
 constexpr auto operator*(const Polynomial<CT, order>& p, DT x) {
   return x * p;
 }
 
 /// Division with scalar
 template <typename CT, UInt order, typename DT, typename = compatible<DT, CT>>
 constexpr auto operator/(const Polynomial<CT, order>& p, DT x) {
   Polynomial<decltype(std::declval<CT>() / x), order> res(p);
-  for (auto&& v : res.coeffs)
-    v /= x;
+  detail::for_scalar<res.terms-1>::div(res.coeffs, x);
   return res;
 }
 
 /// Division with constant
 template <typename CT, UInt order, typename CT2, typename = compatible<CT2, CT>>
 constexpr auto operator/(const Polynomial<CT, order>& p,
                          const Constant<CT2>& x) {
   Polynomial<decltype(std::declval<CT>() / x.coeffs.front()), order> res(p);
-  for (auto&& v : res.coeffs)
-    v /= x.coeffs.front();
+  detail::for_scalar<res.terms-1>::div(res.coeffs, std::get<0>(x.coeffs));
   return res;
 }
 
+namespace detail {
+template <std::size_t i, std::size_t j>
+struct coeff_product {
+  template <typename CT1, typename CT2, typename CT3, std::size_t n,
+            std::size_t m, std::size_t l>
+  static constexpr void prod(std::array<CT1, n>& res,
+                             const std::array<CT2, m>& a,
+                             const std::array<CT3, l>& b) {
+    std::get<i + j>(res) += a[i] * b[j];
+    coeff_product<i, j - 1>::prod(res, a, b);
+  }
+};
+
+template <std::size_t i>
+struct coeff_product<i, 0> {
+  template <typename CT1, typename CT2, typename CT3, std::size_t n,
+            std::size_t m, std::size_t l>
+  static constexpr void prod(std::array<CT1, n>& res,
+                             const std::array<CT2, m>& a,
+                             const std::array<CT3, l>& b) {
+    std::get<i>(res) += a[i] * b[0];
+    coeff_product<i - 1, l - 1>::prod(res, a, b);
+  }
+};
+
+template <>
+struct coeff_product<0, 0> {
+  template <typename CT1, typename CT2, typename CT3, std::size_t n,
+            std::size_t m, std::size_t l>
+  static constexpr void prod(std::array<CT1, n>& res,
+                             const std::array<CT2, m>& a,
+                             const std::array<CT3, l>& b) {
+    std::get<0>(res) += a[0] * b[0];
+  }
+};
+}  // namespace detail
+
 /// Product of polynomials
 template <typename CT1, typename CT2, UInt o1, UInt o2>
 constexpr auto operator*(const Polynomial<CT1, o1>& p1,
                          const Polynomial<CT2, o2>& p2) {
   Polynomial<decltype(std::declval<CT1>() * std::declval<CT2>()), o1 + o2> res;
-  for (UInt i = 0; i < p1.terms; ++i)
-    for (UInt j = 0; j < p2.terms; ++j)
-      res.coeffs[i + j] += p1.coeffs[i] * p2.coeffs[j];
+  detail::coeff_product<p1.terms - 1, p2.terms - 1>::prod(res.coeffs, p1.coeffs,
+                                                          p2.coeffs);
   return res;
 }
 
+namespace detail {
+template <std::size_t i>
+struct coeff_sum {
+  template <typename CT1, typename CT2, typename CT3, std::size_t n,
+            std::size_t m, std::size_t l>
+  static constexpr void sum(std::array<CT1, n>& res,
+                            const std::array<CT2, m>& a,
+                            const std::array<CT3, l>& b) {
+    std::get<i>(res) += ((i < m) ? a[i] : 0) + ((i < l) ? b[i] : 0);
+    coeff_sum<i-1>::sum(res, a, b);
+  }
+};
+
+template <>
+struct coeff_sum<0> {
+  template <typename CT1, typename CT2, typename CT3, std::size_t n,
+            std::size_t m, std::size_t l>
+  static constexpr void sum(std::array<CT1, n>& res,
+                            const std::array<CT2, m>& a,
+                            const std::array<CT3, l>& b) {
+    std::get<0>(res) += a[0] + b[0];
+  }
+};
+}  // namespace detail
+
 /// Sum of polynomials
 template <typename CT1, typename CT2, UInt o1, UInt o2>
 constexpr auto operator+(const Polynomial<CT1, o1>& p1,
                          const Polynomial<CT2, o2>& p2) {
   Polynomial<decltype(std::declval<CT1>() + std::declval<CT2>()),
              std::max(o1, o2)>
       res;
-  for (UInt i = 0; i < res.terms; ++i)
-    res.coeffs[i] += ((i < p1.terms) ? p1.coeffs[i] : 0) +
-                     ((i < p2.terms) ? p2.coeffs[i] : 0);
+  detail::coeff_sum<res.terms - 1>::sum(res.coeffs, p1.coeffs, p2.coeffs);
   return res;
 }
 
 /* -------------------------------------------------------------------------- */
 /* Simplification rules */
 /* -------------------------------------------------------------------------- */
 
 #define POLY_SIMPLIFICATION(opsign, std_op)                                    \
   template <typename CT1, typename CT2, UInt o1, UInt o2, typename Derived>    \
   constexpr auto operator opsign(                                              \
       const Polynomial<CT1, o1>& p1,                                           \
       const BinaryOperation<Polynomial<CT2, o2>, Derived, std_op>& p2) {       \
     return p2.op(p2.op(p1, p2.operands.first), p2.operands.second);            \
   }                                                                            \
                                                                                \
   template <typename CT1, typename CT2, UInt o1, UInt o2, typename Derived>    \
   constexpr auto operator opsign(                                              \
       const BinaryOperation<Polynomial<CT2, o2>, Derived, std_op>& p2,         \
       const Polynomial<CT1, o1>& p1) {                                         \
     return p2.op(p1, p2);                                                      \
   }                                                                            \
                                                                                \
   template <typename CT1, typename CT2, UInt o1, UInt o2, typename Derived>    \
   constexpr auto operator opsign(                                              \
       const Polynomial<CT1, o1>& p1,                                           \
       const SBinaryOperation<Polynomial<CT2, o2>, Derived, std_op>& p2) {      \
     return p2.op(p1, p2.commute());                                            \
   }                                                                            \
                                                                                \
   template <typename CT1, typename CT2, UInt o1, UInt o2, typename Derived>    \
   constexpr auto operator opsign(                                              \
       const SBinaryOperation<Polynomial<CT2, o2>, Derived, std_op>& p2,        \
       const Polynomial<CT1, o1>& p1) {                                         \
     return p2.op(p1, p2.commute());                                            \
   }
 
-POLY_SIMPLIFICATION(+, std::plus<>);
-POLY_SIMPLIFICATION(*, std::multiplies<>);
+POLY_SIMPLIFICATION(+, std::plus<>)
+POLY_SIMPLIFICATION(*, std::multiplies<>)
 
 #undef POLY_SIMPLIFICATION
 
 template <typename CT1, UInt o, typename CT2, typename Derived>
 constexpr auto operator/(const Product<Polynomial<CT1, o>, Derived>& p,
-			 const Constant<CT2>& c) {
+                         const Constant<CT2>& c) {
   return p.operands.second * (p.operands.first / c);
 }
 
 template <typename CT1, UInt o, typename CT2, typename Derived>
 constexpr auto operator/(const Product<Derived, Polynomial<CT1, o>>& p,
-			 const Constant<CT2>& c) {
+                         const Constant<CT2>& c) {
   return p.operands.second * (p.operands.first / c);
 }
 
 template <typename CT1, UInt o, typename CT2, typename Derived>
 constexpr auto operator*(const Polynomial<CT1, o>& p,
-			 const Division<Derived, Constant<CT2>>& d) {
+                         const Division<Derived, Constant<CT2>>& d) {
   return (p / d.operands.second) * d.operands.first;
 }
 
 template <typename CT1, UInt o, typename CT2, typename Derived>
 constexpr auto operator*(const Division<Derived, Constant<CT2>>& d,
                          const Polynomial<CT1, o>& p) {
   return (p / d.operands.second) * d.operands.first;
 }
 
 /* -------------------------------------------------------------------------- */
 /* Output */
 /* -------------------------------------------------------------------------- */
 
 template <typename CT, UInt order>
 std::ostream& operator<<(std::ostream& os, const Polynomial<CT, order>& p) {
   os << p.coeffs[0];
 
   if (p.terms >= 2) {
     os << " + " << p.coeffs[1] << "*z";
     for (UInt i = 2; i < p.terms; ++i)
       os << " + " << p.coeffs[i] << "*z^" << i;
   }
   return os;
 }
 
 }  // namespace expolit
 #endif
diff --git a/tests/SConstruct b/tests/SConstruct
index 2601785..ff4fe1b 100644
--- a/tests/SConstruct
+++ b/tests/SConstruct
@@ -1,7 +1,7 @@
 import os
 
-env = DefaultEnvironment(ENV=os.environ)
-env.MergeFlags('-std=c++17 -g -O0')
+env = DefaultEnvironment(ENV=os.environ, CXX='g++-6')
+env.MergeFlags('-std=c++14 -g -O0 -Wall -Wextra -pedantic')
 env.Append(LIBS=['pthread', 'gtest', 'gtest_main'],
            CPPPATH=['../include'])
 env.Program('tests_gtest', ['tests_gtest.cc'])