From 59588481457e92772e17300e3fc5d92654eca4b1 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Wed, 27 Mar 2019 19:22:35 +0100
Subject: [PATCH] Dynare++: improve comments using Unicode (first batch)

[skip ci]
---
 dynare++/kord/faa_di_bruno.hh        |  4 +-
 dynare++/kord/korder.hh              | 14 +++---
 dynare++/kord/korder_stoch.hh        |  4 +-
 dynare++/sylv/cc/BlockDiagonal.cc    | 18 ++++----
 dynare++/sylv/cc/GeneralMatrix.cc    |  8 ++--
 dynare++/sylv/cc/GeneralMatrix.hh    | 67 ++++++++++++++--------------
 dynare++/sylv/cc/GeneralSylvester.cc |  2 +-
 dynare++/sylv/cc/GeneralSylvester.hh |  4 +-
 dynare++/sylv/cc/KronUtils.cc        | 18 +++-----
 dynare++/sylv/cc/KronUtils.hh        | 11 +++--
 dynare++/sylv/cc/QuasiTriangular.cc  | 22 ++++-----
 dynare++/sylv/cc/QuasiTriangular.hh  | 24 +++++-----
 dynare++/sylv/cc/SchurDecomp.hh      |  2 +-
 dynare++/sylv/cc/SchurDecompEig.cc   | 40 ++++++++---------
 dynare++/sylv/cc/SimilarityDecomp.cc | 15 ++++---
 dynare++/sylv/cc/SylvMatrix.cc       |  6 +--
 dynare++/sylv/cc/SylvMatrix.hh       | 30 ++++++-------
 dynare++/sylv/cc/SylvesterSolver.hh  |  4 +-
 dynare++/sylv/cc/SymSchurDecomp.cc   |  4 +-
 dynare++/sylv/cc/SymSchurDecomp.hh   | 15 +++----
 dynare++/sylv/cc/Vector.hh           | 18 ++++----
 dynare++/tl/cc/equivalence.cc        | 44 +++++++++---------
 dynare++/tl/cc/equivalence.hh        | 62 ++++++++++++-------------
 dynare++/tl/cc/fine_container.cc     |  2 +-
 dynare++/tl/cc/fine_container.hh     | 38 ++++++++--------
 dynare++/tl/cc/fs_tensor.cc          | 36 ++++++++-------
 dynare++/tl/cc/fs_tensor.hh          | 25 ++++++-----
 dynare++/tl/cc/ps_tensor.hh          |  4 +-
 dynare++/tl/cc/pyramid_prod.hh       |  4 +-
 dynare++/tl/cc/stack_container.cc    |  2 +-
 dynare++/tl/cc/stack_container.hh    | 41 ++++++++++-------
 dynare++/tl/cc/t_container.cc        |  2 +-
 dynare++/tl/cc/t_container.hh        |  6 +--
 dynare++/utils/cc/pascal_triangle.cc | 11 +++--
 dynare++/utils/cc/pascal_triangle.hh |  3 ++
 35 files changed, 307 insertions(+), 303 deletions(-)

diff --git a/dynare++/kord/faa_di_bruno.hh b/dynare++/kord/faa_di_bruno.hh
index 15726eca2..8dfffde3a 100644
--- a/dynare++/kord/faa_di_bruno.hh
+++ b/dynare++/kord/faa_di_bruno.hh
@@ -1,8 +1,8 @@
 // Copyright 2005, Ondra Kamenik
 
-// Faa Di Bruno evaluator
+// Faà Di Bruno evaluator
 
-/* This defines a class which implements Faa Di Bruno Formula
+/* This defines a class which implements Faà Di Bruno Formula
    $$\left[B_{s^k}\right]_{\alpha_1\ldots\alpha_l}=\left[f_{z^l}\right]_{\beta_1\ldots\beta_l}
    \sum_{c\in M_{l,k}}\prod_{m=1}^l\left[z_{s^k(c_m)}\right]^{\beta_m}_{c_m(\alpha)}$$
    where $s^k$ is a general symmetry of dimension $k$ and $z$ is a stack of
diff --git a/dynare++/kord/korder.hh b/dynare++/kord/korder.hh
index 9c7aee1ee..5c8a910e8 100644
--- a/dynare++/kord/korder.hh
+++ b/dynare++/kord/korder.hh
@@ -214,9 +214,9 @@ public:
      containers
    - |sylvesterSolve|: solve the sylvester equation (templated fold, and
      unfold)
-   - |faaDiBrunoZ|: calculates derivatives of $F$ by Faa Di Bruno for the
+   - |faaDiBrunoZ|: calculates derivatives of $F$ by Faà Di Bruno for the
      sparse container of system derivatives and $Z$ stack container
-   - |faaDiBrunoG|: calculates derivatives of $G$ by Faa Di Bruno for the
+   - |faaDiBrunoG|: calculates derivatives of $G$ by Faà Di Bruno for the
      dense container $g^{**}$ and $G$ stack
    - |recover_y|: recovers $g_{y^{*i}}$
    - |recover_yu|: recovers $g_{y^{*i}u^j}$
@@ -387,7 +387,7 @@ KOrder::insertDerivative(std::unique_ptr<typename ctraits<t>::Ttensor> der)
                                                                  ypart.nyss(), *der_ptr));
 }
 
-/* Here we implement Faa Di Bruno formula
+/* Here we implement Faà Di Bruno formula
    $$\sum_{l=1}^k\left[f_{z^l}\right]_{\gamma_1\ldots\gamma_l}
    \sum_{c\in M_{l,k}}\prod_{m=1}^l\left[z_{s(c_m)}\right]^{\gamma_m},
    $$
@@ -399,7 +399,7 @@ std::unique_ptr<typename ctraits<t>::Ttensor>
 KOrder::faaDiBrunoZ(const Symmetry &sym) const
 {
   JournalRecordPair pa(journal);
-  pa << "Faa Di Bruno Z container for " << sym << endrec;
+  pa << "Faà Di Bruno Z container for " << sym << endrec;
   auto res = std::make_unique<typename ctraits<t>::Ttensor>(ny, TensorDimens(sym, nvs));
   FaaDiBruno bruno(journal);
   bruno.calculate(Zstack<t>(), f, *res);
@@ -414,7 +414,7 @@ std::unique_ptr<typename ctraits<t>::Ttensor>
 KOrder::faaDiBrunoG(const Symmetry &sym) const
 {
   JournalRecordPair pa(journal);
-  pa << "Faa Di Bruno G container for " << sym << endrec;
+  pa << "Faà Di Bruno G container for " << sym << endrec;
   TensorDimens tdims(sym, nvs);
   auto res = std::make_unique<typename ctraits<t>::Ttensor>(ypart.nyss(), tdims);
   FaaDiBruno bruno(journal);
@@ -762,9 +762,9 @@ KOrder::calcE_k(int k) const
 }
 
 /* Here is the core routine. It calls methods recovering derivatives in
-   the right order. Recall, that the code, namely Faa Di Bruno's formula,
+   the right order. Recall, that the code, namely Faà Di Bruno's formula,
    is implemented as to be run conditionally on the current contents of
-   containers. So, if some call of Faa Di Bruno evaluates derivatives,
+   containers. So, if some call of Faà Di Bruno evaluates derivatives,
    and some derivatives are not present in the container, then it is
    considered to be zero. So, we have to be very careful to put
    everything in the right order. The order here can be derived from
diff --git a/dynare++/kord/korder_stoch.hh b/dynare++/kord/korder_stoch.hh
index 2adf37cd0..27a1589d0 100644
--- a/dynare++/kord/korder_stoch.hh
+++ b/dynare++/kord/korder_stoch.hh
@@ -479,7 +479,7 @@ std::unique_ptr<typename ctraits<t>::Ttensor>
 KOrderStoch::faaDiBrunoZ(const Symmetry &sym) const
 {
   JournalRecordPair pa(journal);
-  pa << "Faa Di Bruno ZX container for " << sym << endrec;
+  pa << "Faà Di Bruno ZX container for " << sym << endrec;
   auto res = std::make_unique<typename ctraits<t>::Ttensor>(ypart.ny(), TensorDimens(sym, nvs));
   FaaDiBruno bruno(journal);
   bruno.calculate(Zstack<t>(), f, *res);
@@ -494,7 +494,7 @@ std::unique_ptr<typename ctraits<t>::Ttensor>
 KOrderStoch::faaDiBrunoG(const Symmetry &sym) const
 {
   JournalRecordPair pa(journal);
-  pa << "Faa Di Bruno GX container for " << sym << endrec;
+  pa << "Faà Di Bruno GX container for " << sym << endrec;
   TensorDimens tdims(sym, nvs);
   auto res = std::make_unique<typename ctraits<t>::Ttensor>(ypart.nyss(), tdims);
   FaaDiBruno bruno(journal);
diff --git a/dynare++/sylv/cc/BlockDiagonal.cc b/dynare++/sylv/cc/BlockDiagonal.cc
index 30a4a96c6..54597e97d 100644
--- a/dynare++/sylv/cc/BlockDiagonal.cc
+++ b/dynare++/sylv/cc/BlockDiagonal.cc
@@ -36,8 +36,8 @@ BlockDiagonal::BlockDiagonal(int p, const BlockDiagonal &b)
 {
 }
 
-/* put zeroes to right upper submatrix whose first column is defined
- * by 'edge' */
+/* Put zeroes to right upper submatrix whose first column is defined
+   by ‘edge’ */
 void
 BlockDiagonal::setZerosToRU(diag_iter edge)
 {
@@ -47,14 +47,14 @@ BlockDiagonal::setZerosToRU(diag_iter edge)
       get(i, j) = 0.0;
 }
 
-/* Updates row_len and col_len so that there are zeroes in upper right part, this
- * |T1 0 |
- * |0  T2|. The first column of T2 is given by diagonal iterator 'edge'.
+/* Updates row_len and col_len so that there are zeroes in upper right part, i.e.
+   ⎡ T1  0 ⎤
+   ⎣  0 T2 ⎦. The first column of T2 is given by diagonal iterator ‘edge’.
 
- * Note the semantics of row_len and col_len. row_len[i] is distance
- * of the right-most non-zero element of i-th row from the left, and
- * col_len[j] is distance of top-most non-zero element of j-th column
- * to the top. (First element has distance 1).
+   Note the semantics of row_len and col_len. row_len[i] is distance
+   of the right-most non-zero element of i-th row from the left, and
+   col_len[j] is distance of top-most non-zero element of j-th column
+   to the top. (First element has distance 1).
  */
 void
 BlockDiagonal::setZeroBlockEdge(diag_iter edge)
