diff --git a/dynare++/tl/testing/monoms.cc b/dynare++/tl/testing/monoms.cc
index a8a2c269d3fa118005f340a246fe671a3f3387fb..737fed31c47518c937984e180d4d660da822d977 100644
--- a/dynare++/tl/testing/monoms.cc
+++ b/dynare++/tl/testing/monoms.cc
@@ -5,9 +5,7 @@
 #include "tl_exception.hh"
 #include "fs_tensor.hh"
 
-#include <cstdlib>
-#include <cmath>
-#include <cstdio>
+#include <iostream>
 
 IntGenerator intgen;
 
@@ -17,23 +15,23 @@ IntGenerator::init(int nf, int ny, int nv, int nw, int nu,
 {
   maxim = mx;
   probab = prob;
-  long int seed = nf;
+  decltype(mtgen)::result_type seed = nf;
   seed = 256*seed + ny;
   seed = 256*seed + nv;
   seed = 256*seed + nw;
   seed = 256*seed + nu;
-  srand48(seed);
+  mtgen.seed(seed);
 }
 
 int
-IntGenerator::get() const
+IntGenerator::get()
 {
-  double d = drand48();
-  auto num_inter = (int) (((double) 2*maxim)/(1.0-probab));
+  double d = dis(mtgen);
+  auto num_inter = static_cast<int>(static_cast<double>(2*maxim)/(1.0-probab));
   int num_zero_inter = num_inter - 2*maxim;
-  if (d < ((double) num_zero_inter)/num_inter)
+  if (d < static_cast<double>(num_zero_inter)/num_inter)
     return 0;
-  return (int) (d*num_inter)-num_zero_inter-maxim;
+  return static_cast<int>((d*num_inter)-num_zero_inter-maxim);
 }
 
 Monom::Monom(int len)
@@ -81,28 +79,17 @@ Monom::multiplyWith(int ex, const Monom &m)
 void
 Monom::print() const
 {
-  printf("[");
+  std::cout << '[';
   for (int i = 0; i < size(); i++)
-    printf("%3d", operator[](i));
-  printf("]");
+    std::cout << operator[](i);
+  std::cout << ']';
 }
 
 Monom1Vector::Monom1Vector(int nxx, int l)
-  : nx(nxx), len(l), x(new Monom *[len])
+  : nx(nxx), len(l)
 {
   for (int i = 0; i < len; i++)
-    {
-      x[i] = new Monom(nx);
-    }
-}
-
-Monom1Vector::~Monom1Vector()
-{
-  for (int i = 0; i < len; i++)
-    {
-      delete x[i];
-    }
-  delete [] x;
+    x.emplace_back(nx);
 }
 
 void
@@ -112,9 +99,7 @@ Monom1Vector::deriv(const IntSequence &c, Vector &out) const
               "Wrong length of output vector in Monom1Vector::deriv");
 
   for (int i = 0; i < len; i++)
-    {
-      out[i] = x[i]->deriv(c);
-    }
+    out[i] = x[i].deriv(c);
 }
 
 FGSTensor *
@@ -133,60 +118,46 @@ Monom1Vector::deriv(int dim) const
 void
 Monom1Vector::print() const
 {
-  printf("Variables: x(%d)\n", nx);
-  printf("Rows: %d\n", len);
+  std::cout << "Variables: x(" << nx << ")\n"
+            << "Rows: " << len << '\n';
   for (int i = 0; i < len; i++)
     {
-      printf("%2d: ", i);
-      x[i]->print();
-      printf("\n");
+      std::cout << i << ": ";
+      x[i].print();
+      std::cout << '\n';
     }
 }
 
 Monom2Vector::Monom2Vector(int nyy, int nuu, int l)
-  : ny(nyy), nu(nuu), len(l), y(new Monom *[len]), u(new Monom *[len])
+  : ny(nyy), nu(nuu), len(l)
 {
   for (int i = 0; i < len; i++)
     {
-      y[i] = new Monom(ny);
-      u[i] = new Monom(nu);
+      y.emplace_back(ny);
+      u.emplace_back(nu);
     }
 }
 
 Monom2Vector::Monom2Vector(const Monom1Vector &g, const Monom2Vector &xmon)
-  : ny(xmon.ny), nu(xmon.nu), len(g.len),
-    y(new Monom *[len]), u(new Monom *[len])
+  : ny(xmon.ny), nu(xmon.nu), len(g.len)
 {
   TL_RAISE_IF(xmon.len != g.nx,
               "Wrong number of x's in Monom2Vector constructor");
 
   for (int i = 0; i < len; i++)
     {
-      y[i] = new Monom(ny, 0);
-      u[i] = new Monom(nu, 0);
+      y.emplace_back(ny, 0);
+      u.emplace_back(nu, 0);
     }
 
   for (int i = 0; i < len; i++)
-    {
-      // multiply from xmon
-      for (int j = 0; j < g.nx; j++)
-        {
-          int ex = g.x[i]->operator[](j);
-          y[i]->multiplyWith(ex, *(xmon.y[j]));
-          u[i]->multiplyWith(ex, *(xmon.u[j]));
-        }
-    }
-}
-
-Monom2Vector::~Monom2Vector()
-{
-  for (int i = 0; i < len; i++)
-    {
-      delete y[i];
-      delete u[i];
-    }
-  delete [] y;
-  delete [] u;
+    // multiply from xmon
+    for (int j = 0; j < g.nx; j++)
+      {
+        int ex = g.x[i].operator[](j);
+        y[i].multiplyWith(ex, xmon.y[j]);
+        u[i].multiplyWith(ex, xmon.u[j]);
+      }
 }
 
 void
