diff --git a/cmake/PugsDoc.cmake b/cmake/PugsDoc.cmake
index c4af5758153ab8c75e0ad4bb82a51af8d89090eb..099eda180d7b37d3824fbfac370151c018718952 100644
--- a/cmake/PugsDoc.cmake
+++ b/cmake/PugsDoc.cmake
@@ -57,7 +57,7 @@ if (EMACS)
     pugsdoc-download-elpa
     ${ORG_GENERATOR_FILES}
     WORKING_DIRECTORY ${PUGS_BINARY_DIR}/doc
-    COMMENT "Built user documentation in doc/userdoc.html"
+    COMMENT "Building user documentation in doc/userdoc.html"
     VERBATIM)
 
   add_custom_target(userdoc-html DEPENDS pugsdoc-dir "${PUGS_BINARY_DIR}/doc/userdoc.html")
@@ -83,7 +83,7 @@ if (EMACS)
       pugsdoc-download-elpa
       ${ORG_GENERATOR_FILES}
       WORKING_DIRECTORY ${PUGS_BINARY_DIR}/doc
-      COMMENT "Built user documentation in doc/userdoc.pdf"
+      COMMENT "Building user documentation in doc/userdoc.pdf"
       VERBATIM)
 
     add_custom_target(userdoc-pdf DEPENDS pugsdoc-dir "${PUGS_BINARY_DIR}/doc/userdoc.pdf" )
diff --git a/doc/userdoc.org b/doc/userdoc.org
index 730486dd17ed76e0a30e2a935807ea36228317d2..fdb22f61fad976e3cda8e0a32448cc5c716dfa0a 100644
--- a/doc/userdoc.org
+++ b/doc/userdoc.org
@@ -661,6 +661,20 @@ expressions.
   123e-2;
   12.3E-1;
 #+END_SRC
+- Small vectors values (types ~R^1~, ~R^2~ and ~R^3~) are written through
+  brackets. Each component of vector values must be a scalar.
+#+BEGIN_SRC pugs :exports both
+  [1];     // R^1 value
+  [1,2];   // R^2 value
+  [1,2,3]; // R^3 value
+#+END_SRC
+- Small matrices values (types ~R^1x1~, ~R^2x2~ and ~R^3x3~) are written by
+  lines through brackets. Each line is enclosed between brackets.
+#+BEGIN_SRC pugs :exports both
+  [[1]];                     // R^1x1 value
+  [[1,2],[3,4]];             // R^2x2 value
+  [[1,2,3],[4,5,6],[7,8,9]]; // R^3x3 value
+#+END_SRC
 - ~string~ values are defined as the set of characters enclosed between
   two double quotes ( ~"~ ). The string /Hello world!/ would be simply
   written as ~"Hello world!"~. Strings support the following escape
@@ -1666,12 +1680,6 @@ the output is
 This is completely equivalent to declaring the variables one after the
 other.
 
-#+BEGIN_note
-Observe that there is no implicit conversion when dealing with
-compound types. The ~=~ operators are used sequentially to set the
-different variables.
-#+END_note
-
 **** Compound definition
 
 One can also use the following definition instruction
@@ -1712,6 +1720,12 @@ to change (allow to use the fresh value of ~x~ in the definition of ~y~),
 but this make the code unclear and this is not the purpose of compound
 types.
 
+#+BEGIN_note
+Observe that there is no implicit conversion when dealing with
+compound types. The ~=~ operators are used sequentially to set the
+different variables.
+#+END_note
+
 **** Compound affectation
 
 The last way to use compound types is illustrated by the following
@@ -1934,7 +1948,7 @@ One gets
 
 *** ~do...while~ loops
 
-The second kind of loops that can be supported is the ~do...while~
+The second kind of loops that is supported is the ~do...while~
 construction which executes /at least/ one time the enclosed ~statement~.
 #+BEGIN_SRC pugs :exports code
   do statement while (condition);
diff --git a/src/analysis/CubeGaussQuadrature.cpp b/src/analysis/CubeGaussQuadrature.cpp
index 8f49fccbeb62b5e4e06da65fbccc0050bb0a0ce8..d7b7afff49b643787eb3e111903885bd324c9fa1 100644
--- a/src/analysis/CubeGaussQuadrature.cpp
+++ b/src/analysis/CubeGaussQuadrature.cpp
@@ -1,5 +1,7 @@
 #include <analysis/CubeGaussQuadrature.hpp>
+
 #include <utils/Exceptions.hpp>
+#include <utils/Stringify.hpp>
 
 void
 CubeGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
@@ -495,7 +497,7 @@ CubeGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
     // LCOV_EXCL_START
   default: {
     throw UnexpectedError("Gauss quadrature formulae handle degrees up to " +
-                          std::to_string(CubeGaussQuadrature::max_degree) + " on cubes");
+                          stringify(CubeGaussQuadrature::max_degree) + " on cubes");
   }
     // LCOV_EXCL_STOP
   }
diff --git a/src/analysis/PrismGaussQuadrature.cpp b/src/analysis/PrismGaussQuadrature.cpp
index 623d3f87864b08fc7b65586ebf09453568dabcf1..e1ae88d38623f1f63e43793e8be914e8ea42e3dd 100644
--- a/src/analysis/PrismGaussQuadrature.cpp
+++ b/src/analysis/PrismGaussQuadrature.cpp
@@ -1,5 +1,7 @@
 #include <analysis/PrismGaussQuadrature.hpp>
+
 #include <utils/Exceptions.hpp>
+#include <utils/Stringify.hpp>
 
 void
 PrismGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
@@ -799,7 +801,7 @@ PrismGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
     // LCOV_EXCL_START
   default: {
     throw UnexpectedError("Gauss quadrature formulae handle degrees up to " +
-                          std::to_string(PrismGaussQuadrature::max_degree) + "on prisms");
+                          stringify(PrismGaussQuadrature::max_degree) + "on prisms");
   }
     // LCOV_EXCL_STOP
   }
diff --git a/src/analysis/PyramidGaussQuadrature.cpp b/src/analysis/PyramidGaussQuadrature.cpp
index 48eb4ec3421f1f6498420c9c851994520dc1b486..679f3ded5b85764bca43ed062f379d8171cd5c65 100644
--- a/src/analysis/PyramidGaussQuadrature.cpp
+++ b/src/analysis/PyramidGaussQuadrature.cpp
@@ -1,5 +1,7 @@
 #include <analysis/PyramidGaussQuadrature.hpp>
+
 #include <utils/Exceptions.hpp>
+#include <utils/Stringify.hpp>
 
 void
 PyramidGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
@@ -971,7 +973,7 @@ PyramidGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
     // LCOV_EXCL_START
   default: {
     throw UnexpectedError("Gauss quadrature formulae handle degrees up to " +
-                          std::to_string(PyramidGaussQuadrature::max_degree) + " on pyramids");
+                          stringify(PyramidGaussQuadrature::max_degree) + " on pyramids");
   }
     // LCOV_EXCL_STOP
   }
diff --git a/src/analysis/QuadratureManager.hpp b/src/analysis/QuadratureManager.hpp
index 04655f291f8ff75e8c1207c4b6171a7b7f461a05..36ea0db4818661d7554d3579bad9c54004fe11ca 100644
--- a/src/analysis/QuadratureManager.hpp
+++ b/src/analysis/QuadratureManager.hpp
@@ -7,6 +7,7 @@
 #include <utils/Exceptions.hpp>
 #include <utils/PugsAssert.hpp>
 #include <utils/PugsMacros.hpp>