diff --git a/dynare++/sylv/cc/GeneralMatrix.cc b/dynare++/sylv/cc/GeneralMatrix.cc
index 26e06f199..8c1910d69 100644
--- a/dynare++/sylv/cc/GeneralMatrix.cc
+++ b/dynare++/sylv/cc/GeneralMatrix.cc
@@ -329,13 +329,13 @@ GeneralMatrix::gemm_partial_right(const std::string &trans, const ConstGeneralMa
 ConstGeneralMatrix::ConstGeneralMatrix(const GeneralMatrix &m, int i, int j, int nrows, int ncols)
   : data(m.getData(), j*m.getLD()+i, (ncols-1)*m.getLD()+nrows), rows(nrows), cols(ncols), ld(m.getLD())
 {
-  // can check that the submatirx is fully in the matrix
+  // FIXME: check that the submatrix is fully in the matrix
 }
 
 ConstGeneralMatrix::ConstGeneralMatrix(const ConstGeneralMatrix &m, int i, int j, int nrows, int ncols)
   : data(m.getData(), j*m.getLD()+i, (ncols-1)*m.getLD()+nrows), rows(nrows), cols(ncols), ld(m.getLD())
 {
-  // can check that the submatirx is fully in the matrix
+  // FIXME: check that the submatrix is fully in the matrix
 }
 
 ConstGeneralMatrix::ConstGeneralMatrix(const GeneralMatrix &m)
@@ -559,8 +559,8 @@ SVDDecomp::solve(const ConstGeneralMatrix &B, GeneralMatrix &X) const
   if (B.numRows() != U.numRows())
     throw SYLV_MES_EXCEPTION("Incompatible number of rows ");
 
-  // reciprocal condition number for determination of zeros in the
-  // end of sigma
+  /* Reciprocal condition number for determination of zeros in the
+     end of sigma */
   constexpr double rcond = 1e-13;
 
   // determine nz=number of zeros in the end of sigma
diff --git a/dynare++/sylv/cc/GeneralMatrix.hh b/dynare++/sylv/cc/GeneralMatrix.hh
index d87808015..715ccbdb7 100644
--- a/dynare++/sylv/cc/GeneralMatrix.hh
+++ b/dynare++/sylv/cc/GeneralMatrix.hh
@@ -282,7 +282,7 @@ public:
     place(ConstGeneralMatrix(m), i, j);
   }
 
-  /* this = a*b */
+  // this = a·b
   void mult(const ConstGeneralMatrix &a, const ConstGeneralMatrix &b);
   void
   mult(const GeneralMatrix &a, const GeneralMatrix &b)
@@ -290,7 +290,7 @@ public:
     mult(ConstGeneralMatrix(a), ConstGeneralMatrix(b));
   }
 
-  /* this = this + scalar*a*b */
+  // this = this + scalar·a·b
   void multAndAdd(const ConstGeneralMatrix &a, const ConstGeneralMatrix &b,
                   double mult = 1.0);
   void
@@ -300,7 +300,7 @@ public:
     multAndAdd(ConstGeneralMatrix(a), ConstGeneralMatrix(b), mult);
   }
 
-  /* this = this + scalar*a*b' */
+  // this = this + scalar·a·bᵀ
   void multAndAdd(const ConstGeneralMatrix &a, const ConstGeneralMatrix &b,
                   const std::string &dum, double mult = 1.0);
   void
@@ -310,7 +310,7 @@ public:
     multAndAdd(ConstGeneralMatrix(a), ConstGeneralMatrix(b), dum, mult);
   }
 
-  /* this = this + scalar*a'*b */
+  // this = this + scalar·aᵀ·b
   void multAndAdd(const ConstGeneralMatrix &a, const std::string &dum, const ConstGeneralMatrix &b,
                   double mult = 1.0);
   void
@@ -320,7 +320,7 @@ public:
     multAndAdd(ConstGeneralMatrix(a), dum, ConstGeneralMatrix(b), mult);
   }
 
-  /* this = this + scalar*a'*b' */
+  // this = this + scalar·aᵀ·bᵀ
   void multAndAdd(const ConstGeneralMatrix &a, const std::string &dum1,
                   const ConstGeneralMatrix &b, const std::string &dum2, double mult = 1.0);
   void
@@ -330,7 +330,7 @@ public:
     multAndAdd(ConstGeneralMatrix(a), dum1, ConstGeneralMatrix(b), dum2, mult);
   }
 
-  /* this = this + scalar*a*a' */
+  // this = this + scalar·a·aᵀ
   void addOuter(const ConstVector &a, double mult = 1.0);
   void
   addOuter(const Vector &a, double mult = 1.0)
@@ -338,7 +338,7 @@ public:
     addOuter(ConstVector(a), mult);
   }
 
-  /* this = this * m */
+  // this = this·m
   void multRight(const ConstGeneralMatrix &m);
   void
   multRight(const GeneralMatrix &m)
@@ -346,7 +346,7 @@ public:
     multRight(ConstGeneralMatrix(m));
   }
 
-  /* this = m * this */
+  // this = m·this
   void multLeft(const ConstGeneralMatrix &m);
   void
   multLeft(const GeneralMatrix &m)
@@ -354,7 +354,7 @@ public:
     multLeft(ConstGeneralMatrix(m));
   }
 
-  /* this = this * m' */
+  // this = this·mᵀ
   void multRightTrans(const ConstGeneralMatrix &m);
   void
   multRightTrans(const GeneralMatrix &m)
@@ -362,7 +362,7 @@ public:
     multRightTrans(ConstGeneralMatrix(m));
   }
 
-  /* this = m' * this */
+  // this = mᵀ·this
   void multLeftTrans(const ConstGeneralMatrix &m);
   void
   multLeftTrans(const GeneralMatrix &m)
@@ -370,64 +370,64 @@ public:
     multLeftTrans(ConstGeneralMatrix(m));
   }
 
-  /* x = scalar(a)*x + scalar(b)*this*d */
+  // x = scalar(a)·x + scalar(b)·this·d
   void
   multVec(double a, Vector &x, double b, const ConstVector &d) const
   {
     ConstGeneralMatrix(*this).multVec(a, x, b, d);
   }
 
-  /* x = scalar(a)*x + scalar(b)*this'*d */
+  // x = scalar(a)·x + scalar(b)·thisᵀ·d
   void
   multVecTrans(double a, Vector &x, double b, const ConstVector &d) const
   {
     ConstGeneralMatrix(*this).multVecTrans(a, x, b, d);
   }
 
-  /* x = x + this*d */
+  // x = x + this·d
   void
   multaVec(Vector &x, const ConstVector &d) const
   {
     ConstGeneralMatrix(*this).multaVec(x, d);
   }
 
-  /* x = x + this'*d */
+  // x = x + thisᵀ·d */
   void
   multaVecTrans(Vector &x, const ConstVector &d) const
   {
     ConstGeneralMatrix(*this).multaVecTrans(x, d);
   }
 
-  /* x = x - this*d */
+  // x = x - this·d
   void
   multsVec(Vector &x, const ConstVector &d) const
   {
     ConstGeneralMatrix(*this).multsVec(x, d);
   }
 
-  /* x = x - this'*d */
+  // x = x - thisᵀ·d
   void
   multsVecTrans(Vector &x, const ConstVector &d) const
   {
     ConstGeneralMatrix(*this).multsVecTrans(x, d);
   }
 
-  /* this = zero */
+  // this = zero
   void zeros();
 
-  /** this = unit (on main diagonal) */
+  // this = unit (on main diagonal)
   void unit();
 
-  /* this = NaN */
+  // this = NaN
   void nans();
 
-  /* this = Inf */
+  // this = ∞
   void infs();
 
-  /* this = scalar*this */
+  // this = scalar·this
   void mult(double a);
 
-  /* this = this + scalar*m */
+  // this = this + scalar·m
   void add(double a, const ConstGeneralMatrix &m);
   void
   add(double a, const GeneralMatrix &m)
@@ -435,7 +435,7 @@ public:
     add(a, ConstGeneralMatrix(m));
   }
 
-  /* this = this + scalar*m' */
+  // this = this + scalar·mᵀ
   void add(double a, const ConstGeneralMatrix &m, const std::string &dum);
   void
   add(double a, const GeneralMatrix &m, const std::string &dum)
@@ -490,7 +490,7 @@ private:
     gemm_partial_right(trans, ConstGeneralMatrix(m), alpha, beta);
   }
 
-  /* this = op(m) *this (without whole copy of this) */
+  // this = op(m)·this (without whole copy of ‘this’)
   void gemm_partial_left(const std::string &trans, const ConstGeneralMatrix &m,
                          double alpha, double beta);
   void
@@ -504,7 +504,7 @@ private:
   static constexpr int md_length = 23;
 };
 
-// Computes a*b
+// Computes a·b
 inline GeneralMatrix
 operator*(const ConstGeneralMatrix &a, const ConstGeneralMatrix &b)
 {
@@ -513,7 +513,7 @@ operator*(const ConstGeneralMatrix &a, const ConstGeneralMatrix &b)
   return m;
 }
 
-// Computes a*b'
+// Computes a·bᵀ
 template<class T>
 GeneralMatrix
 operator*(const ConstGeneralMatrix &a, const TransposedMatrix<T> &b)
@@ -523,7 +523,7 @@ operator*(const ConstGeneralMatrix &a, const TransposedMatrix<T> &b)
   return m;
 }
 
-// Computes a'*b
+// Computes aᵀ·b
 template<class T>
 GeneralMatrix
 operator*(const TransposedMatrix<T> &a, const ConstGeneralMatrix &b)
@@ -533,7 +533,7 @@ operator*(const TransposedMatrix<T> &a, const ConstGeneralMatrix &b)
   return m;
 }
 
-// Computes a'*b'
+// Computes aᵀ·bᵀ
 template<class T1, class T2>
 GeneralMatrix
 operator*(const TransposedMatrix<T1> &a, const TransposedMatrix<T2> &b)