@@ -202,15 +173,13 @@ Monom2Vector::deriv(const Symmetry &s, const IntSequence &c,
   IntSequence cy(c, 0, s[0]);
   IntSequence cu(c, s[0], s.dimen());
   for (int i = 0; i < len; i++)
-    {
-      out[i] = y[i]->deriv(cy) * u[i]->deriv(cu);
-    }
+    out[i] = y[i].deriv(cy) * u[i].deriv(cu);
 }
 
 FGSTensor *
 Monom2Vector::deriv(const Symmetry &s) const
 {
-  IntSequence nvs(2); nvs[0] = ny; nvs[1] = nu;
+  IntSequence nvs{ny, nu};
   FGSTensor *t = new FGSTensor(len, TensorDimens(s, nvs));
   for (Tensor::index it = t->begin(); it != t->end(); ++it)
     {
@@ -225,96 +194,63 @@ Monom2Vector::deriv(int maxdim) const
 {
   auto *res = new FGSContainer(2);
   for (int dim = 1; dim <= maxdim; dim++)
-    {
-      for (int ydim = 0; ydim <= dim; ydim++)
-        {
-          int udim = dim - ydim;
-          Symmetry s{ydim, udim};
-          res->insert(deriv(s));
-        }
-    }
+    for (int ydim = 0; ydim <= dim; ydim++)
+      {
+        int udim = dim - ydim;
+        Symmetry s{ydim, udim};
+        res->insert(deriv(s));
+      }
   return res;
 }
 
 void
 Monom2Vector::print() const
 {
-  printf("Variables: y(%d) u(%d)\n", ny, nu);
-  printf("Rows: %d\n", len);
+  std::cout << "Variables: y(" << ny << ") u(" << nu << ")\n"
+            << "Rows: " << len << '\n';
   for (int i = 0; i < len; i++)
     {
-      printf("%2d: ", i);
-      y[i]->print();
-      printf("    ");
-      u[i]->print();
-      printf("\n");
+      std::cout << i;
+      y[i].print();
+      std::cout << "    ";
+      u[i].print();
+      std::cout << '\n';
     }
 }
 
-Monom4Vector::~Monom4Vector()
-{
-  for (int i = 0; i < len; i++)
-    {
-      delete x1[i];
-      delete x2[i];
-      delete x3[i];
-      delete x4[i];
-    }
-  delete [] x1;
-  delete [] x2;
-  delete [] x3;
-  delete [] x4;
-}
-
 void
 Monom4Vector::init_random()
 {
   for (int i = 0; i < len; i++)
     {
-      x1[i] = new Monom(nx1);
-      x2[i] = new Monom(nx2);
-      x3[i] = new Monom(nx3);
-      x4[i] = new Monom(nx4);
+      x1.emplace_back(nx1);
+      x2.emplace_back(nx2);
+      x3.emplace_back(nx3);
+      x4.emplace_back(nx4);
     }
 }
 
 Monom4Vector::Monom4Vector(int l, int ny, int nu)
-  : len(l), nx1(ny), nx2(nu), nx3(0), nx4(1),
-    x1(new Monom *[len]),
-    x2(new Monom *[len]),
-    x3(new Monom *[len]),
-    x4(new Monom *[len])
+  : len(l), nx1(ny), nx2(nu), nx3(0), nx4(1)
 {
   init_random();
 }
 
 Monom4Vector::Monom4Vector(int l, int ny, int nu, int nup)
-  : len(l), nx1(ny), nx2(nu), nx3(nup), nx4(1),
-    x1(new Monom *[len]),
-    x2(new Monom *[len]),
-    x3(new Monom *[len]),
-    x4(new Monom *[len])
+  : len(l), nx1(ny), nx2(nu), nx3(nup), nx4(1)
 {
   init_random();
 }
 
 Monom4Vector::Monom4Vector(int l, int nbigg, int ng, int ny, int nu)
-  : len(l), nx1(nbigg), nx2(ng), nx3(ny), nx4(nu),
-    x1(new Monom *[len]),
-    x2(new Monom *[len]),
-    x3(new Monom *[len]),
-    x4(new Monom *[len])
+  : len(l), nx1(nbigg), nx2(ng), nx3(ny), nx4(nu)
 {
   init_random();
 }
 
 Monom4Vector::Monom4Vector(const Monom4Vector &f, const Monom4Vector &bigg,
                            const Monom4Vector &g)
-  : len(f.len), nx1(bigg.nx1), nx2(bigg.nx2), nx3(bigg.nx3), nx4(1),
-    x1(new Monom *[len]),
-    x2(new Monom *[len]),
-    x3(new Monom *[len]),
-    x4(new Monom *[len])
+  : len(f.len), nx1(bigg.nx1), nx2(bigg.nx2), nx3(bigg.nx3), nx4(1)
 {
   TL_RAISE_IF(!(bigg.nx1 == g.nx1 && bigg.nx2 == g.nx2 && g.nx3 == 0
                 && bigg.nx4 == 1 && g.nx4 == 1),
@@ -325,10 +261,10 @@ Monom4Vector::Monom4Vector(const Monom4Vector &f, const Monom4Vector &bigg,
 
   for (int i = 0; i < len; i++)
     {
-      x1[i] = new Monom(nx1, 0);
-      x2[i] = new Monom(nx2, 0);
-      x3[i] = new Monom(nx3, 0);
-      x4[i] = new Monom(nx4, 0);
+      x1.emplace_back(nx1, 0);
+      x2.emplace_back(nx2, 0);
+      x3.emplace_back(nx3, 0);
+      x4.emplace_back(nx4, 0);
     }
 
   for (int i = 0; i < len; i++)
@@ -336,24 +272,24 @@ Monom4Vector::Monom4Vector(const Monom4Vector &f, const Monom4Vector &bigg,
       // multiply from G (first argument)
       for (int j = 0; j < f.nx1; j++)
         {
-          int ex = f.x1[i]->operator[](j);
-          x1[i]->multiplyWith(ex, *(bigg.x1[j]));
-          x2[i]->multiplyWith(ex, *(bigg.x2[j]));
-          x3[i]->multiplyWith(ex, *(bigg.x3[j]));
-          x4[i]->multiplyWith(ex, *(bigg.x4[j]));
+          int ex = f.x1[i].operator[](j);
+          x1[i].multiplyWith(ex, bigg.x1[j]);
+          x2[i].multiplyWith(ex, bigg.x2[j]);
+          x3[i].multiplyWith(ex, bigg.x3[j]);
+          x4[i].multiplyWith(ex, bigg.x4[j]);
         }
       // multiply from g (second argument)
       for (int j = 0; j < f.nx2; j++)
         {
-          int ex = f.x2[i]->operator[](j);
-          x1[i]->multiplyWith(ex, *(g.x1[j]));
-          x2[i]->multiplyWith(ex, *(g.x2[j]));
-          x4[i]->multiplyWith(ex, *(g.x4[j]));
+          int ex = f.x2[i].operator[](j);
+          x1[i].multiplyWith(ex, g.x1[j]);
+          x2[i].multiplyWith(ex, g.x2[j]);
+          x4[i].multiplyWith(ex, g.x4[j]);
         }
       // add y as third argument of f
-      x1[i]->add(1, *(f.x3[i]));
+      x1[i].add(1, f.x3[i]);
       // add u as fourth argument of f
-      x2[i]->add(1, *(f.x4[i]));
+      x2[i].add(1, f.x4[i]);
     }
 }
 
@@ -372,23 +308,20 @@ Monom4Vector::deriv(const Symmetry &s, const IntSequence &coor,
     {
       out[i] = 1;
       int off = 0;
-      out[i] *= x1[i]->deriv(IntSequence(coor, off, off+s[0]));
+      out[i] *= x1[i].deriv(IntSequence(coor, off, off+s[0]));
       off += s[0];
-      out[i] *= x2[i]->deriv(IntSequence(coor, off, off+s[1]));
+      out[i] *= x2[i].deriv(IntSequence(coor, off, off+s[1]));
       off += s[1];
-      out[i] *= x3[i]->deriv(IntSequence(coor, off, off+s[2]));
+      out[i] *= x3[i].deriv(IntSequence(coor, off, off+s[2]));
       off += s[2];
-      out[i] *= x4[i]->deriv(IntSequence(coor, off, off+s[3]));
+      out[i] *= x4[i].deriv(IntSequence(coor, off, off+s[3]));
     }
 }
 
 FGSTensor *
 Monom4Vector::deriv(const Symmetry &s) const
 {
-  IntSequence nvs(4);
-  nvs[0] = nx1; nvs[1] = nx2;
-  nvs[2] = nx3; nvs[3] = nx4;
-
+  IntSequence nvs{nx1, nx2, nx3, nx4};
   FGSTensor *res = new FGSTensor(len, TensorDimens(s, nvs));
   for (Tensor::index run = res->begin(); run != res->end(); ++run)
     {
@@ -401,8 +334,7 @@ Monom4Vector::deriv(const Symmetry &s) const
 FSSparseTensor *
 Monom4Vector::deriv(int dim) const
 {
-  IntSequence cum(4);
-  cum[0] = 0; cum[1] = nx1; cum[2] = nx1+nx2; cum[3] = nx1+nx2+nx3;
+  IntSequence cum{0, nx1, nx1+nx2, nx1+nx2+nx3};
 
   auto *res = new FSSparseTensor(dim, nx1+nx2+nx3+nx4, len);
 
@@ -423,12 +355,8 @@ Monom4Vector::deriv(int dim) const
       Vector col(len);
       deriv(ind_sym, ind, col);
       for (int i = 0; i < len; i++)
-        {
-          if (col[i] != 0.0)
-            {
-              res->insert(run.getCoor(), i, col[i]);
-            }
-        }
+        if (col[i] != 0.0)
+          res->insert(run.getCoor(), i, col[i]);
     }
 
   return res;
@@ -437,27 +365,26 @@ Monom4Vector::deriv(int dim) const
 void
 Monom4Vector::print() const
 {
-  printf("Variables: x1(%d) x2(%d) x3(%d) x4(%d)\n",
-         nx1, nx2, nx3, nx4);
-  printf("Rows: %d\n", len);
+  std::cout << "Variables: x1(" << nx1 << ") x2(" << nx2
+            << ") x3(" << nx3 << ") x4(" << nx4 << ")\n"
+            << "Rows: " << len << '\n';
   for (int i = 0; i < len; i++)
     {
-      printf("%2d: ", i);
-      x1[i]->print();
-      printf("    ");
-      x2[i]->print();
-      printf("    ");
-      x3[i]->print();
-      printf("    ");
-      x4[i]->print();
-      printf("\n");
+      std::cout << i << ": ";
+      x1[i].print();
+      std::cout << "    ";
+      x2[i].print();
+      std::cout << "    ";
+      x3[i].print();
+      std::cout << "    ";
+      x4[i].print();
+      std::cout << '\n';
     }
 }
 
-SparseDerivGenerator::SparseDerivGenerator(
-                                           int nf, int ny, int nu, int nup, int nbigg, int ng,
+SparseDerivGenerator::SparseDerivGenerator(int nf, int ny, int nu, int nup, int nbigg, int ng,
                                            int mx, double prob, int maxdim)
-  : maxdimen(maxdim), ts(new FSSparseTensor *[maxdimen])
+  : maxdimen(maxdim), bigg(4), g(4), rcont(4), ts(new FSSparseTensor *[maxdimen])
 {
   intgen.init(nf, ny, nu, nup, nbigg, mx, prob);
 
@@ -465,18 +392,15 @@ SparseDerivGenerator::SparseDerivGenerator(
   Monom4Vector g_m(ng, ny, nu);
   Monom4Vector f(nf, nbigg, ng, ny, nu);
   Monom4Vector r(f, bigg_m, g_m);
-  bigg = new FGSContainer(4);
-  g = new FGSContainer(4);
-  rcont = new FGSContainer(4);
 
   for (int dim = 1; dim <= maxdimen; dim++)
     {
       for (auto &si : SymmetrySet(dim, 4))
         {
-          bigg->insert(bigg_m.deriv(si));
-          rcont->insert(r.deriv(si));
+          bigg.insert(bigg_m.deriv(si));
+          rcont.insert(r.deriv(si));
           if (si[2] == 0)
-            g->insert(g_m.deriv(si));
+            g.insert(g_m.deriv(si));
         }
 
       ts[dim-1] = f.deriv(dim);
@@ -485,9 +409,6 @@ SparseDerivGenerator::SparseDerivGenerator(
 
 SparseDerivGenerator::~SparseDerivGenerator()
 {
-  delete bigg;
-  delete g;
-  delete rcont;
   for (int i = 0; i < maxdimen; i++)
     delete ts[i];
   delete [] ts;
@@ -517,9 +438,7 @@ DenseDerivGenerator::unfold()
 {
   uxcont = new UGSContainer(*xcont);
   for (int i = 0; i < maxdimen; i++)
-    {
-      uts[i] = new UGSTensor(*(ts[i]));
-    }
+    uts[i] = new UGSTensor(*(ts[i]));
 }
 
 DenseDerivGenerator::~DenseDerivGenerator()
diff --git a/dynare++/tl/testing/monoms.hh b/dynare++/tl/testing/monoms.hh
index 59c30773c9fea48544a2b62757b0fb021e91ad17..d6716bfd5e99ddf7a74c440fb224cae3f9a10f09 100644
--- a/dynare++/tl/testing/monoms.hh
+++ b/dynare++/tl/testing/monoms.hh
@@ -4,6 +4,9 @@
 #ifndef MONOMS_H
 #define MONOMS_H
 
+#include <vector>
+#include <random>
+
 #include "int_sequence.hh"
 #include "gs_tensor.hh"
 #include "t_container.hh"
@@ -14,12 +17,12 @@ class IntGenerator
 {
   int maxim{5};
   double probab{0.3};
+  std::mt19937 mtgen;
+  std::uniform_real_distribution<> dis;
 public:
-  IntGenerator()
-     
-  = default;
+  IntGenerator() = default;
   void init(int nf, int ny, int nv, int nw, int nu, int mx, double prob);
-  int get() const;
+  int get();
 };
 
 extern IntGenerator intgen;
@@ -41,10 +44,10 @@ class Monom1Vector
   friend class Monom2Vector;
   int nx;
   int len;
-  Monom **const x;
+  std::vector<Monom> x;
 public:
   Monom1Vector(int nxx, int l);
-  ~Monom1Vector();
+  ~Monom1Vector() = default;
   void deriv(const IntSequence &c, Vector &out) const;
   FGSTensor *deriv(int dim) const;
   void print() const;
@@ -53,17 +56,15 @@ public:
 //class Monom3Vector;
 class Monom2Vector
 {
-  int ny;
-  int nu;
+  int ny, nu;
   int len;
-  Monom **const y;
-  Monom **const u;
+  std::vector<Monom> y, u;
 public:
   // generate random vector of monom two vector
   Monom2Vector(int nyy, int nuu, int l);
   // calculate g(x(y,u))
   Monom2Vector(const Monom1Vector &g, const Monom2Vector &xmon);
-  ~Monom2Vector();
+  ~Monom2Vector() = default;
   void deriv(const Symmetry &s, const IntSequence &c, Vector &out) const;
   FGSTensor *deriv(const Symmetry &s) const;
   FGSContainer *deriv(int maxdim) const;
@@ -73,14 +74,8 @@ public:
 class Monom4Vector
 {
   int len;
-  int nx1;
-  int nx2;
-  int nx3;
-  int nx4;
-  Monom **const x1;
-  Monom **const x2;
-  Monom **const x3;
-  Monom **const x4;
+  int nx1, nx2, nx3, nx4;
+  std::vector<Monom> x1, x2, x3, x4;
 public:
   /* random for g(y,u,sigma) */
   Monom4Vector(int l, int ny, int nu);
@@ -91,7 +86,7 @@ public:
   /* substitution f(G(y,u,u',sigma),g(y,u,sigma),y,u) */
   Monom4Vector(const Monom4Vector &f, const Monom4Vector &bigg,
                const Monom4Vector &g);
-  ~Monom4Vector();
+  ~Monom4Vector() = default;
   FSSparseTensor *deriv(int dim) const;
   FGSTensor *deriv(const Symmetry &s) const;
   void deriv(const Symmetry &s, const IntSequence &coor, Vector &out) const;
@@ -103,9 +98,9 @@ protected:
 struct SparseDerivGenerator
 {
   int maxdimen;
-  FGSContainer *bigg;
-  FGSContainer *g;
-  FGSContainer *rcont;
+  FGSContainer bigg;
+  FGSContainer g;
+  FGSContainer rcont;
   FSSparseTensor **const ts;
   SparseDerivGenerator(int nf, int ny, int nu, int nup, int nbigg, int ng,
                        int mx, double prob, int maxdim);
diff --git a/dynare++/tl/testing/tests.cc b/dynare++/tl/testing/tests.cc
index 3bb93092d4a51842eb042c0535ee268fa0299a52..2b19ffa39bf054e697aad01bd939096410f4c639 100644
--- a/dynare++/tl/testing/tests.cc
+++ b/dynare++/tl/testing/tests.cc
@@ -13,28 +13,28 @@
 #include "ps_tensor.hh"
 #include "tl_static.hh"
 
-#include <cstdio>
-#include <cstring>
+#include <string>
+#include <utility>
+#include <memory>
+#include <vector>
 #include <ctime>
+#include <cstdlib>
+#include <iostream>
+#include <iomanip>
 
 class TestRunnable
 {
-  char name[100];
 public:
+  const std::string name;
   int dim; // dimension of the solved problem
   int nvar; // number of variable of the solved problem
-  TestRunnable(const char *n, int d, int nv)
-    : dim(d), nvar(nv)
+  TestRunnable(std::string name_arg, int d, int nv)
+    : name{std::move(name_arg)}, dim(d), nvar(nv)
   {
-    strncpy(name, n, 100);
   }
+  virtual ~TestRunnable() = default;
   bool test() const;
   virtual bool run() const = 0;
-  const char *
-  getName() const
-  {
-    return name;
-  }
 protected:
   template<class _Ttype>
   static bool index_forward(const Symmetry &s, const IntSequence &nvs);
@@ -92,22 +92,17 @@ protected:
 bool
 TestRunnable::test() const
 {
-  printf("Running test <%s>\n", name);
+  std::cout << "Running test <" << name << ">" << std::endl;
   clock_t start = clock();
   bool passed = run();
   clock_t end = clock();
-  printf("CPU time %8.4g (CPU seconds)..................",
-         ((double) (end-start))/CLOCKS_PER_SEC);
+  std::cout << "CPU time " << static_cast<double>(end-start)/CLOCKS_PER_SEC
+            << " (CPU seconds)..................";
   if (passed)
-    {
-      printf("passed\n\n");
-      return passed;
-    }
+    std::cout << "passed\n\n";
   else
-    {
-      printf("FAILED\n\n");
-      return passed;
-    }
+    std::cout << "FAILED\n\n";
+  return passed;
 }
 
 /****************************************************/
@@ -137,10 +132,10 @@ TestRunnable::index_forward(const Symmetry &s, const IntSequence &nvs)
     }
   while (run != dummy.begin());
 
-  printf("\tnumber of columns    = %d\n", dummy.ncols());
-  printf("\tnumber of increments = %d\n", nincr);
-  printf("\tnumber of decrements = %d\n", ndecr);
-  printf("\tnumber of failures   = %d\n", fails);
+  std::cout << "\tnumber of columns    = " << dummy.ncols() << '\n'
+            << "\tnumber of increments = " << nincr << '\n'
+            << "\tnumber of decrements = " << ndecr << '\n'
+            << "\tnumber of failures   = " << fails << '\n';
 
   return fails == 0;
 }
@@ -168,10 +163,10 @@ TestRunnable::index_backward(const Symmetry &s, const IntSequence &nvs)
       nincr++;
     }
 
-  printf("\tnumber of columns    = %d\n", dummy.ncols());
-  printf("\tnumber of increments = %d\n", nincr);
-  printf("\tnumber of decrements = %d\n", ndecr);
-  printf("\tnumber of failures   = %d\n", fails);
+  std::cout << "\tnumber of columns    = " << dummy.ncols() << '\n'
+            << "\tnumber of increments = " << nincr << '\n'
+            << "\tnumber of decrements = " << ndecr << '\n'
+            << "\tnumber of failures   = " << fails << '\n';
 
   return fails == 0;
 }
@@ -191,9 +186,9 @@ TestRunnable::index_offset(const Symmetry &s, const IntSequence &nvs)
         fails++;
     }
 
-  printf("\tnumber of columns    = %d\n", dummy.ncols());
-  printf("\tnumber of increments = %d\n", nincr);
-  printf("\tnumber of failures   = %d\n", fails);
+  std::cout << "\tnumber of columns    = " << dummy.ncols() << '\n'
+            << "\tnumber of increments = " << nincr << '\n'
+            << "\tnumber of failures   = " << fails << '\n';
 
   return fails == 0;
 }
@@ -206,10 +201,10 @@ TestRunnable::fold_unfold(const FTensor *folded)
   folded2->add(-1.0, *folded);
   double normInf = folded2->getNormInf();
   double norm1 = folded2->getNorm1();
-  printf("\tfolded size:       (%d, %d)\n", folded->nrows(), folded->ncols());
-  printf("\tunfolded size:     (%d, %d)\n", unfolded->nrows(), unfolded->ncols());
-  printf("\tdifference normInf: %8.4g\n", normInf);
-  printf("\tdifference norm1:   %8.4g\n", norm1);
+  std::cout << "\tfolded size:       (" << folded->nrows() << ", " << folded->ncols() << ")\n"
+            << "\tunfolded size:     (" << unfolded->nrows() << ", " << unfolded->ncols() << ")\n"
+            << "\tdifference normInf: " << normInf << '\n'
+            << "\tdifference norm1:   " << norm1 << '\n';
 
   delete folded;
   delete unfolded;
@@ -247,15 +242,12 @@ TestRunnable::dense_prod(const Symmetry &bsym, const IntSequence &bnvs,
   double norm1 = btmp.getNorm1();
   double normInf = btmp.getNormInf();
 
-  printf("\ttime for folded product:     %8.4g\n",
-         ((double) (s2-s1))/CLOCKS_PER_SEC);
-  printf("\ttime for unfolded product:   %8.4g\n",
-         ((double) (s5-s4))/CLOCKS_PER_SEC);
-  printf("\ttime for container convert:  %8.4g\n",
-         ((double) (s3-s2))/CLOCKS_PER_SEC);
-  printf("\tunfolded difference normMax: %10.6g\n", norm);
-  printf("\tunfolded difference norm1:   %10.6g\n", norm1);
-  printf("\tunfolded difference normInf: %10.6g\n", normInf);
+  std::cout << "\ttime for folded product:     " << static_cast<double>(s2-s1)/CLOCKS_PER_SEC << '\n'
+            << "\ttime for unfolded product:   " << static_cast<double>(s5-s4)/CLOCKS_PER_SEC << '\n'
+            << "\ttime for container convert:  " << static_cast<double>(s3-s2)/CLOCKS_PER_SEC << '\n'
+            << "\tunfolded difference normMax: " << norm << '\n'
+            << "\tunfolded difference norm1:   " << norm1 << '\n'
+            << "\tunfolded difference normInf: " << normInf << '\n';
 
   delete cont;
   delete fh;
@@ -269,28 +261,27 @@ TestRunnable::folded_monomial(int ng, int nx, int ny, int nu, int dim)
   clock_t gen_time = clock();
   DenseDerivGenerator gen(ng, nx, ny, nu, 5, 0.3, dim);
   gen_time = clock()-gen_time;
-  printf("\ttime for monom generation: %8.4g\n",
-         ((double) gen_time)/CLOCKS_PER_SEC);
-  IntSequence nvs(2); nvs[0] = ny; nvs[1] = nu;
+  std::cout << "\ttime for monom generation: "
+            << static_cast<double>(gen_time)/CLOCKS_PER_SEC << '\n';
+  IntSequence nvs{ny, nu};
   double maxnorm = 0;
   for (int ydim = 0; ydim <= dim; ydim++)
     {
       Symmetry s{ydim, dim-ydim};
-      printf("\tSymmetry: "); s.print();
+      std::cout << "\tSymmetry: ";
+      s.print();
       FGSTensor res(ng, TensorDimens(s, nvs));
       res.getData().zeros();
       clock_t stime = clock();
       for (int d = 1; d <= dim; d++)
-        {
-          gen.xcont->multAndAdd(*(gen.ts[d-1]), res);
-        }
+        gen.xcont->multAndAdd(*(gen.ts[d-1]), res);
       stime = clock() - stime;
-      printf("\t\ttime for symmetry: %8.4g\n",
-             ((double) stime)/CLOCKS_PER_SEC);
+      std::cout << "\t\ttime for symmetry: "
+                << static_cast<double>(stime)/CLOCKS_PER_SEC << '\n';
       const FGSTensor *mres = gen.rcont->get(s);
       res.add(-1.0, *mres);
       double normtmp = res.getData().getMax();
-      printf("\t\terror normMax:     %10.6g\n", normtmp);
+      std::cout << "\t\terror normMax:     " << normtmp << '\n';
       if (normtmp > maxnorm)
         maxnorm = normtmp;
     }
@@ -303,34 +294,33 @@ TestRunnable::unfolded_monomial(int ng, int nx, int ny, int nu, int dim)
   clock_t gen_time = clock();
   DenseDerivGenerator gen(ng, nx, ny, nu, 5, 0.3, dim);
   gen_time = clock()-gen_time;
-  printf("\ttime for monom generation: %8.4g\n",
-         ((double) gen_time)/CLOCKS_PER_SEC);
+  std::cout << "\ttime for monom generation: "
+            << static_cast<double>(gen_time)/CLOCKS_PER_SEC << '\n';
   clock_t u_time = clock();
   gen.unfold();
   u_time = clock() - u_time;
-  printf("\ttime for monom unfolding:  %8.4g\n",
-         ((double) u_time)/CLOCKS_PER_SEC);
-  IntSequence nvs(2); nvs[0] = ny; nvs[1] = nu;
+  std::cout << "\ttime for monom unfolding:  "
+            << static_cast<double>(u_time)/CLOCKS_PER_SEC << '\n';
+  IntSequence nvs{ny, nu};
   double maxnorm = 0;
   for (int ydim = 0; ydim <= dim; ydim++)
     {
       Symmetry s{ydim, dim-ydim};
-      printf("\tSymmetry: "); s.print();
+      std::cout << "\tSymmetry: ";
+      s.print();
       UGSTensor res(ng, TensorDimens(s, nvs));
       res.getData().zeros();
       clock_t stime = clock();
       for (int d = 1; d <= dim; d++)
-        {
-          gen.uxcont->multAndAdd(*(gen.uts[d-1]), res);
-        }
+        gen.uxcont->multAndAdd(*(gen.uts[d-1]), res);
       stime = clock() - stime;
-      printf("\t\ttime for symmetry: %8.4g\n",
-             ((double) stime)/CLOCKS_PER_SEC);
+      std::cout << "\t\ttime for symmetry: "
+                << static_cast<double>(stime)/CLOCKS_PER_SEC << '\n';
       const FGSTensor *mres = gen.rcont->get(s);
       FGSTensor foldres(res);
       foldres.add(-1.0, *mres);
       double normtmp = foldres.getData().getMax();
-      printf("\t\terror normMax:     %10.6g\n", normtmp);
+      std::cout << "\t\terror normMax:     " << normtmp << '\n';
       if (normtmp > maxnorm)
         maxnorm = normtmp;
     }
@@ -346,42 +336,40 @@ TestRunnable::fold_zcont(int nf, int ny, int nu, int nup, int nbigg,
                           5, 0.55, dim);
   gen_time = clock()-gen_time;
   for (int d = 1; d <= dim; d++)
-    {
-      printf("\tfill of dim=%d tensor:     %3.2f %%\n",
-             d, 100*dg.ts[d-1]->getFillFactor());
-    }
-  printf("\ttime for monom generation: %8.4g\n",
-         ((double) gen_time)/CLOCKS_PER_SEC);
-
-  IntSequence nvs(4);
-  nvs[0] = ny; nvs[1] = nu; nvs[2] = nup; nvs[3] = 1;
+    std::cout << "\tfill of dim=" << d << " tensor:     "
+              << std::setprecision(2) << std::fixed << std::setw(6)
+              << 100*dg.ts[d-1]->getFillFactor()
+              << std::setprecision(6) << std::defaultfloat << " %\n";
+  std::cout << "\ttime for monom generation: "
+            << static_cast<double>(gen_time)/CLOCKS_PER_SEC << '\n';
+
+  IntSequence nvs{ny, nu, nup, 1};
   double maxnorm = 0.0;
 
   // form ZContainer
-  FoldedZContainer zc(dg.bigg, nbigg, dg.g, ng, ny, nu);
+  FoldedZContainer zc(&dg.bigg, nbigg, &dg.g, ng, ny, nu);
 
   for (int d = 2; d <= dim; d++)
-    {
-      for (auto &si : SymmetrySet(d, 4))
-        {
-          printf("\tSymmetry: ");
-          si.print();
-          FGSTensor res(nf, TensorDimens(si, nvs));
-          res.getData().zeros();
-          clock_t stime = clock();
-          for (int l = 1; l <= si.dimen(); l++)
-            zc.multAndAdd(*(dg.ts[l-1]), res);
-          stime = clock() - stime;
-          printf("\t\ttime for symmetry: %8.4g\n",
-                 ((double) stime)/CLOCKS_PER_SEC);
-          const FGSTensor *mres = dg.rcont->get(si);
-          res.add(-1.0, *mres);
-          double normtmp = res.getData().getMax();
-          printf("\t\terror normMax:     %10.6g\n", normtmp);
-          if (normtmp > maxnorm)
-            maxnorm = normtmp;
-        }
-    }
+    for (auto &si : SymmetrySet(d, 4))
+      {
+        std::cout << "\tSymmetry: ";
+        si.print();
+        FGSTensor res(nf, TensorDimens(si, nvs));
+        res.getData().zeros();
+        clock_t stime = clock();
+        for (int l = 1; l <= si.dimen(); l++)
+          zc.multAndAdd(*(dg.ts[l-1]), res);
+        stime = clock() - stime;
+        std::cout << "\t\ttime for symmetry: "
+                  << static_cast<double>(stime)/CLOCKS_PER_SEC << '\n';
+        const FGSTensor *mres = dg.rcont.get(si);
+        res.add(-1.0, *mres);
+        double normtmp = res.getData().getMax();
+        std::cout << "\t\terror normMax:     " << normtmp << '\n';
+        if (normtmp > maxnorm)
+          maxnorm = normtmp;
+      }
+
   return maxnorm < 1.0e-10;
 }
 
@@ -394,50 +382,48 @@ TestRunnable::unfold_zcont(int nf, int ny, int nu, int nup, int nbigg,
                           5, 0.55, dim);
   gen_time = clock()-gen_time;
   for (int d = 1; d <= dim; d++)
-    {
-      printf("\tfill of dim=%d tensor:     %3.2f %%\n",
-             d, 100*dg.ts[d-1]->getFillFactor());
-    }
-  printf("\ttime for monom generation: %8.4g\n",
-         ((double) gen_time)/CLOCKS_PER_SEC);
+    std::cout << "\tfill of dim=" << d << " tensor:     "
+              << std::setprecision(2) << std::fixed << std::setw(6)
+              << 100*dg.ts[d-1]->getFillFactor()
+              << std::setprecision(6) << std::defaultfloat << " %\n";
+  std::cout << "\ttime for monom generation: "
+            << static_cast<double>(gen_time)/CLOCKS_PER_SEC << '\n';
 
   clock_t con_time = clock();
-  UGSContainer uG_cont(*(dg.bigg));
-  UGSContainer ug_cont(*(dg.g));
+  UGSContainer uG_cont(dg.bigg);
+  UGSContainer ug_cont(dg.g);
   con_time = clock()-con_time;
-  printf("\ttime for container unfold: %8.4g\n",
-         ((double) con_time)/CLOCKS_PER_SEC);
+  std::cout << "\ttime for container unfold: "
+            << static_cast<double>(con_time)/CLOCKS_PER_SEC << '\n';
 
-  IntSequence nvs(4);
-  nvs[0] = ny; nvs[1] = nu; nvs[2] = nup; nvs[3] = 1;
+  IntSequence nvs{ny, nu, nup, 1};
   double maxnorm = 0.0;
 
   // form ZContainer
   UnfoldedZContainer zc(&uG_cont, nbigg, &ug_cont, ng, ny, nu);
 
   for (int d = 2; d <= dim; d++)
-    {
-      for (auto &si : SymmetrySet(d, 4))
-        {
-          printf("\tSymmetry: ");
-          si.print();
-          UGSTensor res(nf, TensorDimens(si, nvs));
-          res.getData().zeros();
-          clock_t stime = clock();
-          for (int l = 1; l <= si.dimen(); l++)
-            zc.multAndAdd(*(dg.ts[l-1]), res);
-          stime = clock() - stime;
-          printf("\t\ttime for symmetry: %8.4g\n",
-                 ((double) stime)/CLOCKS_PER_SEC);
-          FGSTensor fold_res(res);
-          const FGSTensor *mres = dg.rcont->get(si);
-          fold_res.add(-1.0, *mres);
-          double normtmp = fold_res.getData().getMax();
-          printf("\t\terror normMax:     %10.6g\n", normtmp);
-          if (normtmp > maxnorm)
-            maxnorm = normtmp;
-        }
-    }
+    for (auto &si : SymmetrySet(d, 4))
+      {
+        std::cout << "\tSymmetry: ";
+        si.print();
+        UGSTensor res(nf, TensorDimens(si, nvs));
+        res.getData().zeros();
+        clock_t stime = clock();
+        for (int l = 1; l <= si.dimen(); l++)
+          zc.multAndAdd(*(dg.ts[l-1]), res);
+        stime = clock() - stime;
+        std::cout << "\t\ttime for symmetry: "
+                  << static_cast<double>(stime)/CLOCKS_PER_SEC << '\n';
+        FGSTensor fold_res(res);
+        const FGSTensor *mres = dg.rcont.get(si);
+        fold_res.add(-1.0, *mres);
+        double normtmp = fold_res.getData().getMax();
+        std::cout << "\t\terror normMax:     " << normtmp << '\n';
+        if (normtmp > maxnorm)
+          maxnorm = normtmp;
+      }
+
   return maxnorm < 1.0e-10;
 }
 
@@ -470,12 +456,12 @@ TestRunnable::folded_contraction(int r, int nv, int dim)
   utime = clock() - utime;
 
   v.add(-1.0, res);
-  printf("\ttime for folded contraction: %8.4g\n",
-         ((double) ctime)/CLOCKS_PER_SEC);
-  printf("\ttime for unfolded power:     %8.4g\n",
-         ((double) utime)/CLOCKS_PER_SEC);
-  printf("\terror normMax:     %10.6g\n", v.getMax());
-  printf("\terror norm1:       %10.6g\n", v.getNorm1());
+  std::cout << "\ttime for folded contraction: "
+            << static_cast<double>(ctime)/CLOCKS_PER_SEC << '\n'
+            << "\ttime for unfolded power:     "
+            << static_cast<double>(utime)/CLOCKS_PER_SEC << '\n'
+            << "\terror normMax:     " << v.getMax() << '\n'
+            << "\terror norm1:       " << v.getNorm1() << '\n';
 
   delete f;
   delete x;
@@ -513,12 +499,12 @@ TestRunnable::unfolded_contraction(int r, int nv, int dim)
   utime = clock() - utime;
 
   v.add(-1.0, res);
-  printf("\ttime for unfolded contraction: %8.4g\n",
-         ((double) ctime)/CLOCKS_PER_SEC);
-  printf("\ttime for unfolded power:       %8.4g\n",
-         ((double) utime)/CLOCKS_PER_SEC);
-  printf("\terror normMax:     %10.6g\n", v.getMax());
-  printf("\terror norm1:       %10.6g\n", v.getNorm1());
+  std::cout << "\ttime for unfolded contraction: "
+            << static_cast<double>(ctime)/CLOCKS_PER_SEC << '\n'
+            << "\ttime for unfolded power:       "
+            << static_cast<double>(utime)/CLOCKS_PER_SEC << '\n'
+            << "\terror normMax:     " << v.getMax() << '\n'
+            << "\terror norm1:       " << v.getNorm1() << '\n';
 
   delete u;
   delete x;
@@ -532,10 +518,14 @@ TestRunnable::poly_eval(int r, int nv, int maxdim)
   Factory fact;
   Vector *x = fact.makeVector(nv);
 
-  Vector out_ft(r); out_ft.zeros();
-  Vector out_fh(r); out_fh.zeros();
-  Vector out_ut(r); out_ut.zeros();
-  Vector out_uh(r); out_uh.zeros();
+  Vector out_ft(r);
+  out_ft.zeros();
+  Vector out_fh(r);
+  out_fh.zeros();
+  Vector out_ut(r);
+  out_ut.zeros();
+  Vector out_uh(r);
+  out_uh.zeros();
 
   UTensorPolynomial *up;
   {
@@ -544,14 +534,14 @@ TestRunnable::poly_eval(int r, int nv, int maxdim)
     clock_t ft_cl = clock();
     fp->evalTrad(out_ft, *x);
     ft_cl = clock() - ft_cl;
-    printf("\ttime for folded power eval:    %8.4g\n",
-           ((double) ft_cl)/CLOCKS_PER_SEC);
+    std::cout << "\ttime for folded power eval:    "
+              << static_cast<double>(ft_cl)/CLOCKS_PER_SEC << '\n';
 
     clock_t fh_cl = clock();
     fp->evalHorner(out_fh, *x);
     fh_cl = clock() - fh_cl;
-    printf("\ttime for folded horner eval:   %8.4g\n",
-           ((double) fh_cl)/CLOCKS_PER_SEC);
+    std::cout << "\ttime for folded horner eval:   "
+              << static_cast<double>(fh_cl)/CLOCKS_PER_SEC << '\n';
 
     up = new UTensorPolynomial(*fp);
     delete fp;
@@ -560,14 +550,14 @@ TestRunnable::poly_eval(int r, int nv, int maxdim)
   clock_t ut_cl = clock();
   up->evalTrad(out_ut, *x);
   ut_cl = clock() - ut_cl;
-  printf("\ttime for unfolded power eval:  %8.4g\n",
-         ((double) ut_cl)/CLOCKS_PER_SEC);
+  std::cout << "\ttime for unfolded power eval:  "
+            << static_cast<double>(ut_cl)/CLOCKS_PER_SEC << '\n';
 
   clock_t uh_cl = clock();
   up->evalHorner(out_uh, *x);
   uh_cl = clock() - uh_cl;
-  printf("\ttime for unfolded horner eval: %8.4g\n",
-         ((double) uh_cl)/CLOCKS_PER_SEC);
+  std::cout << "\ttime for unfolded horner eval: "
+            << static_cast<double>(uh_cl)/CLOCKS_PER_SEC << '\n';
 
   out_ft.add(-1.0, out_ut);
   double max_ft = out_ft.getMax();
@@ -576,9 +566,9 @@ TestRunnable::poly_eval(int r, int nv, int maxdim)
   out_uh.add(-1.0, out_ut);
   double max_uh = out_uh.getMax();
 
-  printf("\tfolded power error norm max:     %10.6g\n", max_ft);
-  printf("\tfolded horner error norm max:    %10.6g\n", max_fh);
-  printf("\tunfolded horner error norm max:  %10.6g\n", max_uh);
+  std::cout << "\tfolded power error norm max:     " << max_ft << '\n'
+            << "\tfolded horner error norm max:    " << max_fh << '\n'
+            << "\tunfolded horner error norm max:  " << max_uh << '\n';
 
   delete up;
   delete x;
@@ -599,7 +589,7 @@ public:
   run() const override
   {
     Symmetry s{2, 3};
-    IntSequence nvs(2); nvs[0] = 4; nvs[1] = 2;
+    IntSequence nvs{4, 2};
     return index_forward<FGSTensor>(s, nvs);
   }
 };
@@ -615,7 +605,7 @@ public:
   run() const override
   {
     Symmetry s{2, 3};
-    IntSequence nvs(2); nvs[0] = 4; nvs[1] = 2;
+    IntSequence nvs{4, 2};
     return index_forward<UGSTensor>(s, nvs);
   }
 };
@@ -631,7 +621,7 @@ public:
   run() const override
   {
     Symmetry s{2, 3, 2};
-    IntSequence nvs(3); nvs[0] = 5; nvs[1] = 2; nvs[2] = 2;
+    IntSequence nvs{5, 2, 2};
     return index_forward<FGSTensor>(s, nvs);
   }
 };
@@ -647,7 +637,7 @@ public:
   run() const override
   {
     Symmetry s{2, 3, 2};
-    IntSequence nvs(3); nvs[0] = 5; nvs[1] = 2; nvs[2] = 2;
+    IntSequence nvs{5, 2, 2};
     return index_forward<UGSTensor>(s, nvs);
   }
 };
@@ -663,7 +653,7 @@ public:
   run() const override
   {
     Symmetry s{1, 1, 3};
-    IntSequence nvs(3); nvs[0] = 3; nvs[1] = 3; nvs[2] = 2;
+    IntSequence nvs{3, 3, 2};
     return index_backward<FGSTensor>(s, nvs);
   }
 };
@@ -679,7 +669,7 @@ public:
   run() const override
   {
     Symmetry s{2, 3, 2};
-    IntSequence nvs(3); nvs[0] = 4; nvs[1] = 2; nvs[2] = 4;
+    IntSequence nvs{4, 2, 4};
     return index_backward<FGSTensor>(s, nvs);
   }
 };
@@ -695,7 +685,7 @@ public:
   run() const override
   {
     Symmetry s{1, 1, 3};
-    IntSequence nvs(3); nvs[0] = 3; nvs[1] = 3; nvs[2] = 2;
+    IntSequence nvs{3, 3, 2};
     return index_backward<UGSTensor>(s, nvs);
   }
 };
@@ -711,7 +701,7 @@ public:
   run() const override
   {
     Symmetry s{2, 3, 2};
-    IntSequence nvs(3); nvs[0] = 4; nvs[1] = 2; nvs[2] = 4;
+    IntSequence nvs{4, 2, 4};
     return index_backward<UGSTensor>(s, nvs);
   }
 };
@@ -727,7 +717,7 @@ public:
   run() const override
   {
     Symmetry s{2, 3};
-    IntSequence nvs(2); nvs[0] = 4; nvs[1] = 2;
+    IntSequence nvs{4, 2};
     return index_offset<FGSTensor>(s, nvs);
   }
 };
@@ -743,7 +733,7 @@ public:
   run() const override
   {
     Symmetry s{2, 3};
-    IntSequence nvs(2); nvs[0] = 4; nvs[1] = 2;
+    IntSequence nvs{4, 2};
     return index_offset<UGSTensor>(s, nvs);
   }
 };