+#include <utils/Stringify.hpp>
 
 class QuadratureManager
 {
@@ -49,7 +50,7 @@ class QuadratureManager
   {
     if (quadrature_descriptor.degree() > maxLineDegree(quadrature_descriptor.type())) {
       throw NormalError(::name(quadrature_descriptor.type()) + " quadrature formulae handle degrees up to " +
-                        std::to_string(maxLineDegree(quadrature_descriptor.type())) + " on lines");
+                        stringify(maxLineDegree(quadrature_descriptor.type())) + " on lines");
     }
     switch (quadrature_descriptor.type()) {
     case QuadratureType::Gauss:
@@ -73,7 +74,7 @@ class QuadratureManager
   {
     if (quadrature_descriptor.degree() > maxTriangleDegree(quadrature_descriptor.type())) {
       throw NormalError(::name(quadrature_descriptor.type()) + " quadrature formulae handle degrees up to " +
-                        std::to_string(maxTriangleDegree(quadrature_descriptor.type())) + " on triangles");
+                        stringify(maxTriangleDegree(quadrature_descriptor.type())) + " on triangles");
     }
     switch (quadrature_descriptor.type()) {
     case QuadratureType::Gauss: {
@@ -92,7 +93,7 @@ class QuadratureManager
   {
     if (quadrature_descriptor.degree() > maxSquareDegree(quadrature_descriptor.type())) {
       throw NormalError(::name(quadrature_descriptor.type()) + " quadrature formulae handle degrees up to " +
-                        std::to_string(maxSquareDegree(quadrature_descriptor.type())) + " on squares");
+                        stringify(maxSquareDegree(quadrature_descriptor.type())) + " on squares");
     }
     switch (quadrature_descriptor.type()) {
     case QuadratureType::Gauss: {
@@ -117,7 +118,7 @@ class QuadratureManager
   {
     if (quadrature_descriptor.degree() > maxTetrahedronDegree(quadrature_descriptor.type())) {
       throw NormalError(::name(quadrature_descriptor.type()) + " quadrature formulae handle degrees up to " +
-                        std::to_string(maxTetrahedronDegree(quadrature_descriptor.type())) + " on tetrahedra");
+                        stringify(maxTetrahedronDegree(quadrature_descriptor.type())) + " on tetrahedra");
     }
     switch (quadrature_descriptor.type()) {
     case QuadratureType::Gauss: {
@@ -136,7 +137,7 @@ class QuadratureManager
   {
     if (quadrature_descriptor.degree() > maxPrismDegree(quadrature_descriptor.type())) {
       throw NormalError(::name(quadrature_descriptor.type()) + " quadrature formulae handle degrees up to " +
-                        std::to_string(maxPrismDegree(quadrature_descriptor.type())) + " on prisms");
+                        stringify(maxPrismDegree(quadrature_descriptor.type())) + " on prisms");
     }
     switch (quadrature_descriptor.type()) {
     case QuadratureType::Gauss: {
@@ -155,7 +156,7 @@ class QuadratureManager
   {
     if (quadrature_descriptor.degree() > maxPyramidDegree(quadrature_descriptor.type())) {
       throw NormalError(::name(quadrature_descriptor.type()) + " quadrature formulae handle degrees up to " +
-                        std::to_string(maxPyramidDegree(quadrature_descriptor.type())) + " on pyramids");
+                        stringify(maxPyramidDegree(quadrature_descriptor.type())) + " on pyramids");
     }
     switch (quadrature_descriptor.type()) {
     case QuadratureType::Gauss: {
@@ -174,7 +175,7 @@ class QuadratureManager
   {
     if (quadrature_descriptor.degree() > maxCubeDegree(quadrature_descriptor.type())) {
       throw NormalError(::name(quadrature_descriptor.type()) + " quadrature formulae handle degrees up to " +
-                        std::to_string(maxCubeDegree(quadrature_descriptor.type())) + " on cubes");
+                        stringify(maxCubeDegree(quadrature_descriptor.type())) + " on cubes");
     }
     switch (quadrature_descriptor.type()) {
     case QuadratureType::Gauss: {
diff --git a/src/analysis/SquareGaussQuadrature.cpp b/src/analysis/SquareGaussQuadrature.cpp
index 2597e013b1b4f7ede3b69981fb43bfd4e901aff7..07ad592bcbc205620e1d45a8ade14c1467caf153 100644
--- a/src/analysis/SquareGaussQuadrature.cpp
+++ b/src/analysis/SquareGaussQuadrature.cpp
@@ -1,5 +1,7 @@
 #include <analysis/SquareGaussQuadrature.hpp>
+
 #include <utils/Exceptions.hpp>
+#include <utils/Stringify.hpp>
 
 void
 SquareGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
@@ -322,7 +324,7 @@ SquareGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
     // LCOV_EXCL_START
   default: {
     throw UnexpectedError("Gauss quadrature formulae handle degrees up to " +
-                          std::to_string(SquareGaussQuadrature::max_degree) + " on squares");
+                          stringify(SquareGaussQuadrature::max_degree) + " on squares");
   }
     // LCOV_EXCL_STOP
   }
diff --git a/src/analysis/TensorialGaussLegendreQuadrature.cpp b/src/analysis/TensorialGaussLegendreQuadrature.cpp
index a306f44135b7a4cf31ddc6e906f0352ca1bde067..7e58462c997236039ec08e90de890b2160013f47 100644
--- a/src/analysis/TensorialGaussLegendreQuadrature.cpp
+++ b/src/analysis/TensorialGaussLegendreQuadrature.cpp
@@ -1,5 +1,6 @@
 #include <analysis/TensorialGaussLegendreQuadrature.hpp>
 #include <utils/Exceptions.hpp>
+#include <utils/Stringify.hpp>
 
 template <>
 void
@@ -232,7 +233,7 @@ TensorialGaussLegendreQuadrature<1>::_buildPointAndWeightLists(const size_t degr
     // LCOV_EXCL_START
   default: {
     throw UnexpectedError("Gauss-Legendre quadrature formulae handle degrees up to " +
-                          std::to_string(TensorialGaussLegendreQuadrature<1>::max_degree));
+                          stringify(TensorialGaussLegendreQuadrature<1>::max_degree));
   }
     // LCOV_EXCL_STOP
   }
diff --git a/src/analysis/TensorialGaussLobattoQuadrature.cpp b/src/analysis/TensorialGaussLobattoQuadrature.cpp
index f53b49434335643684be80e20d5d750dde1c4e93..45bd293df34e30ede40920e36e1c601dca5732fd 100644
--- a/src/analysis/TensorialGaussLobattoQuadrature.cpp
+++ b/src/analysis/TensorialGaussLobattoQuadrature.cpp
@@ -1,5 +1,6 @@
 #include <analysis/TensorialGaussLobattoQuadrature.hpp>
 #include <utils/Exceptions.hpp>
+#include <utils/Stringify.hpp>
 
 template <>
 void
@@ -121,7 +122,7 @@ TensorialGaussLobattoQuadrature<1>::_buildPointAndWeightLists(const size_t degre
     // LCOV_EXCL_START
   default: {
     throw UnexpectedError("Gauss-Lobatto quadrature formulae handle degrees up to " +
-                          std::to_string(TensorialGaussLobattoQuadrature<1>::max_degree));
+                          stringify(TensorialGaussLobattoQuadrature<1>::max_degree));
   }
     // LCOV_EXCL_STOP
   }
diff --git a/src/analysis/TetrahedronGaussQuadrature.cpp b/src/analysis/TetrahedronGaussQuadrature.cpp
index 49fd5efe685c765aedccae3ee1b2c59d82915dd3..b6f6bc795fbf997720c7c479459362f6b798957a 100644
--- a/src/analysis/TetrahedronGaussQuadrature.cpp
+++ b/src/analysis/TetrahedronGaussQuadrature.cpp
@@ -1,5 +1,7 @@
 #include <analysis/TetrahedronGaussQuadrature.hpp>
+
 #include <utils/Exceptions.hpp>
+#include <utils/Stringify.hpp>
 
 void
 TetrahedronGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
@@ -683,7 +685,7 @@ TetrahedronGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
     // LCOV_EXCL_START
   default: {
     throw UnexpectedError("Gauss quadrature formulae handle degrees up to " +
-                          std::to_string(TetrahedronGaussQuadrature::max_degree) + " on tetrahedra");
+                          stringify(TetrahedronGaussQuadrature::max_degree) + " on tetrahedra");
   }
     // LCOV_EXCL_STOP
   }
diff --git a/src/analysis/TriangleGaussQuadrature.cpp b/src/analysis/TriangleGaussQuadrature.cpp
index 627edc4ff8293a7a00b689b5602fa2f8d20c3911..c00392d02bd1cf851f5ca8dd73eac7c87e39ddd4 100644
--- a/src/analysis/TriangleGaussQuadrature.cpp
+++ b/src/analysis/TriangleGaussQuadrature.cpp
@@ -1,5 +1,7 @@
 #include <analysis/TriangleGaussQuadrature.hpp>
+
 #include <utils/Exceptions.hpp>
+#include <utils/Stringify.hpp>
 
 void
 TriangleGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
@@ -485,7 +487,7 @@ TriangleGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
     // LCOV_EXCL_START
   default: {
     throw UnexpectedError("Gauss quadrature formulae handle degrees up to " +
-                          std::to_string(TriangleGaussQuadrature::max_degree) + " on triangles");
+                          stringify(TriangleGaussQuadrature::max_degree) + " on triangles");
   }
     // LCOV_EXCL_STOP
   }
diff --git a/src/language/modules/MeshModule.cpp b/src/language/modules/MeshModule.cpp
index b0f1a2313a6db1a5164ea46862d4fd84d0ddd6a1..f10c4a7f30cdb6385d86a7c5e858a5cdd89cf576 100644
--- a/src/language/modules/MeshModule.cpp
+++ b/src/language/modules/MeshModule.cpp
@@ -114,8 +114,8 @@ MeshModule::MeshModule()
                                 constexpr uint64_t dimension = 1;
 
                                 if (box_sizes.size() != dimension) {
-                                  throw NormalError("expecting " + std::to_string(dimension) +
-                                                    " dimensions, provided " + std::to_string(box_sizes.size()));
+                                  throw NormalError("expecting " + stringify(dimension) + " dimensions, provided " +
+                                                    stringify(box_sizes.size()));
                                 }
 
                                 const TinyVector<dimension, uint64_t> sizes = [&]() {
@@ -141,8 +141,8 @@ MeshModule::MeshModule()
                                 constexpr uint64_t dimension = 2;
 
                                 if (box_sizes.size() != dimension) {
-                                  throw NormalError("expecting " + std::to_string(dimension) +
-                                                    " dimensions, provided " + std::to_string(box_sizes.size()));
+                                  throw NormalError("expecting " + stringify(dimension) + " dimensions, provided " +
+                                                    stringify(box_sizes.size()));
                                 }
 
                                 const TinyVector<dimension, uint64_t> sizes = [&]() {
@@ -168,8 +168,8 @@ MeshModule::MeshModule()
                                 constexpr uint64_t dimension = 3;
 
                                 if (box_sizes.size() != dimension) {
-                                  throw NormalError("expecting " + std::to_string(dimension) +
-                                                    " dimensions, provided " + std::to_string(box_sizes.size()));
+                                  throw NormalError("expecting " + stringify(dimension) + " dimensions, provided " +
+                                                    stringify(box_sizes.size()));
                                 }
 
                                 const TinyVector<dimension, uint64_t> sizes = [&]() {
diff --git a/src/language/node_processor/AffectationProcessor.hpp b/src/language/node_processor/AffectationProcessor.hpp
index a318413fe7ff7b06e820c25b1c5ef268ac3c5c40..f4c0d731a7a11ab91ffe7cb005bf14ebb5525dcf 100644
--- a/src/language/node_processor/AffectationProcessor.hpp
+++ b/src/language/node_processor/AffectationProcessor.hpp
@@ -7,6 +7,7 @@
 #include <language/utils/SymbolTable.hpp>
 #include <utils/Exceptions.hpp>
 #include <utils/PugsTraits.hpp>
+#include <utils/Stringify.hpp>
 
 #include <exception>
 
@@ -16,22 +17,13 @@ struct AffOp;
 template <>
 struct AffOp<language::multiplyeq_op>
 {
-  template <typename T>
-  std::string
-  _stringify(const T& value)
-  {
-    std::ostringstream os;
-    os << std::boolalpha << value;
-    return os.str();
-  }
-
   template <typename A, typename B>
   PUGS_INLINE void
   eval(A& a, const B& b)
   {
     if constexpr (std::is_same_v<uint64_t, A> and std::is_same_v<int64_t, B>) {
       if (b < 0) {
-        throw std::domain_error("trying to affect negative value (" + _stringify(b) + ")");
+        throw std::domain_error("trying to affect negative value (" + stringify(b) + ")");
       }
     }
     a *= b;
@@ -41,22 +33,13 @@ struct AffOp<language::multiplyeq_op>
 template <>
 struct AffOp<language::divideeq_op>
 {
-  template <typename T>
-  std::string
-  _stringify(const T& value)
-  {
-    std::ostringstream os;
-    os << std::boolalpha << value;
-    return os.str();
-  }
-
   template <typename A, typename B>
   PUGS_INLINE void
   eval(A& a, const B& b)
   {
     if constexpr (std::is_same_v<uint64_t, A> and std::is_same_v<int64_t, B>) {
       if (b < 0) {
-        throw std::domain_error("trying to affect negative value (" + _stringify(b) + ")");
+        throw std::domain_error("trying to affect negative value (" + stringify(b) + ")");
       }
     }
     a /= b;
@@ -66,22 +49,13 @@ struct AffOp<language::divideeq_op>
 template <>
 struct AffOp<language::pluseq_op>
 {
-  template <typename T>
-  std::string
-  _stringify(const T& value)
-  {
-    std::ostringstream os;
-    os << std::boolalpha << value;
-    return os.str();
-  }
-
   template <typename A, typename B>
   PUGS_INLINE void
   eval(A& a, const B& b)
   {
     if constexpr (std::is_same_v<uint64_t, A> and std::is_same_v<int64_t, B>) {
       if (static_cast<int64_t>(a + b) < 0) {
-        throw std::domain_error("trying to affect negative value (lhs: " + _stringify(a) + " rhs: " + _stringify(b) +
+        throw std::domain_error("trying to affect negative value (lhs: " + stringify(a) + " rhs: " + stringify(b) +
                                 ")");
       }
     }
@@ -92,22 +66,13 @@ struct AffOp<language::pluseq_op>
 template <>
 struct AffOp<language::minuseq_op>
 {
-  template <typename T>
-  std::string
-  _stringify(const T& value)
-  {
-    std::ostringstream os;
-    os << std::boolalpha << value;
-    return os.str();
-  }
-
   template <typename A, typename B>
   PUGS_INLINE void
   eval(A& a, const B& b)
   {
     if constexpr (std::is_same_v<uint64_t, A> and std::is_same_v<int64_t, B>) {
       if (static_cast<int64_t>(a - b) < 0) {
-        throw std::domain_error("trying to affect negative value (lhs: " + _stringify(a) + " rhs: " + _stringify(b) +
+        throw std::domain_error("trying to affect negative value (lhs: " + stringify(a) + " rhs: " + stringify(b) +
                                 ")");
       }
     }
@@ -143,15 +108,6 @@ class AffectationExecutor final : public IAffectationExecutor
     return true;
   }()};
 
-  template <typename T>
-  std::string
-  _stringify(const T& value)
-  {
-    std::ostringstream os;
-    os << std::boolalpha << value;
-    return os.str();
-  }
-
  public:
   AffectationExecutor(ASTNode& node, DataVariant& lhs) : m_lhs(lhs), m_node{node}
   {
@@ -172,14 +128,14 @@ class AffectationExecutor final : public IAffectationExecutor
             if constexpr (std::is_same_v<std::string, DataT>) {
               m_lhs = std::get<DataT>(rhs);
             } else {
-              m_lhs = std::move(_stringify(std::get<DataT>(rhs)));
+              m_lhs = std::move(stringify(std::get<DataT>(rhs)));
             }
           } else {
             ValueT& lhs = std::get<ValueT>(m_lhs);
             if constexpr (std::is_same_v<std::string, DataT>) {
               lhs += std::get<std::string>(rhs);
             } else {
-              lhs += std::move(_stringify(std::get<DataT>(rhs)));
+              lhs += std::move(stringify(std::get<DataT>(rhs)));
             }
           }
         } else {
@@ -188,7 +144,7 @@ class AffectationExecutor final : public IAffectationExecutor
               const DataT& value = std::get<DataT>(rhs);
               if constexpr (std::is_same_v<uint64_t, ValueT> and std::is_same_v<int64_t, DataT>) {
                 if (value < 0) {
-                  throw std::domain_error("trying to affect negative value (" + _stringify(value) + ")");
+                  throw std::domain_error("trying to affect negative value (" + stringify(value) + ")");
                 }
               }
               m_lhs = static_cast<ValueT>(value);
@@ -236,7 +192,7 @@ class AffectationExecutor final : public IAffectationExecutor
                           if constexpr (std::is_same_v<ValueContentT, uint64_t> and
                                         std::is_same_v<DataContentT, int64_t>) {
                             if (vi < 0) {
-                              throw std::domain_error("trying to affect negative value (" + _stringify(v) + ")");
+                              throw std::domain_error("trying to affect negative value (" + stringify(v) + ")");
                             }
                           }
                           lhs[i] = vi;
@@ -266,9 +222,9 @@ class AffectationExecutor final : public IAffectationExecutor
                         if constexpr (std::is_same_v<typename V_T::value_type, bool>) {
                           // Ugly workaround to allow compilation with libstdc++-9
                           bool v_i = v[i];
-                          lhs[i]   = std::move(_stringify(v_i));
+                          lhs[i]   = std::move(stringify(v_i));
                         } else {
-                          lhs[i] = std::move(_stringify(v[i]));
+                          lhs[i] = std::move(stringify(v[i]));
                         }
                       }
                     } else {
@@ -355,12 +311,12 @@ class AffectationExecutor final : public IAffectationExecutor
                                          std::is_convertible_v<T, ValueContentT>) {
                       if constexpr (std::is_same_v<ValueContentT, uint64_t> and std::is_same_v<T, int64_t>) {
                         if (child_value < 0) {
-                          throw std::domain_error("trying to affect negative value (" + _stringify(child_value) + ")");
+                          throw std::domain_error("trying to affect negative value (" + stringify(child_value) + ")");
                         }
                       }
                       lhs[i] = static_cast<ValueContentT>(child_value);
                     } else if constexpr (std::is_same_v<std::string, ValueContentT>) {
-                      lhs[i] = std::move(_stringify(child_value));
+                      lhs[i] = std::move(stringify(child_value));
                     } else if constexpr (is_tiny_vector_v<ValueContentT>) {
                       if constexpr (std::is_arithmetic_v<T>) {
                         if constexpr (std::is_same_v<ValueContentT, TinyVector<1>>) {
@@ -411,12 +367,12 @@ class AffectationExecutor final : public IAffectationExecutor
                                        std::is_convertible_v<T, ValueContentT>) {
                     if constexpr (std::is_same_v<ValueContentT, uint64_t> and std::is_same_v<T, int64_t>) {
                       if (child_value < 0) {
-                        throw std::domain_error("trying to affect negative value (" + _stringify(child_value) + ")");
+                        throw std::domain_error("trying to affect negative value (" + stringify(child_value) + ")");
                       }
                     }
                     lhs[0] = static_cast<ValueContentT>(child_value);
                   } else if constexpr (std::is_same_v<std::string, ValueContentT>) {
-                    lhs[0] = std::move(_stringify(child_value));
+                    lhs[0] = std::move(stringify(child_value));
                   } else if constexpr (is_tiny_vector_v<ValueContentT>) {
                     if constexpr (std::is_arithmetic_v<T>) {
                       if constexpr (std::is_same_v<ValueContentT, TinyVector<1>>) {
@@ -540,15 +496,6 @@ class AffectationToTupleProcessor final : public AffectationToDataVariantProcess
  private:
   ASTNode& m_rhs_node;
 
-  template <typename T>
-  std::string
-  _stringify(const T& value)
-  {
-    std::ostringstream os;
-    os << std::boolalpha << value;
-    return os.str();
-  }
-
  public:
   DataVariant
   execute(ExecutionPolicy& exec_policy)
@@ -564,12 +511,12 @@ class AffectationToTupleProcessor final : public AffectationToDataVariantProcess
           } else if constexpr (std::is_arithmetic_v<ValueT> and std::is_convertible_v<T, ValueT>) {
             if constexpr (std::is_same_v<uint64_t, ValueT> and std::is_same_v<int64_t, T>) {
               if (v < 0) {
-                throw std::domain_error("trying to affect negative value (" + _stringify(v) + ")");
+                throw std::domain_error("trying to affect negative value (" + stringify(v) + ")");
               }
             }
             *m_lhs = std::vector{std::move(static_cast<ValueT>(v))};
           } else if constexpr (std::is_same_v<std::string, ValueT>) {
-            *m_lhs = std::vector{std::move(_stringify(v))};
+            *m_lhs = std::vector{std::move(stringify(v))};
           } else if constexpr (is_tiny_vector_v<ValueT> or is_tiny_matrix_v<ValueT>) {
             if constexpr (std::is_same_v<ValueT, TinyVector<1>> and std::is_arithmetic_v<T>) {
               *m_lhs = std::vector<TinyVector<1>>{TinyVector<1>{static_cast<double>(v)}};
@@ -610,7 +557,7 @@ class AffectationToTupleFromListProcessor final : public AffectationToDataVarian
 
   template <typename T>
   std::string
-  _stringify(const T& value)
+  stringify(const T& value)
   {
     std::ostringstream os;
     os << std::boolalpha << value;
@@ -631,12 +578,12 @@ class AffectationToTupleFromListProcessor final : public AffectationToDataVarian
             } else if constexpr (std::is_arithmetic_v<ValueT> and std::is_convertible_v<T, ValueT>) {
               if constexpr (std::is_same_v<uint64_t, ValueT> and std::is_same_v<int64_t, T>) {
                 if (child_value < 0) {
-                  throw std::domain_error("trying to affect negative value (" + _stringify(child_value) + ")");
+                  throw std::domain_error("trying to affect negative value (" + stringify(child_value) + ")");
                 }
               }
               tuple_value[i] = static_cast<ValueT>(child_value);
             } else if constexpr (std::is_same_v<std::string, ValueT>) {
-              tuple_value[i] = std::move(_stringify(child_value));
+              tuple_value[i] = std::move(stringify(child_value));
             } else if constexpr (is_tiny_vector_v<ValueT>) {
               if constexpr (std::is_arithmetic_v<T>) {
                 if constexpr (std::is_same_v<ValueT, TinyVector<1>>) {
@@ -694,14 +641,14 @@ class AffectationToTupleFromListProcessor final : public AffectationToDataVarian
         const DataType& vi = values[i];
         if constexpr (std::is_same_v<uint64_t, ValueT> and std::is_same_v<int64_t, DataType>) {
           if (vi < 0) {
-            throw std::domain_error("trying to affect negative value (" + _stringify(vi) + ")");
+            throw std::domain_error("trying to affect negative value (" + stringify(vi) + ")");
           }
         }
         v[i] = static_cast<DataType>(vi);
       }
     } else if constexpr (std::is_same_v<ValueT, std::string>) {
       for (size_t i = 0; i < values.size(); ++i) {
-        v[i] = std::move(_stringify(values[i]));
+        v[i] = std::move(stringify(values[i]));
       }
     } else {
       // LCOV_EXCL_START
diff --git a/src/language/node_processor/ConcatExpressionProcessor.hpp b/src/language/node_processor/ConcatExpressionProcessor.hpp
index 43b785979c986dd0d2525ace5c55f69c30faa259..4653c341c7b322fb338e15c09dada90b98472a8c 100644
--- a/src/language/node_processor/ConcatExpressionProcessor.hpp
+++ b/src/language/node_processor/ConcatExpressionProcessor.hpp
@@ -3,6 +3,7 @@
 
 #include <language/ast/ASTNode.hpp>
 #include <language/node_processor/INodeProcessor.hpp>
+#include <utils/Stringify.hpp>
 
 template <typename A_DataT, typename B_DataT>
 class ConcatExpressionProcessor final : public INodeProcessor
@@ -17,7 +18,7 @@ class ConcatExpressionProcessor final : public INodeProcessor
     if constexpr (std::is_same_v<B_DataT, std::string>) {
       return a + std::get<B_DataT>(b);
     } else if constexpr ((std::is_arithmetic_v<B_DataT>)and(not std::is_same_v<B_DataT, bool>)) {
-      return a + std::to_string(std::get<B_DataT>(b));
+      return a + stringify(std::get<B_DataT>(b));
     } else {
       std::ostringstream os;
       os << a << std::boolalpha << b;
@@ -31,7 +32,7 @@ class ConcatExpressionProcessor final : public INodeProcessor
   {
     static_assert(not std::is_same_v<A_DataT, std::string>, "this case is treated by the other eval function");
     if constexpr ((std::is_arithmetic_v<A_DataT>)and(not std::is_same_v<A_DataT, bool>)) {
-      return std::to_string(std::get<A_DataT>(a)) + b;
+      return stringify(std::get<A_DataT>(a)) + b;
     } else {
       std::ostringstream os;
       os << std::boolalpha << a << b;
diff --git a/src/language/node_processor/FunctionArgumentConverter.hpp b/src/language/node_processor/FunctionArgumentConverter.hpp
index 769bcd9ef92d8379064e9d56d491942e48501f24..bb125b821361989e6e8bb898ba0b56208c32ef84 100644
--- a/src/language/node_processor/FunctionArgumentConverter.hpp
+++ b/src/language/node_processor/FunctionArgumentConverter.hpp
@@ -6,6 +6,7 @@
 #include <utils/Demangle.hpp>
 #include <utils/Exceptions.hpp>
 #include <utils/PugsTraits.hpp>
+#include <utils/Stringify.hpp>
 
 #include <sstream>
 
@@ -35,7 +36,7 @@ class FunctionArgumentToStringConverter final : public IFunctionArgumentConverte
       [&](auto&& v) {
         using T = std::decay_t<decltype(v)>;
         if constexpr (std::is_arithmetic_v<T>) {
-          exec_policy.currentContext()[m_argument_id] = std::to_string(v);
+          exec_policy.currentContext()[m_argument_id] = stringify(v);
         } else if constexpr (std::is_same_v<T, std::string>) {
           exec_policy.currentContext()[m_argument_id] = v;
         } else {
@@ -68,7 +69,7 @@ class FunctionArgumentConverter final : public IFunctionArgumentConverter
       const ProvidedValueType& v = std::get<ProvidedValueType>(value);
       if constexpr (std::is_same_v<uint64_t, ExpectedValueType> and std::is_same_v<int64_t, ProvidedValueType>) {
         if (v < 0) {
-          throw std::domain_error("trying to convert negative value (" + std::to_string(v) + ")");
+          throw std::domain_error("trying to convert negative value (" + stringify(v) + ")");
         }
       }
       exec_policy.currentContext()[m_argument_id] = std::move(static_cast<ExpectedValueType>(v));
@@ -208,7 +209,7 @@ class FunctionTupleArgumentConverter final : public IFunctionArgumentConverter
               for (size_t i = 0; i < v.size(); ++i) {
                 if constexpr (std::is_same_v<ContentType, uint64_t> and std::is_same_v<ContentT, int64_t>) {
                   if (v[i] < 0) {
-                    throw std::domain_error("trying to convert negative value (" + std::to_string(v[i]) + ")");
+                    throw std::domain_error("trying to convert negative value (" + stringify(v[i]) + ")");
                   }
                 }
                 list_value.push_back(static_cast<ContentType>(v[i]));
@@ -224,7 +225,7 @@ class FunctionTupleArgumentConverter final : public IFunctionArgumentConverter
                                not is_tiny_matrix_v<ContentType>) {
             if constexpr (std::is_same_v<ContentType, uint64_t> and std::is_same_v<ValueT, int64_t>) {
               if (v < 0) {
-                throw std::domain_error("trying to convert negative value (" + std::to_string(v) + ")");
+                throw std::domain_error("trying to convert negative value (" + stringify(v) + ")");
               }
             }
             exec_policy.currentContext()[m_argument_id] = std::move(TupleType{static_cast<ContentType>(v)});
@@ -322,7 +323,7 @@ class FunctionListArgumentConverter final : public IFunctionArgumentConverter
                   } else if constexpr (std::is_convertible_v<Vi_T, ContentType>) {
                     if constexpr (std::is_same_v<ContentType, uint64_t> and std::is_same_v<Vi_T, int64_t>) {
                       if (vi < 0) {
-                        throw std::domain_error("trying to convert negative value (" + std::to_string(vi) + ")");
+                        throw std::domain_error("trying to convert negative value (" + stringify(vi) + ")");
                       }
                     }
                     list_value.emplace_back(vi);
diff --git a/src/language/node_processor/FunctionProcessor.hpp b/src/language/node_processor/FunctionProcessor.hpp
index 5947cd663cdd8fc738ab8a30f882024613af9a0c..86742e170992cd50cab13ea6c7495f87447a81c1 100644
--- a/src/language/node_processor/FunctionProcessor.hpp
+++ b/src/language/node_processor/FunctionProcessor.hpp
@@ -7,6 +7,7 @@
 #include <language/node_processor/INodeProcessor.hpp>
 #include <language/utils/FunctionTable.hpp>
 #include <language/utils/SymbolTable.hpp>
+#include <utils/Stringify.hpp>
 
 template <typename ReturnType, typename ExpressionValueType>
 class FunctionExpressionProcessor final : public INodeProcessor
@@ -21,14 +22,14 @@ class FunctionExpressionProcessor final : public INodeProcessor
     if constexpr (std::is_same_v<ReturnType, ExpressionValueType>) {
       return m_function_expression.execute(exec_policy);
     } else if constexpr (std::is_same_v<ReturnType, std::string>) {
-      return std::to_string(std::get<ExpressionValueType>(m_function_expression.execute(exec_policy)));
+      return stringify(std::get<ExpressionValueType>(m_function_expression.execute(exec_policy)));
     } else if constexpr (std::is_same_v<ExpressionValueType, ZeroType>) {
       return ReturnType{ZeroType::zero};
     } else if constexpr (std::is_convertible_v<ExpressionValueType, ReturnType>) {
       const ExpressionValueType& v = std::get<ExpressionValueType>(m_function_expression.execute(exec_policy));
       if constexpr (std::is_same_v<ReturnType, uint64_t> and std::is_same_v<ExpressionValueType, int64_t>) {
         if (v < 0) {
-          throw std::domain_error("trying to convert negative value (" + std::to_string(v) + ")");
+          throw std::domain_error("trying to convert negative value (" + stringify(v) + ")");
         }
       }
       return static_cast<ReturnType>(v);
diff --git a/src/language/utils/ASTNodeDataType.cpp b/src/language/utils/ASTNodeDataType.cpp
index 3c6d4af0f8be9f4002834fb15548316f301b0159..3c4420b12a0615ed1e978638a26fd0f2c72009a5 100644
--- a/src/language/utils/ASTNodeDataType.cpp
+++ b/src/language/utils/ASTNodeDataType.cpp
@@ -6,6 +6,7 @@
 #include <language/utils/ParseError.hpp>
 #include <language/utils/SymbolTable.hpp>
 #include <utils/PugsAssert.hpp>
+#include <utils/Stringify.hpp>
 
 ASTNodeDataType
 getVectorDataType(const ASTNode& type_node)
@@ -180,10 +181,10 @@ dataTypeName(const ASTNodeDataType& data_type)
     name = "R";
     break;
   case ASTNodeDataType::vector_t:
-    name = "R^" + std::to_string(data_type.dimension());
+    name = "R^" + stringify(data_type.dimension());
     break;
   case ASTNodeDataType::matrix_t:
-    name = "R^" + std::to_string(data_type.numberOfRows()) + "x" + std::to_string(data_type.numberOfColumns());
+    name = "R^" + stringify(data_type.numberOfRows()) + "x" + stringify(data_type.numberOfColumns());
     break;
   case ASTNodeDataType::tuple_t:
     name = "(" + dataTypeName(data_type.contentType()) + "...)";
diff --git a/src/language/utils/IntegrateOnCells.hpp b/src/language/utils/IntegrateOnCells.hpp
index fa47bb7856ec695d87d7702d9c561a2ebfe6ca6b..79357544738109bc5c5861e01ca22a1ccaf95af8 100644
--- a/src/language/utils/IntegrateOnCells.hpp
+++ b/src/language/utils/IntegrateOnCells.hpp
@@ -224,7 +224,7 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
           } else {
             // LCOV_EXCL_START
             throw NotImplementedError("integration on pyramid with non-quadrangular base (number of nodes " +
-                                      std::to_string(cell_to_node.size()) + ")");
+                                      stringify(cell_to_node.size()) + ")");
             // LCOV_EXCL_STOP
           }
           break;
@@ -287,7 +287,7 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
           } else {
             // LCOV_EXCL_START
             throw NotImplementedError("integration on diamond with non-quadrangular base (" +
-                                      std::to_string(cell_to_node.size()) + ")");
+                                      stringify(cell_to_node.size()) + ")");
             // LCOV_EXCL_STOP
           }
           break;
@@ -532,7 +532,7 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
           } else {
             // LCOV_EXCL_START
             throw NotImplementedError("integration on diamond with non-quadrangular base (" +
-                                      std::to_string(cell_to_node.size()) + ")");
+                                      stringify(cell_to_node.size()) + ")");
             // LCOV_EXCL_STOP
           }
           break;
diff --git a/src/mesh/DiamondDualConnectivityBuilder.cpp b/src/mesh/DiamondDualConnectivityBuilder.cpp
index 1afa7d40413f2f8c9d88a99004a78b094d2705e5..71e316461d94a0891b721a0fff85e85c145ef26f 100644
--- a/src/mesh/DiamondDualConnectivityBuilder.cpp
+++ b/src/mesh/DiamondDualConnectivityBuilder.cpp
@@ -9,6 +9,7 @@
 #include <mesh/RefId.hpp>
 #include <utils/Array.hpp>
 #include <utils/Messenger.hpp>
+#include <utils/Stringify.hpp>
 
 #include <vector>
 
@@ -445,7 +446,7 @@ DiamondDualConnectivityBuilder::DiamondDualConnectivityBuilder(const IConnectivi
   }
     // LCOV_EXCL_START
   default: {
-    throw UnexpectedError("invalid connectivity dimension: " + std::to_string(connectivity.dimension()));
+    throw UnexpectedError("invalid connectivity dimension: " + stringify(connectivity.dimension()));
   }
     // LCOV_EXCL_STOP
   }
diff --git a/src/mesh/DiamondDualMeshBuilder.cpp b/src/mesh/DiamondDualMeshBuilder.cpp
index 3b3e3d80da6256092ba95bbd49f58e424cde8a56..21e413ff4e94519b79fbc26e87d742e67c0436c2 100644
--- a/src/mesh/DiamondDualMeshBuilder.cpp
+++ b/src/mesh/DiamondDualMeshBuilder.cpp
@@ -7,6 +7,7 @@
 #include <mesh/MeshData.hpp>
 #include <mesh/MeshDataManager.hpp>
 #include <mesh/PrimalToDiamondDualConnectivityDataMapper.hpp>
+#include <utils/Stringify.hpp>
 
 template <size_t Dimension>
 void
@@ -55,7 +56,7 @@ DiamondDualMeshBuilder::DiamondDualMeshBuilder(const IMesh& i_mesh)
   }
     // LCOV_EXCL_START
   default: {
-    throw UnexpectedError("invalid mesh dimension: " + std::to_string(i_mesh.dimension()));
+    throw UnexpectedError("invalid mesh dimension: " + stringify(i_mesh.dimension()));
   }
     // LCOV_EXCL_STOP
   }
diff --git a/src/mesh/Dual1DConnectivityBuilder.cpp b/src/mesh/Dual1DConnectivityBuilder.cpp
index ec0389d33ce17ba32141f7fcea8e251eba08ab70..30c465b8b07375fce9f82fb76a26ac24848f3c67 100644
--- a/src/mesh/Dual1DConnectivityBuilder.cpp
+++ b/src/mesh/Dual1DConnectivityBuilder.cpp
@@ -9,6 +9,7 @@
 #include <mesh/RefId.hpp>
 #include <utils/Array.hpp>
 #include <utils/Messenger.hpp>
+#include <utils/Stringify.hpp>
 
 #include <vector>
 
@@ -219,7 +220,7 @@ Dual1DConnectivityBuilder::Dual1DConnectivityBuilder(const IConnectivity& connec
     this->_buildConnectivityFrom(connectivity);
   } else {
     // LCOV_EXCL_START
-    throw UnexpectedError("invalid connectivity dimension: " + std::to_string(connectivity.dimension()));
+    throw UnexpectedError("invalid connectivity dimension: " + stringify(connectivity.dimension()));
     // LCOV_EXCL_STOP
   }
 }
diff --git a/src/mesh/Dual1DMeshBuilder.cpp b/src/mesh/Dual1DMeshBuilder.cpp
index ed558e2fc3c6c82861f93bc4364eb8dddf4000be..db47aa6d441619b1d1c5e2061a9838209851d5d9 100644
--- a/src/mesh/Dual1DMeshBuilder.cpp
+++ b/src/mesh/Dual1DMeshBuilder.cpp
@@ -8,6 +8,7 @@
 #include <mesh/MeshDataManager.hpp>
 #include <mesh/PrimalToDiamondDualConnectivityDataMapper.hpp>
 #include <mesh/PrimalToDual1DConnectivityDataMapper.hpp>
+#include <utils/Stringify.hpp>
 
 void
 Dual1DMeshBuilder::_buildDual1DMeshFrom(const IMesh& i_mesh)
@@ -46,7 +47,7 @@ Dual1DMeshBuilder::Dual1DMeshBuilder(const IMesh& i_mesh)
     this->_buildDual1DMeshFrom(i_mesh);
   } else {
     // LCOV_EXCL_START
-    throw UnexpectedError("invalid mesh dimension: " + std::to_string(i_mesh.dimension()));
+    throw UnexpectedError("invalid mesh dimension: " + stringify(i_mesh.dimension()));
     // LCOV_EXCL_STOP
   }
 }
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index f12e6ab3c6b17e3ae31110b91789ba8188b53d7c..6a2995687cbff74907766d694b5443ed95551b9a 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -10,6 +10,7 @@
 #include <mesh/Mesh.hpp>
 #include <mesh/RefItemList.hpp>
 #include <utils/Exceptions.hpp>
+#include <utils/Stringify.hpp>
 
 #include <rang.hpp>
 
@@ -710,16 +711,16 @@ GmshReader::GmshReader(const std::string& filename) : m_filename(filename)
     case MESHFORMAT: {
       double fileVersion = this->_getReal();
       if (fileVersion != 2.2) {
-        throw NormalError("Cannot read Gmsh format '" + std::to_string(fileVersion) + "'");
+        throw NormalError("Cannot read Gmsh format '" + stringify(fileVersion) + "'");
       }
       int fileType = this->_getInteger();
       if ((fileType < 0) or (fileType > 1)) {
-        throw NormalError("Cannot read Gmsh file type '" + std::to_string(fileType) + "'");
+        throw NormalError("Cannot read Gmsh file type '" + stringify(fileType) + "'");
       }
 
       int dataSize = this->_getInteger();
       if (dataSize != sizeof(double)) {
-        throw NormalError("Data size not supported '" + std::to_string(dataSize) + "'");
+        throw NormalError("Data size not supported '" + stringify(dataSize) + "'");
       }
 
       kw = this->__nextKeyword();
@@ -776,7 +777,7 @@ GmshReader::GmshReader(const std::string& filename) : m_filename(filename)
       break;
     }
     default: {
-      throw NormalError("invalid mesh dimension" + std::to_string((mesh_dimension)));
+      throw NormalError("invalid mesh dimension" + stringify((mesh_dimension)));
     }
     }
   }
@@ -822,8 +823,7 @@ GmshReader::__readElements2_2()
     const int elementType          = this->_getInteger() - 1;
 
     if ((elementType < 0) or (elementType > 14)) {
-      throw NormalError("reading file '" + m_filename + "': unknown element type '" + std::to_string(elementType) +
-                        "'");
+      throw NormalError("reading file '" + m_filename + "': unknown element type '" + stringify(elementType) + "'");
     }
     m_mesh_data.__elementType[i] = elementType;
     const int numberOfTags       = this->_getInteger();
@@ -1257,7 +1257,7 @@ GmshReader::__proceedData()
       const int a     = m_mesh_data.__verticesCorrepondance[elementVertices[0]];
       const int b     = m_mesh_data.__verticesCorrepondance[elementVertices[1]];
       if ((a < 0) or (b < 0)) {
-        throw NormalError("reading file '" + m_filename + "': error reading element " + std::to_string(i) +
+        throw NormalError("reading file '" + m_filename + "': error reading element " + stringify(i) +
                           " [bad vertices definition]");
       }
       m_mesh_data.__edges[edgeNumber]        = GmshData::Edge(a, b);
@@ -1273,7 +1273,7 @@ GmshReader::__proceedData()
       const int b = m_mesh_data.__verticesCorrepondance[elementVertices[1]];
       const int c = m_mesh_data.__verticesCorrepondance[elementVertices[2]];
       if ((a < 0) or (b < 0) or (c < 0)) {
-        throw NormalError("reading file '" + m_filename + "': error reading element " + std::to_string(i) +
+        throw NormalError("reading file '" + m_filename + "': error reading element " + stringify(i) +
                           " [bad vertices definition]");
       }
       m_mesh_data.__triangles[triangleNumber]        = GmshData::Triangle(a, b, c);
@@ -1290,7 +1290,7 @@ GmshReader::__proceedData()
       const int c = m_mesh_data.__verticesCorrepondance[elementVertices[2]];
       const int d = m_mesh_data.__verticesCorrepondance[elementVertices[3]];
       if ((a < 0) or (b < 0) or (c < 0) or (d < 0)) {
-        throw NormalError("reading file '" + m_filename + "': error reading element " + std::to_string(i) +
+        throw NormalError("reading file '" + m_filename + "': error reading element " + stringify(i) +
                           " [bad vertices definition]");
       }
       m_mesh_data.__quadrangles[quadrilateralNumber]        = GmshData::Quadrangle(a, b, c, d);
@@ -1307,7 +1307,7 @@ GmshReader::__proceedData()
       const int c = m_mesh_data.__verticesCorrepondance[elementVertices[2]];
       const int d = m_mesh_data.__verticesCorrepondance[elementVertices[3]];
       if ((a < 0) or (b < 0) or (c < 0) or (d < 0)) {
-        throw NormalError("reading file '" + m_filename + "': error reading element " + std::to_string(i) +
+        throw NormalError("reading file '" + m_filename + "': error reading element " + stringify(i) +
                           " [bad vertices definition]");
       }
       m_mesh_data.__tetrahedra[tetrahedronNumber]        = GmshData::Tetrahedron(a, b, c, d);
@@ -1328,7 +1328,7 @@ GmshReader::__proceedData()
       const int g = m_mesh_data.__verticesCorrepondance[elementVertices[6]];
       const int h = m_mesh_data.__verticesCorrepondance[elementVertices[7]];
       if ((a < 0) or (b < 0) or (c < 0) or (d < 0) or (e < 0) or (f < 0) or (g < 0) or (h < 0)) {
-        throw NormalError("reading file '" + m_filename + "': error reading element " + std::to_string(i) +
+        throw NormalError("reading file '" + m_filename + "': error reading element " + stringify(i) +
                           " [bad vertices definition]");
       }
       m_mesh_data.__hexahedra[hexahedronNumber]        = GmshData::Hexahedron(a, b, c, d, e, f, g, h);
@@ -1346,7 +1346,7 @@ GmshReader::__proceedData()
       const int e       = m_mesh_data.__verticesCorrepondance[elementVertices[4]];
       const int f       = m_mesh_data.__verticesCorrepondance[elementVertices[5]];
       if ((a < 0) or (b < 0) or (c < 0) or (d < 0) or (e < 0) or (f < 0)) {
-        throw NormalError("reading file '" + m_filename + "': error reading element " + std::to_string(i) +
+        throw NormalError("reading file '" + m_filename + "': error reading element " + stringify(i) +
                           " [bad vertices definition]");
       }
       m_mesh_data.__prisms[prism_number]        = GmshData::Prism(a, b, c, d, e, f);
@@ -1364,7 +1364,7 @@ GmshReader::__proceedData()
       const int d = m_mesh_data.__verticesCorrepondance[elementVertices[3]];
       const int e = m_mesh_data.__verticesCorrepondance[elementVertices[4]];
       if ((a < 0) or (b < 0) or (c < 0) or (d < 0) or (e < 0)) {
-        throw NormalError("reading file '" + m_filename + "': error reading element " + std::to_string(i) +
+        throw NormalError("reading file '" + m_filename + "': error reading element " + stringify(i) +
                           " [bad vertices definition]");
       }
       m_mesh_data.__pyramids[pyramid_number]        = GmshData::Pyramid(a, b, c, d, e);
diff --git a/src/mesh/MedianDualConnectivityBuilder.cpp b/src/mesh/MedianDualConnectivityBuilder.cpp
index 8ba8046c074355d702b6c3ec599d13232e080a0c..c6770d05ee1d93891d560bd1b1373fd2e68da5eb 100644
--- a/src/mesh/MedianDualConnectivityBuilder.cpp
+++ b/src/mesh/MedianDualConnectivityBuilder.cpp
@@ -9,6 +9,7 @@
 #include <mesh/RefId.hpp>
 #include <utils/Array.hpp>
 #include <utils/Messenger.hpp>
+#include <utils/Stringify.hpp>
 
 #include <vector>
 
@@ -417,7 +418,7 @@ MedianDualConnectivityBuilder::MedianDualConnectivityBuilder(const IConnectivity
   }
     // LCOV_EXCL_START
   default: {
-    throw UnexpectedError("invalid connectivity dimension: " + std::to_string(connectivity.dimension()));
+    throw UnexpectedError("invalid connectivity dimension: " + stringify(connectivity.dimension()));
   }
     // LCOV_EXCL_STOP
   }
diff --git a/src/mesh/MedianDualMeshBuilder.cpp b/src/mesh/MedianDualMeshBuilder.cpp
index aae78311ed7ab548714f7f6d90a96196632dcbb9..79d71dadc0e619457dadb205db053935283c170a 100644
--- a/src/mesh/MedianDualMeshBuilder.cpp
+++ b/src/mesh/MedianDualMeshBuilder.cpp
@@ -7,6 +7,7 @@
 #include <mesh/MeshData.hpp>
 #include <mesh/MeshDataManager.hpp>
 #include <mesh/PrimalToMedianDualConnectivityDataMapper.hpp>
+#include <utils/Stringify.hpp>
 
 template <size_t Dimension>
 void
@@ -54,7 +55,7 @@ MedianDualMeshBuilder::MedianDualMeshBuilder(const IMesh& i_mesh)
   }
     // LCOV_EXCL_START
   default: {
-    throw UnexpectedError("invalid mesh dimension: " + std::to_string(i_mesh.dimension()));
+    throw UnexpectedError("invalid mesh dimension: " + stringify(i_mesh.dimension()));
   }
     // LCOV_EXCL_STOP
   }
diff --git a/src/mesh/RefId.hpp b/src/mesh/RefId.hpp
index 307324ee8bb0515fe79f046cfbf32c1eea89ec44..010be2eb1d5a3b2198261a55d93742261cb5b0f5 100644
--- a/src/mesh/RefId.hpp
+++ b/src/mesh/RefId.hpp
@@ -1,6 +1,8 @@
 #ifndef REF_ID_HPP
 #define REF_ID_HPP
 
+#include <utils/Stringify.hpp>
+
 #include <iostream>
 #include <string>
 
@@ -18,7 +20,7 @@ class RefId
   friend std::ostream&
   operator<<(std::ostream& os, const RefId& ref_id)
   {
-    if (std::to_string(ref_id.m_tag_number) != ref_id.m_tag_name) {
+    if (stringify(ref_id.m_tag_number) != ref_id.m_tag_name) {
       os << ref_id.m_tag_name << '(' << ref_id.m_tag_number << ')';
     } else {
       os << ref_id.m_tag_number;
@@ -62,7 +64,7 @@ class RefId
     ;
   }
 
-  explicit RefId(TagNumberType tag_number) : m_tag_number(tag_number), m_tag_name(std::to_string(tag_number))
+  explicit RefId(TagNumberType tag_number) : m_tag_number(tag_number), m_tag_name(stringify(tag_number))
   {
     ;
   }
diff --git a/src/output/GnuplotWriter.cpp b/src/output/GnuplotWriter.cpp
index 324620fd67f3330bd92958fd8df8a8ccadfc5088..cf85b52eb9e020d165bdc8a30d86609ea1a17a3c 100644
--- a/src/output/GnuplotWriter.cpp
+++ b/src/output/GnuplotWriter.cpp
@@ -8,6 +8,7 @@
 #include <utils/Messenger.hpp>
 #include <utils/PugsTraits.hpp>
 #include <utils/RevisionInfo.hpp>
+#include <utils/Stringify.hpp>
 
 #include <utils/Demangle.hpp>
 
@@ -277,7 +278,7 @@ GnuplotWriter::writeMesh(const std::shared_ptr<const IMesh>& p_mesh) const
     break;
   }
   default: {
-    throw NormalError("gnuplot format is not available in dimension " + std::to_string(p_mesh->dimension()));
+    throw NormalError("gnuplot format is not available in dimension " + stringify(p_mesh->dimension()));
   }
   }
 }
@@ -300,7 +301,7 @@ GnuplotWriter::write(const std::vector<std::shared_ptr<const NamedDiscreteFuncti
     break;
   }
   default: {
-    throw NormalError("gnuplot format is not available in dimension " + std::to_string(mesh->dimension()));
+    throw NormalError("gnuplot format is not available in dimension " + stringify(mesh->dimension()));
   }
   }
 }
diff --git a/src/output/GnuplotWriter1D.cpp b/src/output/GnuplotWriter1D.cpp
index a57405d3aabc4ef23aefd6bd2b5f74f2f1681fa3..65bb4039a1c5296bf028a869d3c5d1774b1ddd6c 100644
--- a/src/output/GnuplotWriter1D.cpp
+++ b/src/output/GnuplotWriter1D.cpp
@@ -8,6 +8,7 @@
 #include <utils/Messenger.hpp>
 #include <utils/PugsTraits.hpp>
 #include <utils/RevisionInfo.hpp>
+#include <utils/Stringify.hpp>
 
 #include <utils/Demangle.hpp>
 
@@ -315,13 +316,13 @@ GnuplotWriter1D::write(const std::vector<std::shared_ptr<const NamedDiscreteFunc
   }
   case 2: {
     std::ostringstream errorMsg;
-    errorMsg << "gnuplot_1d_writer is not available in dimension " << std::to_string(mesh->dimension()) << '\n'
+    errorMsg << "gnuplot_1d_writer is not available in dimension " << stringify(mesh->dimension()) << '\n'
              << rang::style::bold << "note:" << rang::style::reset << " one can use " << rang::fgB::blue
-             << "gnuplot_writer" << rang::style::reset << " in dimension 2";
+             << "gnuplot_writer" << rang::fg::reset << " in dimension 2";
     throw NormalError(errorMsg.str());
   }
   default: {
-    throw NormalError("gnuplot format is not available in dimension " + std::to_string(mesh->dimension()));
+    throw NormalError("gnuplot format is not available in dimension " + stringify(mesh->dimension()));
   }
   }
 }
diff --git a/src/output/VTKWriter.cpp b/src/output/VTKWriter.cpp
index 98173a1402b588ceb9c9584bf5dbe84c136fc8e3..de7be9f664133bc7de2833893ddc6b19a99e774d 100644
--- a/src/output/VTKWriter.cpp
+++ b/src/output/VTKWriter.cpp
@@ -6,6 +6,7 @@
 #include <mesh/MeshDataManager.hpp>
 #include <utils/Messenger.hpp>
 #include <utils/RevisionInfo.hpp>
+#include <utils/Stringify.hpp>
 #include <utils/pugs_config.hpp>
 
 #include <ctime>
@@ -211,12 +212,12 @@ struct VTKWriter::VTKType
 
     if constexpr (std::is_integral_v<DataType>) {
       if constexpr (std::is_unsigned_v<DataType>) {
-        return "UInt" + std::to_string(sizeof(DataType) * 8);
+        return "UInt" + stringify(sizeof(DataType) * 8);
       } else {
-        return "UInt" + std::to_string(sizeof(DataType) * 8);
+        return "UInt" + stringify(sizeof(DataType) * 8);
       }
     } else if constexpr (std::is_floating_point_v<DataType>) {
-      return "Float" + std::to_string(sizeof(DataType) * 8);
+      return "Float" + stringify(sizeof(DataType) * 8);
     }
   }();
 };
diff --git a/src/scheme/CellIntegrator.hpp b/src/scheme/CellIntegrator.hpp
index 6960901d335f894e00a86239772f33f3abc50c67..26544b7c4f742dee16295922d9a9206320847844 100644
--- a/src/scheme/CellIntegrator.hpp
+++ b/src/scheme/CellIntegrator.hpp
@@ -14,6 +14,7 @@
 #include <mesh/Connectivity.hpp>
 #include <mesh/ItemId.hpp>
 #include <utils/Array.hpp>
+#include <utils/Stringify.hpp>
 
 #include <functional>
 
@@ -200,7 +201,7 @@ class CellIntegrator
           } else {
             // LCOV_EXCL_START
             throw NotImplementedError("integration on pyramid with non-quadrangular base (number of nodes " +
-                                      std::to_string(cell_to_node.size()) + ")");
+                                      stringify(cell_to_node.size()) + ")");
             // LCOV_EXCL_STOP
           }
           break;
@@ -255,7 +256,7 @@ class CellIntegrator
           } else {
             // LCOV_EXCL_START
             throw NotImplementedError("integration on diamond with non-quadrangular base (" +
-                                      std::to_string(cell_to_node.size()) + ")");
+                                      stringify(cell_to_node.size()) + ")");
             // LCOV_EXCL_STOP
           }
           break;
@@ -465,7 +466,7 @@ class CellIntegrator
           } else {
             // LCOV_EXCL_START
             throw NotImplementedError("integration on diamond with non-quadrangular base (" +
-                                      std::to_string(cell_to_node.size()) + ")");
+                                      stringify(cell_to_node.size()) + ")");
             // LCOV_EXCL_STOP
           }
           break;
diff --git a/src/scheme/DiscreteFunctionUtils.cpp b/src/scheme/DiscreteFunctionUtils.cpp
index f7057c4abcf920eea37d678eb410e7edf58006b3..fe861c0868b6bf620696cba09d8ac72bcf6aa948 100644
--- a/src/scheme/DiscreteFunctionUtils.cpp
+++ b/src/scheme/DiscreteFunctionUtils.cpp
@@ -4,6 +4,7 @@
 #include <mesh/IMesh.hpp>
 #include <mesh/Mesh.hpp>
 #include <scheme/DiscreteFunctionP0.hpp>
+#include <utils/Stringify.hpp>
 
 template <size_t Dimension, typename DataType>
 std::shared_ptr<const IDiscreteFunction>
@@ -52,8 +53,7 @@ shallowCopy(const std::shared_ptr<const Mesh<Connectivity<Dimension>>>& mesh,
       }
         // LCOV_EXCL_START
       default: {
-        throw UnexpectedError("invalid data vector dimension: " +
-                              std::to_string(discrete_function->dataType().dimension()));
+        throw UnexpectedError("invalid data vector dimension: " + stringify(discrete_function->dataType().dimension()));
       }
         // LCOV_EXCL_STOP
       }
@@ -62,8 +62,8 @@ shallowCopy(const std::shared_ptr<const Mesh<Connectivity<Dimension>>>& mesh,
       if (discrete_function->dataType().numberOfRows() != discrete_function->dataType().numberOfColumns()) {
         // LCOV_EXCL_START
         throw UnexpectedError(
-          "invalid data matrix dimensions: " + std::to_string(discrete_function->dataType().numberOfRows()) + "x" +
-          std::to_string(discrete_function->dataType().numberOfColumns()));
+          "invalid data matrix dimensions: " + stringify(discrete_function->dataType().numberOfRows()) + "x" +
+          stringify(discrete_function->dataType().numberOfColumns()));
         // LCOV_EXCL_STOP
       }
       switch (discrete_function->dataType().numberOfRows()) {
@@ -82,8 +82,8 @@ shallowCopy(const std::shared_ptr<const Mesh<Connectivity<Dimension>>>& mesh,
         // LCOV_EXCL_START
       default: {
         throw UnexpectedError(
-          "invalid data matrix dimensions: " + std::to_string(discrete_function->dataType().numberOfRows()) + "x" +
-          std::to_string(discrete_function->dataType().numberOfColumns()));
+          "invalid data matrix dimensions: " + stringify(discrete_function->dataType().numberOfRows()) + "x" +
+          stringify(discrete_function->dataType().numberOfColumns()));
       }
         // LCOV_EXCL_STOP
       }
diff --git a/src/utils/BuildInfo.cpp b/src/utils/BuildInfo.cpp
index 71d702c99ab86f7737b23a8ef42fa53845875cef..5a8fdb1bfdd66b209714a69b8f48c0de2fba07bd 100644
--- a/src/utils/BuildInfo.cpp
+++ b/src/utils/BuildInfo.cpp
@@ -1,4 +1,5 @@
 #include <utils/BuildInfo.hpp>
+#include <utils/Stringify.hpp>
 #include <utils/pugs_build_info.hpp>
 #include <utils/pugs_config.hpp>
 
@@ -55,8 +56,8 @@ std::string
 BuildInfo::petscLibrary()
 {
 #ifdef PUGS_HAS_PETSC
-  return std::to_string(PETSC_VERSION_MAJOR) + "." + std::to_string(PETSC_VERSION_MINOR) + "." +
-         std::to_string(PETSC_VERSION_SUBMINOR);
+  return stringify(PETSC_VERSION_MAJOR) + "." + stringify(PETSC_VERSION_MINOR) + "." +
+         stringify(PETSC_VERSION_SUBMINOR);
 #else
   return "none";
 #endif   // PUGS_HAS_PETSC
@@ -66,8 +67,8 @@ std::string
 BuildInfo::slepcLibrary()
 {
 #ifdef PUGS_HAS_SLEPC
-  return std::to_string(SLEPC_VERSION_MAJOR) + "." + std::to_string(SLEPC_VERSION_MINOR) + "." +
-         std::to_string(SLEPC_VERSION_SUBMINOR);
+  return stringify(SLEPC_VERSION_MAJOR) + "." + stringify(SLEPC_VERSION_MINOR) + "." +
+         stringify(SLEPC_VERSION_SUBMINOR);
 #else
   return "none";
 #endif   // PUGS_HAS_SLEPC
diff --git a/src/utils/Stringify.hpp b/src/utils/Stringify.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..7ff962c4f2e91313e65bf414127d60e3c0f8bef5
--- /dev/null
+++ b/src/utils/Stringify.hpp
@@ -0,0 +1,18 @@
+#ifndef STRINGIFY_HPP
+#define STRINGIFY_HPP
+
+#include <sstream>
+#include <string>
+
+// This utility is used in replacement of std::to_string. It allows a
+// better coherence of outputs especially in the case of double values
+template <typename T>
+std::string
+stringify(const T& value)
+{
+  std::ostringstream os;
+  os << std::boolalpha << value;
+  return os.str();
+}
+
+#endif   // STRINGIFY_HPP
diff --git a/tests/mpi_test_main.cpp b/tests/mpi_test_main.cpp
index 93f3c7d908c7e4d3840231ba90278123e387a945..51b222c08669e17f2e15e4b458f5210484e126b3 100644
--- a/tests/mpi_test_main.cpp
+++ b/tests/mpi_test_main.cpp
@@ -11,6 +11,7 @@
 #include <utils/Messenger.hpp>
 #include <utils/PETScWrapper.hpp>
 #include <utils/RandomEngine.hpp>
+#include <utils/Stringify.hpp>
 #include <utils/pugs_config.hpp>
 
 #include <MeshDataBaseForTests.hpp>
@@ -48,7 +49,7 @@ main(int argc, char* argv[])
       if (parallel::rank() != 0) {
         // Disable outputs for ranks != 0
         setenv("GCOV_PREFIX", gcov_prefix.string().c_str(), 1);
-        parallel_output /= output_base_name + std::to_string(parallel::rank());
+        parallel_output /= output_base_name + stringify(parallel::rank());
 
         Catch::ConfigData data{session.configData()};
         data.outputFilename = parallel_output.string();
@@ -74,7 +75,7 @@ main(int argc, char* argv[])
 
           for (size_t i_rank = 1; i_rank < parallel::size(); ++i_rank) {
             std::filesystem::path parallel_output = std::filesystem::path{PUGS_BINARY_DIR}.append("tests");
-            parallel_output /= output_base_name + std::to_string(i_rank);
+            parallel_output /= output_base_name + stringify(i_rank);
             session.config().stream() << " - " << rang::fg::green << parallel_output.parent_path().string()
                                       << parallel_output.preferred_separator << rang::style::reset << rang::fgB::green
                                       << parallel_output.filename().string() << rang::style::reset << '\n';
diff --git a/tests/test_BinaryExpressionProcessor_raw.cpp b/tests/test_BinaryExpressionProcessor_raw.cpp
index 22e32de89c715a5fa8e24a38007e0ab792e9dd19..f81ca5a0139e864ee47f358b943383b1d9101e3f 100644
--- a/tests/test_BinaryExpressionProcessor_raw.cpp
+++ b/tests/test_BinaryExpressionProcessor_raw.cpp
@@ -3,6 +3,8 @@
 
 #include <language/node_processor/BinaryExpressionProcessor.hpp>
 #include <language/utils/OFStream.hpp>
+#include <utils/Stringify.hpp>
+
 #include <utils/pugs_config.hpp>
 
 // clazy:excludeall=non-pod-global-static
@@ -53,7 +55,7 @@ TEST_CASE("BinaryExpressionProcessor raw operators", "[language]")
 
   {
     std::filesystem::path path{PUGS_BINARY_DIR};
-    path.append("tests").append(std::string{"binary_expression_processor_shift_left_"} + std::to_string(getpid()));
+    path.append("tests").append(std::string{"binary_expression_processor_shift_left_"} + stringify(getpid()));
 
     std::string filename = path.string();
 
diff --git a/tests/test_BinaryExpressionProcessor_shift.cpp b/tests/test_BinaryExpressionProcessor_shift.cpp
index 010201d122ea543b828a1de7060700aa78313593..71addbac7a1c80aaf49d1b236f1edfed1f182a4e 100644
--- a/tests/test_BinaryExpressionProcessor_shift.cpp
+++ b/tests/test_BinaryExpressionProcessor_shift.cpp
@@ -2,6 +2,8 @@
 #include <catch2/matchers/catch_matchers_all.hpp>
 
 #include <test_BinaryExpressionProcessor_utils.hpp>
+#include <utils/Stringify.hpp>
+
 #include <utils/pugs_config.hpp>
 
 #include <fstream>
@@ -16,7 +18,7 @@ TEST_CASE("BinaryExpressionProcessor shift", "[language]")
   SECTION("<<")
   {
     std::filesystem::path path{PUGS_BINARY_DIR};
-    path.append("tests").append(std::string{"binary_expression_processor_"} + std::to_string(getpid()));
+    path.append("tests").append(std::string{"binary_expression_processor_"} + stringify(getpid()));
 
     std::string filename = path.string();
 
diff --git a/tests/test_BuildInfo.cpp b/tests/test_BuildInfo.cpp
index 3713ef9d39ca5e74fcc62b47816cc51e1dbd41fd..847db9f3cdc8517bcca91b50329fb6ff70fbe3f7 100644
--- a/tests/test_BuildInfo.cpp
+++ b/tests/test_BuildInfo.cpp
@@ -2,6 +2,7 @@
 #include <catch2/matchers/catch_matchers_all.hpp>
 
 #include <utils/BuildInfo.hpp>
+#include <utils/Stringify.hpp>
 #include <utils/pugs_build_info.hpp>
 #include <utils/pugs_config.hpp>
 
@@ -55,8 +56,8 @@ TEST_CASE("BuildInfo", "[utils]")
   SECTION("petsc")
   {
 #ifdef PUGS_HAS_PETSC
-    const std::string petsc_library = std::to_string(PETSC_VERSION_MAJOR) + "." + std::to_string(PETSC_VERSION_MINOR) +
-                                      "." + std::to_string(PETSC_VERSION_SUBMINOR);
+    const std::string petsc_library =
+      stringify(PETSC_VERSION_MAJOR) + "." + stringify(PETSC_VERSION_MINOR) + "." + stringify(PETSC_VERSION_SUBMINOR);
 
     REQUIRE(BuildInfo::petscLibrary() == petsc_library);
 #else
diff --git a/tests/test_ConcatExpressionProcessor.cpp b/tests/test_ConcatExpressionProcessor.cpp
index d432ad0cc005d6568e02558ece048df180977c06..e7779f0bc3296773915de07d529b126b5380674c 100644
--- a/tests/test_ConcatExpressionProcessor.cpp
+++ b/tests/test_ConcatExpressionProcessor.cpp
@@ -11,6 +11,7 @@
 #include <language/ast/ASTSymbolTableBuilder.hpp>
 #include <language/utils/ASTPrinter.hpp>
 #include <utils/Demangle.hpp>
+#include <utils/Stringify.hpp>
 
 #include <pegtl/string_input.hpp>
 
@@ -79,13 +80,13 @@ TEST_CASE("ConcatExpressionProcessor", "[language]")
   SECTION("string + R")
   {
     CHECK_CONCAT_EXPRESSION_RESULT(R"(let s:string, s = "foo_"; s = s+2.4;)", "s",
-                                   std::string{"foo_"} + std::to_string(2.4));
+                                   std::string{"foo_"} + stringify(2.4));
   }
 
   SECTION("R + string")
   {
     CHECK_CONCAT_EXPRESSION_RESULT(R"(let s:string, s = "_foo"; s = 2.4+s;)", "s",
-                                   std::to_string(2.4) + std::string{"_foo"});
+                                   stringify(2.4) + std::string{"_foo"});
   }
 
   SECTION("string + B")
diff --git a/tests/test_CubeGaussQuadrature.cpp b/tests/test_CubeGaussQuadrature.cpp
index d4221fe322dc6dc4f898ac2563dc14d463b2e71f..a475a706c6b00a52a4bc69c43c25c340f0eceb8d 100644
--- a/tests/test_CubeGaussQuadrature.cpp
+++ b/tests/test_CubeGaussQuadrature.cpp
@@ -8,6 +8,7 @@
 #include <analysis/GaussQuadratureDescriptor.hpp>
 #include <analysis/QuadratureManager.hpp>
 #include <utils/Exceptions.hpp>
+#include <utils/Stringify.hpp>
 
 // clazy:excludeall=non-pod-global-static
 
@@ -491,7 +492,7 @@ TEST_CASE("CubeGaussQuadrature", "[analysis]")
     REQUIRE_THROWS_WITH(QuadratureManager::instance().getCubeFormula(
                           GaussQuadratureDescriptor(CubeGaussQuadrature ::max_degree + 1)),
                         "error: Gauss quadrature formulae handle degrees up to " +
-                          std::to_string(CubeGaussQuadrature ::max_degree) + " on cubes");
+                          stringify(CubeGaussQuadrature ::max_degree) + " on cubes");
   }
 
   SECTION("Access functions")
diff --git a/tests/test_FunctionArgumentConverter.cpp b/tests/test_FunctionArgumentConverter.cpp
index 05783b28bf4209e7c8bca30244581e2e386edabc..f02deac4a56de1f0e3f1cc81b027afaa72a3a0f6 100644
--- a/tests/test_FunctionArgumentConverter.cpp
+++ b/tests/test_FunctionArgumentConverter.cpp
@@ -3,6 +3,7 @@
 
 #include <language/node_processor/FunctionArgumentConverter.hpp>
 #include <language/utils/SymbolTable.hpp>
+#include <utils/Stringify.hpp>
 
 // clazy:excludeall=non-pod-global-static
 
@@ -29,7 +30,7 @@ TEST_CASE("FunctionArgumentConverter", "[language]")
 
     REQUIRE(std::get<std::string>(execution_policy.currentContext()[0]) == s);
     REQUIRE(std::get<std::string>(execution_policy.currentContext()[1]) == os_X.str());
-    REQUIRE(std::get<std::string>(execution_policy.currentContext()[2]) == std::to_string(x));
+    REQUIRE(std::get<std::string>(execution_policy.currentContext()[2]) == stringify(x));
   }
 
   SECTION("FunctionArgumentConverter")
diff --git a/tests/test_OFStream.cpp b/tests/test_OFStream.cpp
index edcee297572862e12d62a8a5032757ea84b364e2..fbb5818dbc70451ff212c820edacfa595963e961 100644
--- a/tests/test_OFStream.cpp
+++ b/tests/test_OFStream.cpp
@@ -3,6 +3,7 @@
 
 #include <language/utils/OFStream.hpp>
 #include <utils/Messenger.hpp>
+#include <utils/Stringify.hpp>
 
 #include <filesystem>
 
@@ -13,7 +14,7 @@ TEST_CASE("OFStream", "[language]")
   SECTION("ofstream")
   {
     const std::string basename = std::filesystem::path{PUGS_BINARY_DIR}.append("tests").append("ofstream_");
-    const std::string filename = basename + std::to_string(parallel::rank());
+    const std::string filename = basename + stringify(parallel::rank());
 
     // Ensures that the file is closed after this line
     std::make_shared<OFStream>(filename) << "foo" << 3 << " bar\n";
diff --git a/tests/test_PrismGaussQuadrature.cpp b/tests/test_PrismGaussQuadrature.cpp
index 000edeed6a40556d702bd440a26aab90f37b7e5d..f714ae4702f7feb51e544209d2d8c6666b0a57b6 100644
--- a/tests/test_PrismGaussQuadrature.cpp
+++ b/tests/test_PrismGaussQuadrature.cpp
@@ -7,6 +7,7 @@
 #include <analysis/QuadratureManager.hpp>
 #include <geometry/TetrahedronTransformation.hpp>
 #include <utils/Exceptions.hpp>
+#include <utils/Stringify.hpp>
 
 // clazy:excludeall=non-pod-global-static
 
@@ -685,6 +686,6 @@ TEST_CASE("PrismGaussQuadrature", "[analysis]")
     REQUIRE_THROWS_WITH(QuadratureManager::instance().getPrismFormula(
                           GaussQuadratureDescriptor(PrismGaussQuadrature ::max_degree + 1)),
                         "error: Gauss quadrature formulae handle degrees up to " +
-                          std::to_string(PrismGaussQuadrature ::max_degree) + " on prisms");
+                          stringify(PrismGaussQuadrature ::max_degree) + " on prisms");
   }
 }
diff --git a/tests/test_PyramidGaussQuadrature.cpp b/tests/test_PyramidGaussQuadrature.cpp
index ea022b1360908ccd3f270d3afda53b710a486488..6252b620652776d6ea71713b8faa5348427f9840 100644
--- a/tests/test_PyramidGaussQuadrature.cpp
+++ b/tests/test_PyramidGaussQuadrature.cpp
@@ -7,6 +7,7 @@
 #include <analysis/QuadratureManager.hpp>
 #include <geometry/TetrahedronTransformation.hpp>
 #include <utils/Exceptions.hpp>
+#include <utils/Stringify.hpp>
 
 // clazy:excludeall=non-pod-global-static
 
@@ -569,6 +570,6 @@ TEST_CASE("PyramidGaussQuadrature", "[analysis]")
     REQUIRE_THROWS_WITH(QuadratureManager::instance().getPyramidFormula(
                           GaussQuadratureDescriptor(PyramidGaussQuadrature ::max_degree + 1)),
                         "error: Gauss quadrature formulae handle degrees up to " +
-                          std::to_string(PyramidGaussQuadrature ::max_degree) + " on pyramids");
+                          stringify(PyramidGaussQuadrature ::max_degree) + " on pyramids");
   }
 }
diff --git a/tests/test_SquareGaussQuadrature.cpp b/tests/test_SquareGaussQuadrature.cpp
index 74f442a3154e23e4436c9b170a9649e72a61d4ae..210b8acea70037ed09bf13db8e0312d05e54b5d7 100644
--- a/tests/test_SquareGaussQuadrature.cpp
+++ b/tests/test_SquareGaussQuadrature.cpp
@@ -8,6 +8,7 @@
 #include <analysis/QuadratureManager.hpp>
 #include <analysis/SquareGaussQuadrature.hpp>
 #include <utils/Exceptions.hpp>
+#include <utils/Stringify.hpp>
 
 // clazy:excludeall=non-pod-global-static
 
@@ -462,7 +463,7 @@ TEST_CASE("SquareGaussQuadrature", "[analysis]")
     REQUIRE_THROWS_WITH(QuadratureManager::instance().getSquareFormula(
                           GaussQuadratureDescriptor(SquareGaussQuadrature ::max_degree + 1)),
                         "error: Gauss quadrature formulae handle degrees up to " +
-                          std::to_string(SquareGaussQuadrature ::max_degree) + " on squares");
+                          stringify(SquareGaussQuadrature ::max_degree) + " on squares");
   }
 
   SECTION("Access functions")
diff --git a/tests/test_TensorialGaussLegendreQuadrature.cpp b/tests/test_TensorialGaussLegendreQuadrature.cpp
index 342a4fbaa9ef8c54e274f7ad65800aee30986559..c2dc34b92072f70b7ba47b223a8a2038236f8a64 100644
--- a/tests/test_TensorialGaussLegendreQuadrature.cpp
+++ b/tests/test_TensorialGaussLegendreQuadrature.cpp
@@ -7,6 +7,7 @@
 #include <analysis/QuadratureManager.hpp>
 #include <analysis/TensorialGaussLegendreQuadrature.hpp>
 #include <utils/Exceptions.hpp>
+#include <utils/Stringify.hpp>
 
 // clazy:excludeall=non-pod-global-static
 
@@ -523,12 +524,12 @@ TEST_CASE("TensorialGaussLegendreQuadrature", "[analysis]")
       REQUIRE_THROWS_WITH(QuadratureManager::instance().getLineFormula(
                             GaussLegendreQuadratureDescriptor(TensorialGaussLegendreQuadrature<1>::max_degree + 1)),
                           "error: Gauss-Legendre quadrature formulae handle degrees up to " +
-                            std::to_string(TensorialGaussLegendreQuadrature<1>::max_degree) + " on lines");
+                            stringify(TensorialGaussLegendreQuadrature<1>::max_degree) + " on lines");
 
       REQUIRE_THROWS_WITH(QuadratureManager::instance().getLineFormula(
                             GaussQuadratureDescriptor(TensorialGaussLegendreQuadrature<1>::max_degree + 1)),
                           "error: Gauss quadrature formulae handle degrees up to " +
-                            std::to_string(TensorialGaussLegendreQuadrature<1>::max_degree) + " on lines");
+                            stringify(TensorialGaussLegendreQuadrature<1>::max_degree) + " on lines");
     }
   }
 