@@ -546,16 +546,15 @@ operator*(const TransposedMatrix<T1> &a, const TransposedMatrix<T2> &b)
 class SVDDecomp
 {
 protected:
-  /** Minimum of number of rows and columns of the decomposed
-   * matrix. */
+  // Minimum of number of rows and columns of the decomposed matrix
   const int minmn;
-  /** Singular values. */
+  // Singular values
   Vector sigma;
-  /** Orthogonal matrix U. */
+  // Orthogonal matrix U
   GeneralMatrix U;
-  /** Orthogonal matrix V^T. */
+  // Orthogonal matrix Vᵀ
   GeneralMatrix VT;
-  /** Convered flag. */
+  // Convered flag
   bool conv;
 public:
   SVDDecomp(const GeneralMatrix &A)
diff --git a/dynare++/sylv/cc/GeneralSylvester.cc b/dynare++/sylv/cc/GeneralSylvester.cc
index f1043ee9b..4f40278ca 100644
--- a/dynare++/sylv/cc/GeneralSylvester.cc
+++ b/dynare++/sylv/cc/GeneralSylvester.cc
@@ -107,7 +107,7 @@ GeneralSylvester::check(const ConstVector &ds)
   if (!solved)
     throw SYLV_MES_EXCEPTION("Cannot run check on system, which is not solved yet.");
 
-  // calculate xcheck = AX+BXC^i-D
+  // calculate xcheck = A·X+B·X·⊗ⁱC-D
   SylvMatrix dcheck(d.numRows(), d.numCols());
   dcheck.multLeft(b.numRows()-b.numCols(), b, d);
   dcheck.multRightKron(c, order);
diff --git a/dynare++/sylv/cc/GeneralSylvester.hh b/dynare++/sylv/cc/GeneralSylvester.hh
index 79cceac50..2111af2aa 100644
--- a/dynare++/sylv/cc/GeneralSylvester.hh
+++ b/dynare++/sylv/cc/GeneralSylvester.hh
@@ -24,7 +24,7 @@ class GeneralSylvester
   std::unique_ptr<SimilarityDecomp> cdecomp;
   std::unique_ptr<SylvesterSolver> sylv;
 public:
-  /* construct with my copy of d*/
+  // Construct with my copy of d
   GeneralSylvester(int ord, int n, int m, int zero_cols,
                    const ConstVector &da, const ConstVector &db,
                    const ConstVector &dc, const ConstVector &dd,
@@ -33,7 +33,7 @@ public:
                    const ConstVector &da, const ConstVector &db,
                    const ConstVector &dc, const ConstVector &dd,
                    bool alloc_for_check = false);
-  /* construct with provided storage for d */
+  // Construct with provided storage for d
   GeneralSylvester(int ord, int n, int m, int zero_cols,
                    const ConstVector &da, const ConstVector &db,
                    const ConstVector &dc, Vector &dd,
diff --git a/dynare++/sylv/cc/KronUtils.cc b/dynare++/sylv/cc/KronUtils.cc
index 3bdcbcd3a..1115b38f4 100644
--- a/dynare++/sylv/cc/KronUtils.cc
+++ b/dynare++/sylv/cc/KronUtils.cc
@@ -10,13 +10,11 @@ KronUtils::multAtLevel(int level, const QuasiTriangular &t,
                        KronVector &x)
 {
   if (0 < level && level < x.getDepth())
-    {
-      for (int i = 0; i < x.getM(); i++)
-        {
-          KronVector xi(x, i);
-          multAtLevel(level, t, xi);
-        }
-    }
+    for (int i = 0; i < x.getM(); i++)
+      {
+        KronVector xi(x, i);
+        multAtLevel(level, t, xi);
+      }
   else if (0 == level && 0 < x.getDepth())
     {
       GeneralMatrix tmp(x, x.getN(), power(x.getM(), x.getDepth()));
@@ -61,8 +59,6 @@ KronUtils::multKron(const QuasiTriangular &f, const QuasiTriangular &k,
 {
   multAtLevel(0, k, x);
   if (x.getDepth() > 0)
-    {
-      for (int level = 1; level <= x.getDepth(); level++)
-        multAtLevelTrans(level, f, x);
-    }
+    for (int level = 1; level <= x.getDepth(); level++)
+      multAtLevelTrans(level, f, x);
 }
diff --git a/dynare++/sylv/cc/KronUtils.hh b/dynare++/sylv/cc/KronUtils.hh
index dcd7c2277..c7935e676 100644
--- a/dynare++/sylv/cc/KronUtils.hh
+++ b/dynare++/sylv/cc/KronUtils.hh
@@ -11,17 +11,16 @@
 class KronUtils
 {
 public:
-  /* multiplies I_m\otimes..\I_m\otimes T\otimes I_m...I_m\otimes I_n
-     with given b and returns x. T must be (m,m), number of
-     \otimes is b.getDepth(), level is a number of I_m's between T
-     and I_n plus 1. If level=0, then we multiply
-     \I_m\otimes ..\otimes I_m\otimes T, T is (n,n) */
+  /* Computes x = (Iₘ⊗…⊗Iₘ⊗T⊗Iₘ⊗…⊗Iₘ⊗Iₙ)·x, where x is n×mᵈ.
+     T must be m×m, number of ⊗ is d, level is the number of Iₘ’s
+     between T and Iₙ plus 1. If level=0, then we multiply by Iₘ⊗…⊗Iₘ⊗T,
+     T must be n×n. */
   static void multAtLevel(int level, const QuasiTriangular &t,
                           KronVector &x);
   static void multAtLevelTrans(int level, const QuasiTriangular &t,
                                KronVector &x);
 
-  /* multiplies x=(F'\otimes F'\otimes..\otimes K)x */
+  // Computes x=(Fᵀ⊗Fᵀ⊗…⊗K)·x
   static void multKron(const QuasiTriangular &f, const QuasiTriangular &k,
                        KronVector &x);
 };
diff --git a/dynare++/sylv/cc/QuasiTriangular.cc b/dynare++/sylv/cc/QuasiTriangular.cc
index bd633d7a5..fc9249ff9 100644
--- a/dynare++/sylv/cc/QuasiTriangular.cc
+++ b/dynare++/sylv/cc/QuasiTriangular.cc
@@ -34,9 +34,9 @@ DiagonalBlock::getSize() const
     return std::sqrt(getDeterminant());
 }
 
-// this function makes Diagonal inconsistent, it should only be used
-// on temorary matrices, which will not be used any more, e.g. in
-// QuasiTriangular::solve (we need fast performance)
+/* This function makes Diagonal inconsistent, it should only be used
+   on temorary matrices, which will not be used any more, e.g. in
+   QuasiTriangular::solve() (we need fast performance) */
 void
 DiagonalBlock::setReal()
 {
@@ -64,7 +64,7 @@ DiagonalBlock::checkBlock(const double *d, int d_size)
 
 Diagonal::Diagonal(double *data, int d_size)
 {
-  int nc = getNumComplex(data, d_size); // return nc <= d_size/2
+  int nc = getNumComplex(data, d_size); // return nc ≤ d_size/2
   num_all = d_size - nc;
   num_real = d_size - 2*nc;
 
@@ -184,9 +184,9 @@ Diagonal::getEigenValues(Vector &eig) const
     }
 }
 
-/* swaps logically blocks 'it', and '++it'. remember to move also
- * addresses, alpha, beta1, beta2. This is a dirty (but most
- * effective) way how to do it. */
+/* Swaps logically blocks ‘it’, and ‘++it’. remember to move also
+   addresses, alpha, beta1, beta2. This is a dirty (but most
+   effective) way how to do it. */
 void
 Diagonal::swapLogically(diag_iter it)
 {
@@ -432,7 +432,7 @@ QuasiTriangular::QuasiTriangular(const SchurDecomp &decomp)
 {
 }
 
-/* this pads matrix with intial columns with zeros */
+// This pads matrix with intial columns with zeros
 QuasiTriangular::QuasiTriangular(const SchurDecompZero &decomp)
   : SqSylvMatrix(decomp.getDim())
 {
@@ -590,7 +590,7 @@ QuasiTriangular::solvePreTrans(Vector &x, double &eig_min)
   dtrsv("U", "T", "N", &nn, getData().base(), &lda, x.base(), &incx);
 }
 
-/* calculates x = Tb */
+// Calculates x = T·b
 void
 QuasiTriangular::multVec(Vector &x, const ConstVector &b) const
 {
@@ -639,7 +639,7 @@ QuasiTriangular::multaVecTrans(Vector &x, const ConstVector &b) const
   x.add(1.0, tmp);
 }
 
-/* calculates x=x+(T\otimes I)b, where size of I is given by b (KronVector) */
+// Calculates x=x+(T⊗I)·b, where size of I is given by b (KronVector)
 void
 QuasiTriangular::multaKron(KronVector &x, const ConstKronVector &b) const
 {
@@ -649,7 +649,7 @@ QuasiTriangular::multaKron(KronVector &x, const ConstKronVector &b) const
   x_resh.multAndAdd(b_resh, ConstGeneralMatrix(*this), "trans");
 }
 
-/* calculates x=x+(T'\otimes I)b, where size of I is given by b (KronVector) */
+// Calculates x=x+(T'⊗I)·b, where size of I is given by b (KronVector)
 void
 QuasiTriangular::multaKronTrans(KronVector &x, const ConstKronVector &b) const
 {
diff --git a/dynare++/sylv/cc/QuasiTriangular.hh b/dynare++/sylv/cc/QuasiTriangular.hh
index 1ba51081f..5a1c9dd08 100644
--- a/dynare++/sylv/cc/QuasiTriangular.hh
+++ b/dynare++/sylv/cc/QuasiTriangular.hh
@@ -335,29 +335,29 @@ public:
   diag_iter findClosestDiagBlock(diag_iter start, diag_iter end, double a);
   diag_iter findNextLargerBlock(diag_iter start, diag_iter end, double a);
 
-  /* (I+T)y = x, y-->x  */
+  /* (I+T)·y = x, y→x  */
   virtual void solvePre(Vector &x, double &eig_min);
-  /* (I+T')y = x, y-->x */
+  /* (I+T')·y = x, y→x */
   virtual void solvePreTrans(Vector &x, double &eig_min);
   /* (I+T)x = b */
   virtual void solve(Vector &x, const ConstVector &b, double &eig_min);
   /* (I+T')x = b */
   virtual void solveTrans(Vector &x, const ConstVector &b, double &eig_min);
-  /* x = Tb */
+  /* x = T·b */
   virtual void multVec(Vector &x, const ConstVector &b) const;
-  /* x = T'b */
+  /* x = T'·b */
   virtual void multVecTrans(Vector &x, const ConstVector &b) const;
-  /* x = x + Tb */
+  /* x = x + T·b */
   virtual void multaVec(Vector &x, const ConstVector &b) const;
-  /* x = x + T'b */
+  /* x = x + T'·b */
   virtual void multaVecTrans(Vector &x, const ConstVector &b) const;
-  /* x = (T\otimes I)x */
+  /* x = (T⊗I)·x */
   virtual void multKron(KronVector &x) const;
-  /* x = (T'\otimes I)x */
+  /* x = (T'⊗I)·x */
   virtual void multKronTrans(KronVector &x) const;
-  /* A = T*A */
+  /* A = T·A */
   virtual void multLeftOther(GeneralMatrix &a) const;
-  /* A = T'*A */
+  /* A = T'·A */
   virtual void multLeftOtherTrans(GeneralMatrix &a) const;
 
   const_diag_iter
@@ -417,9 +417,9 @@ protected:
   void addMatrix(double r, const QuasiTriangular &t);
 private:
   void addUnit();
-  /* x = x + (T\otimes I)b */
+  /* x = x + (T⊗I)·b */
   void multaKron(KronVector &x, const ConstKronVector &b) const;
-  /* x = x + (T'\otimes I)b */
+  /* x = x + (T'⊗)·b */
   void multaKronTrans(KronVector &x, const ConstKronVector &b) const;
   /* implementation via iterators, useful for large matrices */
   void setMatrixViaIter(double r, const QuasiTriangular &t);
diff --git a/dynare++/sylv/cc/SchurDecomp.hh b/dynare++/sylv/cc/SchurDecomp.hh
index a8af976f5..efc643c01 100644
--- a/dynare++/sylv/cc/SchurDecomp.hh
+++ b/dynare++/sylv/cc/SchurDecomp.hh
@@ -47,7 +47,7 @@ public:
 
 class SchurDecompZero : public SchurDecomp
 {
-  GeneralMatrix ru; /* right upper matrix */
+  GeneralMatrix ru; // right upper matrix
 public:
   SchurDecompZero(const GeneralMatrix &m);
   ConstGeneralMatrix
diff --git a/dynare++/sylv/cc/SchurDecompEig.cc b/dynare++/sylv/cc/SchurDecompEig.cc
index 452015b00..5a679e511 100644
--- a/dynare++/sylv/cc/SchurDecompEig.cc
+++ b/dynare++/sylv/cc/SchurDecompEig.cc
@@ -5,11 +5,11 @@
 
 #include <vector>
 
-/* bubble diagonal 1-1, or 2-2 block from position 'from' to position
- * 'to'. If an eigenvalue cannot be swapped with its neighbour, the
- * neighbour is bubbled also in front. The method returns a new
- * position 'to', where the original block pointed by 'to' happens to
- * appear at the end. 'from' must be greater than 'to'.
+/* Bubble diagonal 1-1, or 2-2 block from position ‘from’ to position
+   ‘to’. If an eigenvalue cannot be swapped with its neighbour, the
+   neighbour is bubbled also in front. The method returns a new
+   position ‘to’, where the original block pointed by ‘to’ happens to
+   appear at the end. ‘from’ must be greater than ‘to’.
  */
 SchurDecompEig::diag_iter
 SchurDecompEig::bubbleEigen(diag_iter from, diag_iter to)
@@ -22,9 +22,9 @@ SchurDecompEig::bubbleEigen(diag_iter from, diag_iter to)
         ++to;
       else
         {
-          // bubble all eigenvalues from runm(incl.) to run(excl.),
-          // this includes either bubbling generated eigenvalues due
-          // to split, or an eigenvalue which couldn't be swapped
+          /* Bubble all eigenvalues from runm(incl.) to run(excl.),
+             this includes either bubbling generated eigenvalues due
+             to split, or an eigenvalue which couldn't be swapped */
           while (runm != run)
             {
               to = bubbleEigen(runm, to);
@@ -35,16 +35,16 @@ SchurDecompEig::bubbleEigen(diag_iter from, diag_iter to)
   return to;
 }
 
-/* this tries to swap two neighbouring eigenvalues, 'it' and '--it',
- * and returns 'itadd'. If the blocks can be swapped, new eigenvalues
- * can emerge due to possible 2-2 block splits. 'it' then points to
- * the last eigenvalue coming from block pointed by 'it' at the
- * begining, and 'itadd' points to the first. On swap failure, 'it' is
- * not changed, and 'itadd' points to previous eignevalue (which must
- * be moved backwards before). In either case, it is necessary to
- * resolve eigenvalues from 'itadd' to 'it', before the 'it' can be
- * resolved.
- * The success is signaled by returned true.
+/* This tries to swap two neighbouring eigenvalues, ‘it’ and ‘--it’,
+   and returns ‘itadd’. If the blocks can be swapped, new eigenvalues
+   can emerge due to possible 2-2 block splits. ‘it’ then points to
+   the last eigenvalue coming from block pointed by ‘it’ at the
+   begining, and ‘itadd’ points to the first. On swap failure, ‘it’ is
+   not changed, and ‘itadd’ points to previous eignevalue (which must
+   be moved backwards before). In either case, it is necessary to
+   resolve eigenvalues from ‘itadd’ to ‘it’, before the ‘it’ can be
+   resolved.
+   The success is signaled by returned true.
  */
 bool
 SchurDecompEig::tryToSwap(diag_iter &it, diag_iter &itadd)
@@ -66,10 +66,10 @@ SchurDecompEig::tryToSwap(diag_iter &it, diag_iter &itadd)
     {
       // swap successful
       getT().swapDiagLogically(itadd);
-      //check for 2-2 block splits
+      // check for 2-2 block splits
       getT().checkDiagConsistency(it);
       getT().checkDiagConsistency(itadd);
-      // and go back by 'it' in NEW eigenvalue set
+      // and go back by ‘it’ in NEW eigenvalue set
       --it;
       return true;
     }
diff --git a/dynare++/sylv/cc/SimilarityDecomp.cc b/dynare++/sylv/cc/SimilarityDecomp.cc
index 8d0e2df6b..1b7e0e37f 100644
--- a/dynare++/sylv/cc/SimilarityDecomp.cc
+++ b/dynare++/sylv/cc/SimilarityDecomp.cc
@@ -33,9 +33,9 @@ SimilarityDecomp::getXDim(diag_iter start, diag_iter end,
   rows = ei - si;
 }
 
-/* find solution of X for diagonal block given by start(incl.) and
- * end(excl.). If the solution cannot be found, or it is greater than
- * norm, X is not changed and flase is returned.
+/* Find solution of X for diagonal block given by start(incl.) and
+   end(excl.). If the solution cannot be found, or it is greater than
+   norm, X is not changed and flase is returned.
  */
 bool
 SimilarityDecomp::solveX(diag_iter start, diag_iter end,
@@ -68,7 +68,8 @@ SimilarityDecomp::solveX(diag_iter start, diag_iter end,
   return true;
 }
 
-/* multiply Q and invQ with (I -X; 0 I), and (I X; 0 I). This also sets X=-X. */
+/*                         ⎡ I -X ⎤     ⎡ I X ⎤
+  Multiply Q and invQ with ⎣ 0  I ⎦ and ⎣ 0 I ⎦. This also sets X=-X. */
 void
 SimilarityDecomp::updateTransform(diag_iter start, diag_iter end,
                                   GeneralMatrix &X)
@@ -127,21 +128,21 @@ SimilarityDecomp::diagonalize(double norm)
 void
 SimilarityDecomp::check(SylvParams &pars, const GeneralMatrix &m) const
 {
-  // M - Q*B*inv(Q)
+  // M - Q·B·Q⁻¹
   SqSylvMatrix c(getQ() * getB());
   c.multRight(getInvQ());
   c.add(-1.0, m);
   pars.f_err1 = c.getNorm1();
   pars.f_errI = c.getNormInf();
 
-  // I - Q*inv(Q)
+  // I - Q·Q⁻¹
   c.setUnit();
   c.mult(-1);
   c.multAndAdd(getQ(), getInvQ());
   pars.viv_err1 = c.getNorm1();
   pars.viv_errI = c.getNormInf();
 
-  // I - inv(Q)*Q
+  // I - Q⁻¹·Q
   c.setUnit();
   c.mult(-1);
   c.multAndAdd(getInvQ(), getQ());
diff --git a/dynare++/sylv/cc/SylvMatrix.cc b/dynare++/sylv/cc/SylvMatrix.cc
index 736295bfa..5a11477dd 100644
--- a/dynare++/sylv/cc/SylvMatrix.cc
+++ b/dynare++/sylv/cc/SylvMatrix.cc
@@ -42,9 +42,9 @@ SylvMatrix::multLeft(int zero_cols, const GeneralMatrix &a, const GeneralMatrix
       || rows != b.numRows() || cols != b.numCols())
     throw SYLV_MES_EXCEPTION("Wrong matrix dimensions for multLeft.");
 
-  // here we cannot call SylvMatrix::gemm since it would require
-  // another copy of (usually big) b (we are not able to do inplace
-  // submatrix of const GeneralMatrix)
+  /* Here we cannot call SylvMatrix::gemm() since it would require
+     another copy of (usually big) b (we are not able to do inplace
+     submatrix of const GeneralMatrix) */
   if (a.getLD() > 0 && ld > 0)
     {
       blas_int mm = a.numRows();
diff --git a/dynare++/sylv/cc/SylvMatrix.hh b/dynare++/sylv/cc/SylvMatrix.hh
index 8829d3401..387b20293 100644
--- a/dynare++/sylv/cc/SylvMatrix.hh
+++ b/dynare++/sylv/cc/SylvMatrix.hh
@@ -40,23 +40,23 @@ public:
   SylvMatrix &operator=(const SylvMatrix &m) = default;
   SylvMatrix &operator=(SylvMatrix &&m) = default;
 
-  /* this = |I 0|* this
-     |0 m|         */
+  /*        ⎡ I 0 ⎤
+     this = ⎣ 0 m ⎦ · this */
   void multLeftI(const SqSylvMatrix &m);
-  /* this = |I  0|* this
-     |0 m'|         */
+  /*        ⎡ I 0  ⎤
+     this = ⎣ 0 m' ⎦ · this */
   void multLeftITrans(const SqSylvMatrix &m);
-  /* this = |0 a|*b, so that |0 a| is square */
+  // this = [0 a]·b, so that [0 a] is square
   void multLeft(int zero_cols, const GeneralMatrix &a, const GeneralMatrix &b);
-  /* this = this * (m\otimes m..\otimes m) */
+  // this = this·(m⊗m⊗…⊗m)
   void multRightKron(const SqSylvMatrix &m, int order);
-  /* this = this * (m'\otimes m'..\otimes m') */
+  // this = this·(mᵀ⊗mᵀ⊗…⊗mᵀ)
   void multRightKronTrans(const SqSylvMatrix &m, int order);
-  /* this = P*this, x = P*x, where P is gauss transformation setting
-   * a given element to zero */
+  /* this = P·this, x = P·x, where P is gauss transformation setting
+     a given element to zero */
   void eliminateLeft(int row, int col, Vector &x);
-  /* this = this*P, x = P'*x, where P is gauss transformation setting
-   * a given element to zero */
+  /* this = this·P, x = Pᵀ·x, where P is gauss transformation setting
+     a given element to zero */
   void eliminateRight(int row, int col, Vector &x);
 };
 
@@ -84,14 +84,14 @@ public:
   }
   SqSylvMatrix &operator=(const SqSylvMatrix &m) = default;
   SqSylvMatrix &operator=(SqSylvMatrix &&m) = default;
-  /* x = (this \otimes this..\otimes this)*d */
+  // x = (this⊗this⊗…⊗this)·d
   void multVecKron(KronVector &x, const ConstKronVector &d) const;
-  /* x = (this' \otimes this'..\otimes this')*d */
+  // x = (thisᵀ⊗thisᵀ⊗…⊗thisᵀ)·d
   void multVecKronTrans(KronVector &x, const ConstKronVector &d) const;
-  /* a = inv(this)*a, b=inv(this)*b */
+  // a = this⁻¹·a, b=this⁻¹·b */
   void multInvLeft2(GeneralMatrix &a, GeneralMatrix &b,
                     double &rcond1, double &rcondinf) const;
-  /* this = I */
+  // this = I
   void setUnit();
 };
 
diff --git a/dynare++/sylv/cc/SylvesterSolver.hh b/dynare++/sylv/cc/SylvesterSolver.hh
index 622dc266b..be358c847 100644
--- a/dynare++/sylv/cc/SylvesterSolver.hh
+++ b/dynare++/sylv/cc/SylvesterSolver.hh
@@ -20,8 +20,8 @@ protected:
   const std::unique_ptr<const QuasiTriangular> matrixK;
   const std::unique_ptr<const QuasiTriangular> matrixF;
 private:
-  /* return true when it is more efficient to use QuasiTriangular
-   * than QuasiTriangularZero */
+  /* Return true when it is more efficient to use QuasiTriangular
+     than QuasiTriangularZero */
   static bool
   zeroPad(const SchurDecompZero &kdecomp)
   {
diff --git a/dynare++/sylv/cc/SymSchurDecomp.cc b/dynare++/sylv/cc/SymSchurDecomp.cc
index 82ca81dbc..a9ab2478c 100644
--- a/dynare++/sylv/cc/SymSchurDecomp.cc
+++ b/dynare++/sylv/cc/SymSchurDecomp.cc
@@ -78,8 +78,8 @@ SymSchurDecomp::getFactor(GeneralMatrix &f) const
     }
 }
 
-// LAPACK says that eigenvalues are ordered in ascending order, but we
-// do not rely her on it
+/* LAPACK says that eigenvalues are ordered in ascending order, but we
+   do not rely on it */
 bool
 SymSchurDecomp::isPositiveSemidefinite() const
 {
diff --git a/dynare++/sylv/cc/SymSchurDecomp.hh b/dynare++/sylv/cc/SymSchurDecomp.hh
index b9d7e1ecc..04deaf8ac 100644
--- a/dynare++/sylv/cc/SymSchurDecomp.hh
+++ b/dynare++/sylv/cc/SymSchurDecomp.hh
@@ -13,8 +13,8 @@ protected:
   Vector lambda;
   SqSylvMatrix q;
 public:
-  /** Calculates A = Q*Lambda*Q^T, where A is assummed to be
-   * symmetric and Lambda real diagonal, hence a vector. */
+  /* Computes the factorization A = Q·Λ·Qᵀ, where A is assummed to be
+     symmetric and Λ real diagonal, hence a vector. */
   SymSchurDecomp(const ConstGeneralMatrix &a);
   SymSchurDecomp(const SymSchurDecomp &ssd) = default;
   virtual ~SymSchurDecomp() = default;
@@ -28,15 +28,14 @@ public:
   {
     return q;
   }
-  /** Return factor F*F^T = A, raises and exception if A is not
-   * positive semidefinite, F must be square. */
+  /* Return factor F·Fᵀ = A, raises and exception if A is not
+     positive semidefinite, F must be square. */
   void getFactor(GeneralMatrix &f) const;
-  /** Returns true if A is positive semidefinite. */
+  // Returns true if A is positive semidefinite.
   bool isPositiveSemidefinite() const;
-  /** Correct definitness. This sets all eigenvalues between minus
-   * tolerance and zero to zero. */
+  /* Correct definitness. This sets all eigenvalues between minus
+     tolerance and zero to zero. */
   void correctDefinitness(double tol);
-
 };
 
 #endif
diff --git a/dynare++/sylv/cc/Vector.hh b/dynare++/sylv/cc/Vector.hh
index c5632f110..80ae718a3 100644
--- a/dynare++/sylv/cc/Vector.hh
+++ b/dynare++/sylv/cc/Vector.hh
@@ -5,9 +5,9 @@
 #ifndef VECTOR_H
 #define VECTOR_H
 
-/* NOTE! Vector and ConstVector have not common super class in order
- * to avoid running virtual method invokation mechanism. Some
- * members, and methods are thus duplicated */
+/* NOTE: Vector and ConstVector have not common super class in order
+   to avoid running virtual method invokation mechanism. Some
+   members, and methods are thus duplicated */
 
 #include <complex>
 
@@ -23,7 +23,7 @@ class Vector
   friend class ConstVector;
 protected:
   int len{0};
-  int s{1}; // stride (also called "skip" in some places)
+  int s{1}; // stride (also called “skip” in some places)
   double *data;
   bool destroy{true};
 public:
@@ -43,7 +43,7 @@ public:
     v.data = nullptr;
     v.destroy = false;
   }
-  // We don't want implict conversion from ConstVector, since it's expensive
+  // We don't want implict conversion from ConstVector, since it’s expensive
   explicit Vector(const ConstVector &v);
   Vector(double *d, int l)
     : len(l), data{d}, destroy{false}
@@ -90,10 +90,10 @@ public:
     return s;
   }
 
-  /** Exact equality. */
+  // Exact equality.
   bool operator==(const Vector &y) const;
   bool operator!=(const Vector &y) const;
-  /** Lexicographic ordering. */
+  // Lexicographic ordering.
   bool operator<(const Vector &y) const;
   bool operator<=(const Vector &y) const;
   bool operator>(const Vector &y) const;
@@ -125,8 +125,8 @@ public:
   void print() const;
 
   /* multiplies | alpha -beta1|           |b1|   |x1|
-     |             |\otimes I .|  | = |  |
-     | -beta2 alpha|           |b2|   |x2|
+                |             |\otimes I .|  | = |  |
+                | -beta2 alpha|           |b2|   |x2|
   */
   static void mult2(double alpha, double beta1, double beta2,
                     Vector &x1, Vector &x2,
diff --git a/dynare++/tl/cc/equivalence.cc b/dynare++/tl/cc/equivalence.cc
index 3c1a392e2..14080d500 100644
--- a/dynare++/tl/cc/equivalence.cc
+++ b/dynare++/tl/cc/equivalence.cc
@@ -39,7 +39,7 @@ OrdSequence::operator==(const OrdSequence &s) const
   return (i == length());
 }
 
-/* The first |add| adds a given integer to the class, the second
+/* The first add() adds a given integer to the class, the second
    iterates through a given sequence and adds everything found in the
    given class. */
 
@@ -65,7 +65,7 @@ OrdSequence::add(const OrdSequence &s)
     }
 }
 
-/* Answers |true| if a given number is in the class. */
+/* Answers true if a given number is in the class. */
 bool
 OrdSequence::has(int i) const
 {
@@ -183,7 +183,7 @@ Equivalence::findHaving(int i)
   return si;
 }
 
-/* Find $j$-th class for a given $j$. */
+/* Find j-th class for a given j. */
 
 Equivalence::const_seqit
 Equivalence::find(int j) const
@@ -228,7 +228,7 @@ Equivalence::insert(const OrdSequence &s)
 /* Trace the equivalence into the integer sequence. The classes are in
    some order (described earlier), and items within classes are ordered,
    so this implies, that the data can be linearized. This method
-   ``prints'' them to the sequence. We allow for tracing only a given
+   “prints” them to the sequence. We allow for tracing only a given
    number of classes from the beginning. */
 
 void
@@ -275,18 +275,16 @@ Equivalence::print(const std::string &prefix) const
     }
 }
 
-/* Here we construct a set of all equivalences over $n$-element
-   set. The construction proceeds as follows. We maintain a list of added
-   equivalences. At each iteration we pop front of the list, try to add
-   all parents of the popped equivalence. This action adds new
-   equivalences to the object and also to the added list. We finish the
-   iterations when the added list is empty.
+/* Here we construct a set of all equivalences over n-element set. The
+   construction proceeds as follows. We maintain a list of added equivalences.
+   At each iteration we pop front of the list, try to add all parents of the
+   popped equivalence. This action adds new equivalences to the object and also
+   to the added list. We finish the iterations when the added list is empty.
 
-   In the beginning we start with
-   $\{\{0\},\{1\},\ldots,\{n-1\}\}$. Adding of parents is an action which
-   for a given equivalence tries to glue all possible couples and checks
-   whether a new equivalence is already in the equivalence set. This is
-   not effective, but we will do the construction only ones.
+   In the beginning we start with { {0}, {1}, …, {n-1} }. Adding of parents is
+   an action which for a given equivalence tries to glue all possible couples
+   and checks whether a new equivalence is already in the equivalence set. This
+   is not effective, but we will do the construction only ones.
 
    In this way we breath-first search a lattice of all equivalences. Note
    that the lattice is modular, that is why the result of a construction
@@ -311,12 +309,12 @@ EquivalenceSet::EquivalenceSet(int num)
     equis.emplace_back(n, "");
 }
 
-/* This method is used in |addParents| and returns |true| if the object
+/* This method is used in addParents() and returns true if the object
    already has that equivalence. We trace list of equivalences in reverse
    order since equivalences are ordered in the list from the most
    primitive (nothing equivalent) to maximal (all is equivalent). Since
-   we will have much more results of |has| method as |true|, and
-   |operator==| between equivalences is quick if number of classes
+   we will have much more results of has() method as true, and
+   operator==() between equivalences is quick if number of classes
    differ, and in time we will compare with equivalences with less
    classes, then it is more efficient to trace the equivalences from less
    classes to more classes. hence the reverse order. */
@@ -334,7 +332,7 @@ EquivalenceSet::has(const Equivalence &e) const
 
 /* Responsibility of this methods is to try to glue all possible
    couples within a given equivalence and add those which are not in the
-   list yet. These are added also to the |added| list.
+   list yet. These are added also to the ‘added’ list.
 
    If number of classes is 2 or 1, we exit, because there is nothing to
    be added. */
@@ -373,14 +371,14 @@ EquivalenceSet::print(const std::string &prefix) const
     }
 }
 