@@ -759,7 +749,7 @@ public:
   run() const override
   {
     Symmetry s{2, 3, 2};
-    IntSequence nvs(3); nvs[0] = 5; nvs[1] = 2; nvs[2] = 2;
+    IntSequence nvs{5, 2, 2};
     return index_offset<FGSTensor>(s, nvs);
   }
 };
@@ -775,7 +765,7 @@ public:
   run() const override
   {
     Symmetry s{2, 3, 2};
-    IntSequence nvs(3); nvs[0] = 5; nvs[1] = 2; nvs[2] = 2;
+    IntSequence nvs{5, 2, 2};
     return index_offset<UGSTensor>(s, nvs);
   }
 };
@@ -805,7 +795,7 @@ public:
   run() const override
   {
     Symmetry s{1, 2, 2};
-    IntSequence nvs(3); nvs[0] = 3; nvs[1] = 3; nvs[2] = 2;
+    IntSequence nvs{3, 3, 2};
     return gs_fold_unfold(5, s, nvs);
   }
 };
@@ -835,7 +825,7 @@ public:
   run() const override
   {
     Symmetry s{2, 1, 2};
-    IntSequence nvs(3); nvs[0] = 6; nvs[1] = 2; nvs[2] = 6;
+    IntSequence nvs{6, 2, 6};
     return gs_fold_unfold(5, s, nvs);
   }
 };