@@ -587,7 +588,7 @@ TEST_CASE("TensorialGaussLegendreQuadrature", "[analysis]")
       REQUIRE_THROWS_WITH(QuadratureManager::instance().getSquareFormula(
                             GaussLegendreQuadratureDescriptor(TensorialGaussLegendreQuadrature<2>::max_degree + 1)),
                           "error: Gauss-Legendre quadrature formulae handle degrees up to " +
-                            std::to_string(TensorialGaussLegendreQuadrature<2>::max_degree) + " on squares");
+                            stringify(TensorialGaussLegendreQuadrature<2>::max_degree) + " on squares");
     }
   }
 
@@ -652,7 +653,7 @@ TEST_CASE("TensorialGaussLegendreQuadrature", "[analysis]")
       REQUIRE_THROWS_WITH(QuadratureManager::instance().getCubeFormula(
                             GaussLegendreQuadratureDescriptor(TensorialGaussLegendreQuadrature<3>::max_degree + 1)),
                           "error: Gauss-Legendre quadrature formulae handle degrees up to " +
-                            std::to_string(TensorialGaussLegendreQuadrature<3>::max_degree) + " on cubes");
+                            stringify(TensorialGaussLegendreQuadrature<3>::max_degree) + " on cubes");
     }
   }
 }