-/* Construct the bundle. |nmax| is a maximum size of underlying set. */
+/* Construct the bundle. nmax is a maximum size of underlying set. */
 EquivalenceBundle::EquivalenceBundle(int nmax)
 {
   nmax = std::max(nmax, 1);
   generateUpTo(nmax);
 }
 
-/* Remember, that the first item is |EquivalenceSet(1)|. */
+/* Remember, that the first item is EquivalenceSet(1). */
 const EquivalenceSet &
 EquivalenceBundle::get(int n) const
 {
@@ -389,8 +387,8 @@ EquivalenceBundle::get(int n) const
   return bundle[n-1];
 }
 
-/* Get |curmax| which is a maximum size in the bundle, and generate for
-   all sizes from |curmax+1| up to |nmax|. */
+/* Get ‘curmax’ which is a maximum size in the bundle, and generate for
+   all sizes from curmax+1 up to nmax. */
 
 void
 EquivalenceBundle::generateUpTo(int nmax)
diff --git a/dynare++/tl/cc/equivalence.hh b/dynare++/tl/cc/equivalence.hh
index 45838a6a2..6df576be4 100644
--- a/dynare++/tl/cc/equivalence.hh
+++ b/dynare++/tl/cc/equivalence.hh
@@ -2,34 +2,30 @@
 
 // Equivalences.
 