@@ -878,7 +868,7 @@ public:
   bool
   run() const override
   {
-    IntSequence bnvs(2); bnvs[0] = 3; bnvs[1] = 2;
+    IntSequence bnvs{3, 2};
     return dense_prod(Symmetry{1, 2}, bnvs, 2, 3, 2);
   }
 };
@@ -893,7 +883,7 @@ public:
   bool
   run() const override
   {
-    IntSequence bnvs(2); bnvs[0] = 10; bnvs[1] = 7;
+    IntSequence bnvs{10, 7};
     return dense_prod(Symmetry{2, 3}, bnvs, 3, 15, 10);
   }
 };
@@ -908,7 +898,7 @@ public:
   bool
   run() const override
   {
-    IntSequence bnvs(2); bnvs[0] = 13; bnvs[1] = 11;
+    IntSequence bnvs{13, 11};
     return dense_prod(Symmetry{3, 2}, bnvs, 3, 20, 20);
   }
 };
@@ -1116,86 +1106,83 @@ public:
 int
 main()
 {
-  TestRunnable *all_tests[50];
+  std::vector<std::unique_ptr<TestRunnable>> all_tests;
   // fill in vector of all tests
-  int num_tests = 0;
-  all_tests[num_tests++] = new SmallIndexForwardFold();
-  all_tests[num_tests++] = new SmallIndexForwardUnfold();
-  all_tests[num_tests++] = new IndexForwardFold();
-  all_tests[num_tests++] = new IndexForwardUnfold();
-  all_tests[num_tests++] = new SmallIndexBackwardFold();
-  all_tests[num_tests++] = new IndexBackwardFold();
-  all_tests[num_tests++] = new SmallIndexBackwardUnfold();
-  all_tests[num_tests++] = new IndexBackwardUnfold();
-  all_tests[num_tests++] = new SmallIndexOffsetFold();
-  all_tests[num_tests++] = new SmallIndexOffsetUnfold();
-  all_tests[num_tests++] = new IndexOffsetFold();
-  all_tests[num_tests++] = new IndexOffsetUnfold();
-  all_tests[num_tests++] = new SmallFoldUnfoldFS();
-  all_tests[num_tests++] = new SmallFoldUnfoldGS();
-  all_tests[num_tests++] = new FoldUnfoldFS();
-  all_tests[num_tests++] = new FoldUnfoldGS();
-  all_tests[num_tests++] = new SmallFoldUnfoldR();
-  all_tests[num_tests++] = new FoldUnfoldR();
-  all_tests[num_tests++] = new SmallDenseProd();
-  all_tests[num_tests++] = new DenseProd();
-  all_tests[num_tests++] = new BigDenseProd();
-  all_tests[num_tests++] = new SmallFoldedMonomial();
-  all_tests[num_tests++] = new FoldedMonomial();
-  all_tests[num_tests++] = new SmallUnfoldedMonomial();
-  all_tests[num_tests++] = new UnfoldedMonomial();
-  all_tests[num_tests++] = new FoldedContractionSmall();
-  all_tests[num_tests++] = new FoldedContractionBig();
-  all_tests[num_tests++] = new UnfoldedContractionSmall();
-  all_tests[num_tests++] = new UnfoldedContractionBig();
-  all_tests[num_tests++] = new PolyEvalSmall();
-  all_tests[num_tests++] = new PolyEvalBig();
-  all_tests[num_tests++] = new FoldZContSmall();
-  all_tests[num_tests++] = new FoldZCont();
-  all_tests[num_tests++] = new UnfoldZContSmall();
-  all_tests[num_tests++] = new UnfoldZCont();
+  all_tests.push_back(std::make_unique<SmallIndexForwardFold>());
+  all_tests.push_back(std::make_unique<SmallIndexForwardUnfold>());
+  all_tests.push_back(std::make_unique<IndexForwardFold>());
+  all_tests.push_back(std::make_unique<IndexForwardUnfold>());
+  all_tests.push_back(std::make_unique<SmallIndexBackwardFold>());
+  all_tests.push_back(std::make_unique<IndexBackwardFold>());
+  all_tests.push_back(std::make_unique<SmallIndexBackwardUnfold>());
+  all_tests.push_back(std::make_unique<IndexBackwardUnfold>());
+  all_tests.push_back(std::make_unique<SmallIndexOffsetFold>());
+  all_tests.push_back(std::make_unique<SmallIndexOffsetUnfold>());
+  all_tests.push_back(std::make_unique<IndexOffsetFold>());
+  all_tests.push_back(std::make_unique<IndexOffsetUnfold>());
+  all_tests.push_back(std::make_unique<SmallFoldUnfoldFS>());
+  all_tests.push_back(std::make_unique<SmallFoldUnfoldGS>());
+  all_tests.push_back(std::make_unique<FoldUnfoldFS>());
+  all_tests.push_back(std::make_unique<FoldUnfoldGS>());
+  all_tests.push_back(std::make_unique<SmallFoldUnfoldR>());
+  all_tests.push_back(std::make_unique<FoldUnfoldR>());
+  all_tests.push_back(std::make_unique<SmallDenseProd>());
+  all_tests.push_back(std::make_unique<DenseProd>());
+  all_tests.push_back(std::make_unique<BigDenseProd>());
+  all_tests.push_back(std::make_unique<SmallFoldedMonomial>());
+  all_tests.push_back(std::make_unique<FoldedMonomial>());
+  all_tests.push_back(std::make_unique<SmallUnfoldedMonomial>());
+  all_tests.push_back(std::make_unique<UnfoldedMonomial>());
+  all_tests.push_back(std::make_unique<FoldedContractionSmall>());
+  all_tests.push_back(std::make_unique<FoldedContractionBig>());
+  all_tests.push_back(std::make_unique<UnfoldedContractionSmall>());
+  all_tests.push_back(std::make_unique<UnfoldedContractionBig>());
+  all_tests.push_back(std::make_unique<PolyEvalSmall>());
+  all_tests.push_back(std::make_unique<PolyEvalBig>());
+  all_tests.push_back(std::make_unique<FoldZContSmall>());
+  all_tests.push_back(std::make_unique<FoldZCont>());
+  all_tests.push_back(std::make_unique<UnfoldZContSmall>());
+  all_tests.push_back(std::make_unique<UnfoldZCont>());
 
   // find maximum dimension and maximum nvar
   int dmax = 0;
   int nvmax = 0;
-  for (int i = 0; i < num_tests; i++)
+  for (const auto &test : all_tests)
     {
-      if (dmax < all_tests[i]->dim)
-        dmax = all_tests[i]->dim;
-      if (nvmax < all_tests[i]->nvar)
-        nvmax = all_tests[i]->nvar;
+      if (dmax < test->dim)
+        dmax = test->dim;
+      if (nvmax < test->nvar)
+        nvmax = test->nvar;
     }
   tls.init(dmax, nvmax); // initialize library
 
   // launch the tests
   int success = 0;
-  for (int i = 0; i < num_tests; i++)
+  for (const auto &test : all_tests)
     {
       try
         {
-          if (all_tests[i]->test())
+          if (test->test())
             success++;
         }
       catch (const TLException &e)
         {
-          printf("Caugth TL exception in <%s>:\n", all_tests[i]->getName());
+          std::cout << "Caught TL exception in <" << test->name << ">:\n";
           e.print();
         }
       catch (SylvException &e)
         {
-          printf("Caught Sylv exception in <%s>:\n", all_tests[i]->getName());
+          std::cout << "Caught Sylv exception in <" << test->name << ">:\n";
           e.printMessage();
         }
     }
 
-  printf("There were %d tests that failed out of %d tests run.\n",
-         num_tests - success, num_tests);
+  int nfailed = all_tests.size() - success;
+  std::cout << "There were " << nfailed << " tests that failed out of "
+            << all_tests.size() << " tests run." << std::endl;
 
-  // destroy
-  for (int i = 0; i < num_tests; i++)
-    {
-      delete all_tests[i];
-    }
-
-  return 0;
+  if (nfailed)
+    return EXIT_FAILURE;
+  else
+    return EXIT_SUCCESS;
 }