diff --git a/tests/test_TensorialGaussLobattoQuadrature.cpp b/tests/test_TensorialGaussLobattoQuadrature.cpp
index 9ee4028cf0aa9b58d9f7478e91f2146f7f77a675..5c2e8869b6577a7a5ebf4704935ac5e058215a70 100644
--- a/tests/test_TensorialGaussLobattoQuadrature.cpp
+++ b/tests/test_TensorialGaussLobattoQuadrature.cpp
@@ -6,6 +6,7 @@
 #include <analysis/QuadratureManager.hpp>
 #include <analysis/TensorialGaussLobattoQuadrature.hpp>
 #include <utils/Exceptions.hpp>
+#include <utils/Stringify.hpp>
 
 // clazy:excludeall=non-pod-global-static
 
@@ -297,7 +298,7 @@ TEST_CASE("TensorialGaussLobattoQuadrature", "[analysis]")
       REQUIRE_THROWS_WITH(QuadratureManager::instance().getLineFormula(
                             GaussLobattoQuadratureDescriptor(TensorialGaussLobattoQuadrature<1>::max_degree + 1)),
                           "error: Gauss-Lobatto quadrature formulae handle degrees up to " +
-                            std::to_string(TensorialGaussLobattoQuadrature<1>::max_degree) + " on lines");
+                            stringify(TensorialGaussLobattoQuadrature<1>::max_degree) + " on lines");
     }
   }
 