-/* Here we define an equivalence of a set of integers $\{0, 1, \ldots,
-   k-1\}$. The purpose is clear, in the tensor library we often iterate
+/* Here we define an equivalence of a set of integers {0, 1, …, k-1}.
+   The purpose is clear, in the tensor library we often iterate
    through all equivalences and sum matrices. We need an abstraction for
    an equivalence class, equivalence and a set of all equivalences.
 
    The equivalence class (which is basically a set of integers) is here
    implemented as ordered integer sequence. The ordered sequence is not
-   implemented via |IntSequence|, but via |vector<int>| since we need
+   implemented via IntSequence, but via vector<int> since we need
    insertions. The equivalence is implemented as an ordered list of
    equivalence classes, and equivalence set is a list of equivalences.
 
    The ordering of the equivalence classes within an equivalence is very
-   important. For instance, if we iterate through equivalences for $k=5$
-   and pickup some equivalence class, say $\{\{0,4\},\{1,2\},\{3\}\}$, we
+   important. For instance, if we iterate through equivalences for k=5
+   and pickup some equivalence class, say { {0,4}, {1,2}, {3} }, we
    then evaluate something like:
-   $$\left[B_{y^2u^3}\right]_{\alpha_1\alpha_2\beta_1\beta_2\beta_3}=
-   \cdots+\left[g_{y^3}\right]_{\gamma_1\gamma_2\gamma_3}
-   \left[g_{yu}\right]^{\gamma_1}_{\alpha_1\beta_3}
-   \left[g_{yu}\right]^{\gamma_2}_{\alpha_2\beta_1}
-   \left[g_u\right]^{\gamma_3}_{\beta_2}+\cdots
-   $$
+
+   [B_y²u³]_α₁α₂β₁β₂β₃ = … + [f_z³]_γ₁γ₂γ₃ [z_yu]^γ₁_α₁β₃ [z_yu]^γ₂_α₂β₁ [zᵤ]^γ₃_β₂ + …
+
    If the tensors are unfolded, we can evaluate this expression as
-   $$g_{y^3}\cdot\left(g_{yu}\otimes g_{yu}\otimes g_{u}\right)\cdot P,$$
-   where $P$ is a suitable permutation of columns of the expressions,
-   which permutes them so that the index
-   $(\alpha_1,\beta_3,\alpha_2,\beta_1,\beta_2)$ would go to
-   $(\alpha_1,\alpha_2,\beta_1,\beta_2,\beta_3)$.
-   The permutation $P$ can be very ineffective (copying great amount of
+   f_z³·(z_yu ⊗ z_yu ⊗ zᵤ)·P, where P is a suitable permutation of columns of
+   the expressions, which permutes them so that the index
+   (α₁,β₃,α₂,β₁,β₂) would go to (α₁,α₂,β₁,β₂,β₃).
+
+   The permutation P can be very ineffective (copying great amount of
    small chunks of data) if the equivalence class ordering is chosen
    badly. However, we do not provide any heuristic minimizing a total
    time spent in all permutations. We choose an ordering which orders the
@@ -46,10 +42,10 @@
 #include <string>
 
 /* Here is the abstraction for an equivalence class. We implement it as
-   |vector<int>|. We have a constructor for empty class, copy
+   vector<int>. We have a constructor for empty class, copy
    constructor. What is important here is the ordering operator
-   |operator<| and methods for addition of an integer, and addition of
-   another sequence. Also we provide method |has| which returns true if a
+   operator<() and methods for addition of an integer, and addition of
+   another sequence. Also we provide method has() which returns true if a
    given integer is contained. */
 
 class OrdSequence
@@ -81,10 +77,10 @@ private:
 };
 
 /* Here is the abstraction for the equivalence. It is a list of
-   equivalence classes. Also we remember |n|, which is a size of
-   underlying set $\{0, 1, \ldots, n-1\}$.
+   equivalence classes. Also we remember n, which is a size of
+   underlying set {0, 1, …, n-1}.
 
-   Method |trace| ``prints'' the equivalence into the integer sequence. */
+   Method trace() “prints” the equivalence into the integer sequence. */
 
 class Permutation;
 class Equivalence
@@ -96,11 +92,11 @@ public:
   using const_seqit = std::list<OrdSequence>::const_iterator;
   using seqit = std::list<OrdSequence>::iterator;
 
-  // Constructs $\{\{0\},\{1\},\ldots,\{n-1\}\}$
+  // Constructs { {0}, {1}, …, {n-1} }
   explicit Equivalence(int num);
-  // Constructs $\{\{0,1,\ldots,n-1\}\}$
+  // Constructs { {0,1,…,n-1 } }
   Equivalence(int num, const std::string &dummy);
-  // Copy constructor plus gluing |i1| and |i2| in one class
+  // Copy constructor plus gluing i1 and i2 in one class
   Equivalence(const Equivalence &e, int i1, int i2);
 
   bool operator==(const Equivalence &e) const;
@@ -154,7 +150,7 @@ protected:
      given number or we can find an equivalence class of a given index within
      the ordering.
 
-     We have also an |insert| method which inserts a given class
+     We have also an insert() method which inserts a given class
      according to the class ordering. */
   const_seqit findHaving(int i) const;
   seqit findHaving(int i);
@@ -162,12 +158,12 @@ protected:
 
 };
 
-/* The |EquivalenceSet| is a list of equivalences. The unique
-   constructor constructs a set of all equivalences over $n$-element
+/* The EquivalenceSet is a list of equivalences. The unique
+   constructor constructs a set of all equivalences over an n-elements
    set. The equivalences are sorted in the list so that equivalences with
    fewer number of classes are in the end.
 
-   The two methods |has| and |addParents| are useful in the constructor. */
+   The two methods has() and addParents() are useful in the constructor. */
 
 class EquivalenceSet
 {
@@ -191,12 +187,12 @@ private:
   void addParents(const Equivalence &e, std::list<Equivalence> &added);
 };
 
-/* The equivalence bundle class only encapsulates |EquivalenceSet|s
+/* The equivalence bundle class only encapsulates EquivalenceSet·s
    from 1 up to a given number. It is able to retrieve the equivalence set
-   over $n$-element set for a given $n$, and also it can generate some more
+   over n-element set for a given n, and also it can generate some more
    sets on request.
 
-   It is fully responsible for storage needed for |EquivalenceSet|s. */
+   It is fully responsible for storage needed for EquivalenceSet·s. */
 
 class EquivalenceBundle
 {
diff --git a/dynare++/tl/cc/fine_container.cc b/dynare++/tl/cc/fine_container.cc
index 75ccacba4..50c1a7b1c 100644
--- a/dynare++/tl/cc/fine_container.cc
+++ b/dynare++/tl/cc/fine_container.cc
@@ -5,7 +5,7 @@
 #include <cmath>
 
 /* Here we construct the vector of new sizes of containers (before
-   |nc|) and copy all remaining sizes behind |nc|. */
+   nc) and copy all remaining sizes behind nc. */
 
 SizeRefinement::SizeRefinement(const IntSequence &s, int nc, int max)
 {
diff --git a/dynare++/tl/cc/fine_container.hh b/dynare++/tl/cc/fine_container.hh
index 014b3a6b4..33be347db 100644
--- a/dynare++/tl/cc/fine_container.hh
+++ b/dynare++/tl/cc/fine_container.hh
@@ -12,14 +12,14 @@
    of each stack in the refined container. The resulting object is stack
    container, so everything works seamlessly.
 
-   We define here a class for refinement of sizes |SizeRefinement|, this
+   We define here a class for refinement of sizes SizeRefinement, this
    is purely an auxiliary class allowing us to write a code more
-   concisely. The main class of this file is |FineContainer|, which
-   corresponds to refining. The two more classes |FoldedFineContainer|
-   and |UnfoldedFineContainer| are its specializations.
+   concisely. The main class of this file is FineContainer, which
+   corresponds to refining. The two more classes FoldedFineContainer
+   and UnfoldedFineContainer are its specializations.
 
    NOTE: This code was implemented with a hope that it will help to cut
-   down memory allocations during the Faa Di Bruno formula
+   down memory allocations during the Faà Di Bruno formula
    evaluation. However, it seems that this needs to be accompanied with a
    similar thing for tensor multidimensional index. Thus, the abstraction
    is not currently used, but it might be useful in future. */
@@ -32,9 +32,9 @@
 #include <vector>
 #include <memory>
 
-/* This class splits the first |nc| elements of the given sequence |s|
-   to a sequence not having items greater than given |max|. The remaining
-   elements (those behind |nc|) are left untouched. It also remembers the
+/* This class splits the first nc elements of the given sequence s
+   to a sequence not having items greater than given max. The remaining
+   elements (those behind nc) are left untouched. It also remembers the
    mapping, i.e. for a given index in a new sequence, it is able to
    return a corresponding index in old sequence. */
 
@@ -68,7 +68,7 @@ public:
 };
 
 /* This main class of this class refines a given stack container, and
-   inherits from the stack container. It also defines the |getType|
+   inherits from the stack container. It also defines the getType()
    method, which returns a type for a given stack as the type of the
    corresponding (old) stack of the former stack container. */
 