@@ -356,7 +357,7 @@ TEST_CASE("TensorialGaussLobattoQuadrature", "[analysis]")
       REQUIRE_THROWS_WITH(QuadratureManager::instance().getSquareFormula(
                             GaussLobattoQuadratureDescriptor(TensorialGaussLobattoQuadrature<2>::max_degree + 1)),
                           "error: Gauss-Lobatto quadrature formulae handle degrees up to " +
-                            std::to_string(TensorialGaussLobattoQuadrature<2>::max_degree) + " on squares");
+                            stringify(TensorialGaussLobattoQuadrature<2>::max_degree) + " on squares");
     }
   }
 
@@ -421,7 +422,7 @@ TEST_CASE("TensorialGaussLobattoQuadrature", "[analysis]")
       REQUIRE_THROWS_WITH(QuadratureManager::instance().getCubeFormula(
                             GaussLobattoQuadratureDescriptor(TensorialGaussLobattoQuadrature<3>::max_degree + 1)),
                           "error: Gauss-Lobatto quadrature formulae handle degrees up to " +
-                            std::to_string(TensorialGaussLobattoQuadrature<3>::max_degree) + " on cubes");
+                            stringify(TensorialGaussLobattoQuadrature<3>::max_degree) + " on cubes");
     }
   }
 