@@ -82,19 +82,19 @@ protected:
   std::vector<std::unique_ptr<_Ctype>> ref_conts;
   const _Stype &stack_cont;
 public:
-  /* Here we construct the |SizeRefinement| and allocate space for the
+  /* Here we construct the SizeRefinement and allocate space for the
      refined containers. Then, the containers are created and put to
-     |conts| array. Note that the containers do not claim any further
+     conts array. Note that the containers do not claim any further
      space, since all the tensors of the created containers are in-place
      submatrices.
 
-     Here we use a dirty trick of converting |const| pointer to non-|const|
+     Here we use a dirty trick of converting const pointer to non-const
      pointer and passing it to a subtensor container constructor. The
-     containers are stored in |ref_conts| and then in |conts| from
-     |StackContainer|. However, this is safe since neither |ref_conts| nor
-     |conts| are used in non-|const| contexts. For example,
-     |StackContainer| has only a |const| method to return a member of
-     |conts|. */
+     containers are stored in ref_conts and then in conts from
+     StackContainer. However, this is safe since neither ref_conts nor
+     conts are used in non-const contexts. For example,
+     StackContainer has only a const method to return a member of
+     conts. */
 
   FineContainer(const _Stype &sc, int max)
     : SizeRefinement(sc.getStackSizes(), sc.numConts(), max),
@@ -130,7 +130,7 @@ public:
 
 };
 
-/* Here is |FineContainer| specialization for folded tensors. */
+/* Here is FineContainer specialization for folded tensors. */
 class FoldedFineContainer : public FineContainer<FGSTensor>, public FoldedStackContainer
 {
 public:
@@ -140,7 +140,7 @@ public:
   }
 };
 
-/* Here is |FineContainer| specialization for unfolded tensors. */
+/* Here is FineContainer specialization for unfolded tensors. */
 class UnfoldedFineContainer : public FineContainer<UGSTensor>, public UnfoldedStackContainer
 {
 public:
diff --git a/dynare++/tl/cc/fs_tensor.cc b/dynare++/tl/cc/fs_tensor.cc
index 24e097b85..7c207fa07 100644
--- a/dynare++/tl/cc/fs_tensor.cc
+++ b/dynare++/tl/cc/fs_tensor.cc
@@ -8,13 +8,13 @@
 #include "pascal_triangle.hh"
 
 /* This constructs a fully symmetric tensor as given by the contraction:
-   $$\left[g_{y^n}\right]_{\alpha_1\ldots\alpha_n}=
-   \left[t_{y^{n+1}}\right]_{\alpha_1\ldots\alpha_n\beta}[x]^\beta$$
 
-   We go through all columns of output tensor $[g]$ and for each column
+   [g_yⁿ]_α₁…αₙ = [t_yⁿ⁺¹]_α₁…αₙβ [x]^β
+
+   We go through all columns of output tensor [g] and for each column
    we cycle through all variables, insert a variable to the column
-   coordinates obtaining a column of tensor $[t]$. the column is multiplied
-   by an appropriate item of |x| and added to the column of $[g]$ tensor. */
+   coordinates obtaining a column of tensor [t]. The column is multiplied
+   by an appropriate item of x and added to the column of [g] tensor. */
 
 FFSTensor::FFSTensor(const FFSTensor &t, const ConstVector &x)
   : FTensor(indor::along_col, IntSequence(t.dimen()-1, t.nvar()),
@@ -38,8 +38,10 @@ FFSTensor::FFSTensor(const FFSTensor &t, const ConstVector &x)
 }
 
 /* This returns number of indices for folded tensor with full
-   symmetry. Let $n$ be a number of variables |nvar| and $d$ the
-   dimension |dim|. Then the number of indices is $\pmatrix{n+d-1\cr d}$. */
+   symmetry. Let n be a number of variables and d the
+                                                ⎛n+d-1⎞
+   dimension dim. Then the number of indices is ⎝  d  ⎠.
+ */
 
 int
 FFSTensor::calcMaxOffset(int nvar, int d)
@@ -91,7 +93,7 @@ FFSTensor::unfold() const
 }
 
 /* Incrementing is easy. We have to increment by calling static method
-   |UTensor::increment| first. In this way, we have coordinates of
+   UTensor::increment() first. In this way, we have coordinates of
    unfolded tensor. Then we have to skip to the closest folded index
    which corresponds to monotonizeing the integer sequence. */
 
@@ -105,7 +107,7 @@ FFSTensor::increment(IntSequence &v) const
   v.monotone();
 }
 
-/* Decrement calls static |FTensor::decrement|. */
+/* Decrement calls static FTensor::decrement(). */
 
 void
 FFSTensor::decrement(IntSequence &v) const
@@ -143,7 +145,7 @@ FFSTensor::addSubTensor(const FGSTensor &t)
   TL_RAISE_IF(nvar() != t.getDims().getNVS().sum(),
               "Wrong nvs for FFSTensor::addSubTensor");
 
-  // set shift for |addSubTensor|
+  // set shift for addSubTensor()
   /* Code shared with UFSTensor::addSubTensor() */
   IntSequence shift_pre(t.getSym().num(), 0);
   for (int i = 1; i < t.getSym().num(); i++)
@@ -160,8 +162,8 @@ FFSTensor::addSubTensor(const FGSTensor &t)
     }
 }
 
-// |UFSTensor| contraction constructor
-/* This is a bit more straightforward than |@<|FFSTensor| contraction constructor@>|.
+// UFSTensor contraction constructor
+/* This is a bit more straightforward than FFSTensor contraction constructor.
    We do not add column by column but we do it by submatrices due to
    regularity of the unfolded tensor. */
 
@@ -186,7 +188,7 @@ UFSTensor::UFSTensor(const UFSTensor &t, const ConstVector &x)
 }
 
 /* Here we convert folded full symmetry tensor to unfolded. We copy all
-   columns of folded tensor, and then call |unfoldData()|. */
+   columns of folded tensor, and then call unfoldData(). */
 
 UFSTensor::UFSTensor(const FFSTensor &ft)
   : UTensor(indor::along_col, IntSequence(ft.dimen(), ft.nvar()),
@@ -207,8 +209,8 @@ UFSTensor::fold() const
   return std::make_unique<FFSTensor>(*this);
 }
 
-// |UFSTensor| increment and decrement
-/* Here we just call |UTensor| respective static methods. */
+// UFSTensor increment and decrement
+/* Here we just call UTensor respective static methods. */
 void
 UFSTensor::increment(IntSequence &v) const
 {
@@ -236,7 +238,7 @@ UFSTensor::getOffset(const IntSequence &v) const
   return UTensor::getOffset(v, nv);
 }
 
-/* This is very similar to |@<|FFSTensor::addSubTensor| code@>|. The
+/* This is very similar to FFSTensor::addSubTensor(). The
    only difference is the addition. We go through all columns in the full
    symmetry tensor and cancel the shift. If the coordinates after the
    cancellation are positive, we find the column in the general symmetry
@@ -250,7 +252,7 @@ UFSTensor::addSubTensor(const UGSTensor &t)
   TL_RAISE_IF(nvar() != t.getDims().getNVS().sum(),
               "Wrong nvs for UFSTensor::addSubTensor");
 
-  // set shift for |addSubTensor|
+  // set shift for addSubTensor()
   /* Code shared with FFSTensor::addSubTensor() */
   IntSequence shift_pre(t.getSym().num(), 0);
   for (int i = 1; i < t.getSym().num(); i++)
diff --git a/dynare++/tl/cc/fs_tensor.hh b/dynare++/tl/cc/fs_tensor.hh
index b4aaacb4f..7fc926673 100644
--- a/dynare++/tl/cc/fs_tensor.hh
+++ b/dynare++/tl/cc/fs_tensor.hh
@@ -19,21 +19,21 @@ class FRSingleTensor;
 class FSSparseTensor;
 
 /* Folded tensor with full symmetry maintains only information about
-   number of symmetrical variables |nv|. Further, we implement what is
-   left from the super class |FTensor|.
+   number of symmetrical variables nv. Further, we implement what is
+   left from the super class FTensor.
 
-   We implement |getOffset| which should be used with care since
+   We implement getOffset() which should be used with care since
    its complexity.
 
    We implement a method adding a given general symmetry tensor to the
    full symmetry tensor supposing the variables of the general symmetry
    tensor are stacked giving only one variable of the full symmetry
-   tensor. For instance, if $x=[y^T, u^T]^T$, then we can add tensor
-   $\left[g_{y^2u}\right]$ to tensor $g_{x^3}$. This is done in method
-   |addSubTensor|. Consult |@<|FGSTensor| class declaration@>| to know
+   tensor. For instance, if x=[yᵀ, uᵀ]ᵀ, then we can add tensor
+   [g_y²u] to tensor [g_x³]. This is done in method
+   addSubTensor(). Consult FGSTensor class declaration to know
    what is general symmetry tensor.
 
-   Note that the past-the-end index is of the form $(nv, \ldots, nv)$, because
+   Note that the past-the-end index is of the form (nv,…,nv), because
    of the specific implementation of FFSTensor::increment().
 */
 
@@ -53,10 +53,11 @@ public:
   }
 
   /* Constructs a tensor by one-dimensional contraction from the higher
-     dimensional tensor |t|. This is, it constructs a tensor
-     $$\left[g_{y^n}\right]_{\alpha_1\ldots\alpha_n}=
-     \left[t_{y^{n+1}}\right]_{\alpha_1\ldots\alpha_n\beta}[x]^\beta$$ See the
-     implementation for details. */
+     dimensional tensor t. This is, it constructs a tensor
+
+     [g_yⁿ]_α₁…αₙ = [t_yⁿ⁺¹]_α₁…αₙβ [x]^β
+
+     See the implementation for details. */
   FFSTensor(const FFSTensor &t, const ConstVector &x);
 
   /* Converts from sparse tensor (which is fully symmetric and folded by
@@ -95,7 +96,7 @@ public:
 };
 
 /* Unfolded fully symmetric tensor is almost the same in structure as
-   |FFSTensor|, but the method |unfoldData|. It takes columns which also
+   FFSTensor, but the method unfoldData(). It takes columns which also
    exist in folded version and copies them to all their symmetrical
    locations. This is useful when constructing unfolded tensor from
    folded one. */