diff --git a/tests/test_TetrahedronGaussQuadrature.cpp b/tests/test_TetrahedronGaussQuadrature.cpp
index a412c13688596e37916fd50a522b00a6aaae8485..a2e832542f164d02456d8e205d48445cba16d61a 100644
--- a/tests/test_TetrahedronGaussQuadrature.cpp
+++ b/tests/test_TetrahedronGaussQuadrature.cpp
@@ -7,6 +7,7 @@
 #include <analysis/TetrahedronGaussQuadrature.hpp>
 #include <geometry/TetrahedronTransformation.hpp>
 #include <utils/Exceptions.hpp>
+#include <utils/Stringify.hpp>
 
 // clazy:excludeall=non-pod-global-static
 
@@ -682,6 +683,6 @@ TEST_CASE("TetrahedronGaussQuadrature", "[analysis]")
     REQUIRE_THROWS_WITH(QuadratureManager::instance().getTetrahedronFormula(
                           GaussQuadratureDescriptor(TetrahedronGaussQuadrature ::max_degree + 1)),
                         "error: Gauss quadrature formulae handle degrees up to " +
-                          std::to_string(TetrahedronGaussQuadrature ::max_degree) + " on tetrahedra");
+                          stringify(TetrahedronGaussQuadrature ::max_degree) + " on tetrahedra");
   }
 }
diff --git a/tests/test_TriangleGaussQuadrature.cpp b/tests/test_TriangleGaussQuadrature.cpp
index 4ed18d60284861241861c45dd744a8e1b38c794b..6f467386163ffe4a4d5fa8c5d7bef944d53bb0b6 100644
--- a/tests/test_TriangleGaussQuadrature.cpp
+++ b/tests/test_TriangleGaussQuadrature.cpp
@@ -7,6 +7,7 @@
 #include <analysis/TriangleGaussQuadrature.hpp>
 #include <geometry/TriangleTransformation.hpp>
 #include <utils/Exceptions.hpp>
+#include <utils/Stringify.hpp>
 
 // clazy:excludeall=non-pod-global-static
 
@@ -627,6 +628,6 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
     REQUIRE_THROWS_WITH(QuadratureManager::instance().getTriangleFormula(
                           GaussQuadratureDescriptor(TriangleGaussQuadrature::max_degree + 1)),
                         "error: Gauss quadrature formulae handle degrees up to " +
-                          std::to_string(TriangleGaussQuadrature::max_degree) + " on triangles");
+                          stringify(TriangleGaussQuadrature::max_degree) + " on triangles");
   }
 }