diff --git a/dynare++/tl/cc/ps_tensor.hh b/dynare++/tl/cc/ps_tensor.hh
index e8d50dd59..9d4472dbf 100644
--- a/dynare++/tl/cc/ps_tensor.hh
+++ b/dynare++/tl/cc/ps_tensor.hh
@@ -269,11 +269,11 @@ protected:
 };
 
 /* Here we define an abstraction of the permuted symmetry folded
-   tensor. It is needed in context of the Faa Di Bruno formula for folded
+   tensor. It is needed in context of the Faà Di Bruno formula for folded
    stack container multiplied with container of dense folded tensors, or
    multiplied by one full symmetry sparse tensor.
 
-   For example, if we perform the Faa Di Bruno for $F=f(z)$, where
+   For example, if we perform the Faà Di Bruno for $F=f(z)$, where
    $z=[g(x,y,u,v), h(x,y,u), x, y]^T$, we get for one concrete
    equivalence:
    $$
diff --git a/dynare++/tl/cc/pyramid_prod.hh b/dynare++/tl/cc/pyramid_prod.hh
index 9038ebaa9..fe37671ac 100644
--- a/dynare++/tl/cc/pyramid_prod.hh
+++ b/dynare++/tl/cc/pyramid_prod.hh
@@ -2,8 +2,8 @@
 
 // Multiplying tensor columns.
 
-/* In here, we implement the Faa Di Bruno for folded
-   tensors. Recall, that one step of the Faa Di Bruno is a formula:
+/* In here, we implement the Faà Di Bruno for folded
+   tensors. Recall, that one step of the Faà Di Bruno is a formula:
    $$\left[B_{s^k}\right]_{\alpha_1\ldots\alpha_k}=
    [h_{y^l}]_{\gamma_1\ldots\gamma_l}
    \prod_{m=1}^l\left[g_{s^{\vert c_m\vert}}\right]^{\gamma_m}_{c_m(\alpha)}
diff --git a/dynare++/tl/cc/stack_container.cc b/dynare++/tl/cc/stack_container.cc
index 17b89665e..95ac42c24 100644
--- a/dynare++/tl/cc/stack_container.cc
+++ b/dynare++/tl/cc/stack_container.cc
@@ -28,7 +28,7 @@ FoldedStackContainer::multAndAdd(const FSSparseTensor &t,
 }
 
 // |FoldedStackContainer::multAndAdd| dense code
-/* Here we perform the Faa Di Bruno step for a given dimension |dim|, and for
+/* Here we perform the Faà Di Bruno step for a given dimension |dim|, and for
    the dense fully symmetric tensor which is scattered in the container
    of general symmetric tensors. The implementation is pretty the same as
    |@<|UnfoldedStackContainer::multAndAdd| dense code@>|. */
diff --git a/dynare++/tl/cc/stack_container.hh b/dynare++/tl/cc/stack_container.hh
index 28e93401a..d4dcee5db 100644
--- a/dynare++/tl/cc/stack_container.hh
+++ b/dynare++/tl/cc/stack_container.hh
@@ -3,28 +3,35 @@
 // Stack of containers.
 
 /* Here we develop abstractions for stacked containers of tensors. For
-   instance, in perturbation methods for DSGE we need function
-   $$z(y,u,u',\sigma)=\left[\matrix{G(y,u,u',\sigma)\cr g(y,u,\sigma)\cr y\cr u}\right]$$
-   and we need to calculate one step of Faa Di Bruno formula
-   $$\left[B_{s^k}\right]_{\alpha_1\ldots\alpha_l}=\left[f_{z^l}\right]_{\beta_1\ldots\beta_l}
-   \sum_{c\in M_{l,k}}\prod_{m=1}^l\left[z_{s^k(c_m)}\right]^{\beta_m}_{c_m(\alpha)}$$
-   where we have containers for derivatives of $G$ and $g$.
+   instance, in perturbation methods for DSGE models, we need the function:
+
+                  ⎡ G(y*,u,u′,σ) ⎤
+   z(y*,u,u′,σ) = ⎢  g(y*,u,σ)   ⎥
+                  ⎢      y*      ⎥
+                  ⎣      u       ⎦
+
+   and we need to calculate one step of Faà Di Bruno formula:
+
+                                            ₗ
+   [B_sᵏ]_α₁,…,αₗ = [f_zˡ]_β₁,…,βₗ    ∑     ∏  [z_(s^|cₘ|)]_cₘ(α)^βₘ
+                                   c∈ℳₗ,ₖ ᵐ⁼¹
+
+   where we have containers for derivatives of G and g.
 
    The main purpose of this file is to define abstractions for stack of
-   containers and possibly raw variables, and code |multAndAdd| method
-   calculating (one step of) the Faa Di Bruno formula for folded and
-   unfolded tensors. Note also, that tensors $\left[f_{z^l}\right]$ are
-   sparse.
+   containers and possibly raw variables, and code multAndAdd() method
+   calculating (one step of) the Faà Di Bruno formula for folded and
+   unfolded tensors. Note also, that tensors [f_zˡ] are sparse.
 
    The abstractions are built as follows. At the top, there is an
    interface describing stack of columns. It contains pure virtual
    methods needed for manipulating the container stack. For technical
    reasons it is a template. Both versions (folded, and unfolded) provide
-   all interface necessary for implementation of |multAndAdd|. The second
+   all interface necessary for implementation of multAndAdd(). The second
    way of inheritance is first general implementation of the interface
-   |StackContainer|, and then specific (|ZContainer| for our specific
-   $z$). The only method which is virtual also after |StackContainer| is
-   |getType|, which is implemented in the specialization and determines
+   StackContainer, and then specific (ZContainer for our specific z).
+   The only method which is virtual also after StackContainer is
+   getType(), which is implemented in the specialization and determines
    behaviour of the stack. The complete classes are obtained by
    inheriting from the both branches, as it is drawn below:
 
@@ -53,7 +60,7 @@
    {|UnfoldedStackContainer|}{|ZContainer<UGSTensor>|}{|UnfoldedZContainer|}
    }
 
-   We have also two supporting classes |StackProduct| and |KronProdStack|
+   We have also two supporting classes StackProduct and KronProdStack
    and a number of worker classes used as threads. */
 
 #ifndef STACK_CONTAINER_H
@@ -85,8 +92,8 @@
 
    Method |createPackedColumn| returns a vector of stack derivatives with
    respect to the given symmetry and of the given column, where all zeros
-   from zero types, or unit matrices are deleted. See {\tt
-   kron\_prod2.hweb} for explanation. */
+   from zero types, or unit matrices are deleted. See kron_prod.hh for an
+   explanation. */
 
 template <class _Ttype>
 class StackContainerInterface
diff --git a/dynare++/tl/cc/t_container.cc b/dynare++/tl/cc/t_container.cc
index 9d9c44764..ad14e1773 100644
--- a/dynare++/tl/cc/t_container.cc
+++ b/dynare++/tl/cc/t_container.cc
@@ -55,7 +55,7 @@ FGSContainer::FGSContainer(const UGSContainer &c)
 }
 
 // |FGSContainer::multAndAdd| folded code
-/* Here we perform one step of the Faa Di Bruno operation. We call the
+/* Here we perform one step of the Faà Di Bruno operation. We call the
    |multAndAdd| for unfolded tensor. */
 void
 FGSContainer::multAndAdd(const FGSTensor &t, FGSTensor &out) const
diff --git a/dynare++/tl/cc/t_container.hh b/dynare++/tl/cc/t_container.hh
index 53bf58efd..71bea1f6f 100644
--- a/dynare++/tl/cc/t_container.hh
+++ b/dynare++/tl/cc/t_container.hh
@@ -3,7 +3,7 @@
 // Tensor containers.
 
 /* One of primary purposes of the tensor library is to perform one step
-   of the Faa Di Bruno formula:
+   of the Faà Di Bruno formula:
    $$\left[B_{s^k}\right]_{\alpha_1\ldots\alpha_k}=
    [h_{y^l}]_{\gamma_1\ldots\gamma_l}\sum_{c\in M_{l,k}}
    \prod_{m=1}^l\left[g_{s^{\vert c_m\vert}}\right]^{\gamma_m}_{c_m(\alpha)}
@@ -30,14 +30,14 @@
    most one tensor.
 
    The class has two purposes: The first is to provide storage (insert
-   and retrieve). The second is to perform the above step of Faa Di Bruno. This is
+   and retrieve). The second is to perform the above step of Faà Di Bruno. This is
    going through all equivalences with $l$ classes, perform the tensor
    product and add to the result.
 
    We define a template class |TensorContainer|. From different
    instantiations of the template class we will inherit to create concrete
    classes, for example container of unfolded general symmetric
-   tensors. The one step of the Faa Di Bruno (we call it |multAndAdd|) is
+   tensors. The one step of the Faà Di Bruno (we call it |multAndAdd|) is
    implemented in the concrete subclasses, because the implementation
    depends on storage. Note even, that |multAndAdd| has not a template
    common declaration. This is because sparse tensor $h$ is multiplied by
diff --git a/dynare++/utils/cc/pascal_triangle.cc b/dynare++/utils/cc/pascal_triangle.cc
index 90ed80cc9..67d987ee4 100644
--- a/dynare++/utils/cc/pascal_triangle.cc
+++ b/dynare++/utils/cc/pascal_triangle.cc
@@ -11,10 +11,13 @@ PascalRow::setFromPrevious(const PascalRow &prev)
   prolong(prev);
 }
 
-/** This prolongs the PascalRow. If it is empty, we set the first item
- * to k+1, which is noverk(k+1,k) which is the second item in the real
- * pascal row, which starts from noverk(k,k)=1. Then we calculate
- * other items from the provided row which must be the one with k-1.*/
+/* This prolongs the PascalRow. If it is empty, we set the first item
+                    ⎛k+1⎞
+   to k+1, which is ⎝ k ⎠, which is the second item in the real
+                                 ⎛k⎞
+   pascal row, which starts from ⎝k⎠=1. Then we calculate
+   other items from the provided row which must be the one with k-1.
+ */
 void
 PascalRow::prolong(const PascalRow &prev)
 {
diff --git a/dynare++/utils/cc/pascal_triangle.hh b/dynare++/utils/cc/pascal_triangle.hh
index 8b9e2fe08..34f80b1da 100644
--- a/dynare++/utils/cc/pascal_triangle.hh
+++ b/dynare++/utils/cc/pascal_triangle.hh
@@ -24,6 +24,9 @@ public:
 namespace PascalTriangle
 {
   void ensure(int n, int k);
+  /*         ⎛n⎞
+    Computes ⎝k⎠, hence the function name (“n over k”).
+  */
   int noverk(int n, int k);
   void print();
 };