Dynare++: make multinomial coeffs computation a method of IntSequence

Also improve on the comments.

dynare++/tl/cc/int_sequence.cc

@@ -3,6 +3,7 @@
 #include "int_sequence.hh"
 #include "symmetry.hh"
 #include "tl_exception.hh"
+#include "pascal_triangle.hh"
 
 #include <iostream>
 #include <limits>
@@ -262,3 +263,36 @@ IntSequence::print() const
     std::cout << operator[](i) << ' ';
   std::cout << ']' << std::endl;
 }
+
+/* Here we calculate the multinomial coefficients
+    $\left(\matrix{a\cr b_1,\ldots,b_n}\right)$, where $a=b_1+\ldots+b_n$.
+
+   See:
+    https://en.wikipedia.org/wiki/Binomial_coefficient#Generalization_to_multinomials
+    https://en.wikipedia.org/wiki/Multinomial_theorem
+
+   For n=1, the coefficient is equal to 1.
+   For n=2, the multinomial coeffs correspond to the binomial coeffs, i.e. the binomial
+    (a; b) is equal to the multinomial (a; b,a-b).
+   For n>=3, we have the identity
+   $$\left(\matrix{a\cr b_1,\ldots,b_n}\right)=\left(\matrix{b_1+b_2\cr b_1}\right)\cdot
+   \left(\matrix{a\cr b_1+b_2,b_3,\ldots,b_n}\right)$$ (where the first factor
+   on the right hand side is to be interpreted as a binomial coefficient)
+
+   This number is exactly a number of unfolded indices corresponding to one
+   folded index, where the sequence $b_1,\ldots,b_n$ is the symmetry of the
+   index. This can be easily seen if the multinomial coefficient is interpreted
+   as the number of unique permutations of a word, where $a$ is the length of
+   the word, $n$ is the number of distinct letters, and the $b_i$ are the
+   number of repetitions of each letter. For example, for a symmetry of the
+   form $y^4 u^2 v^3$, we want to compute the number of permutations of the word
+   $yyyyuuvvv$. This is equal to the multinomial coefficient (9; 4,2,3). */
+
+int
+IntSequence::noverseq()
+{
+  if (size() == 0 || size() == 1)
+    return 1;
+  data[1] += data[0];
+  return PascalTriangle::noverk(data[1], data[0]) * IntSequence(*this, 1, size()).noverseq();
+}

dynare++/tl/cc/int_sequence.hh

@@ -158,6 +158,10 @@ public:
   bool isConstant() const;
   bool isSorted() const;
   void print() const;
+  // Compute multinomial coefficient (sum(this); this)
+  /* Warning: this operation is destructive; make a copy if you want to keep
+     the original sequence */
+  int noverseq();
 };
 
 #endif

dynare++/tl/cc/sparse_tensor.cc

@@ -82,8 +82,7 @@ SparseTensor::getUnfoldIndexFillFactor() const
   while (start_col != m.end())
     {
       const IntSequence &key = (*start_col).first;
-      Symmetry s(key);
+      cnt += Symmetry(key).noverseq();
-      cnt += Tensor::noverseq(s);
       start_col = m.upper_bound(key);
     }
 

dynare++/tl/cc/tensor.cc

@@ -5,26 +5,6 @@
 #include "tl_static.hh"
 #include "pascal_triangle.hh"
 
-// |Tensor::noverseq_ip| static method
-/* Here we calculate a generalized combination number
-   $\left(\matrix{a\cr b_1,\ldots,b_n}\right)$, where $a=b_1+\ldots+
-   b_n$. We use the identity
-   $$\left(\matrix{a\cr b_1,\ldots,b_n}\right)=\left(\matrix{b_1+b_2\cr b_1}\right)\cdot
-   \left(\matrix{a\cr b_1+b_2,b_3,\ldots,b_n}\right)$$
-
-   This number is exactly a number of unfolded indices corresponding to
-   one folded index, where the sequence $b_1,\ldots,b_n$ is the symmetry
-   of the index. */
-
-int
-Tensor::noverseq_ip(IntSequence &s)
-{
-  if (s.size() == 0 || s.size() == 1)
-    return 1;
-  s[1] += s[0];
-  return PascalTriangle::noverk(s[1], s[0]) * noverseq(IntSequence(s, 1, s.size()));
-}
-
 /* Here we increment a given sequence within full symmetry given by
    |nv|, which is number of variables in each dimension. The underlying
    tensor is unfolded, so we increase the rightmost by one, and if it is

dynare++/tl/cc/tensor.hh

@@ -217,15 +217,6 @@ public:
   {
     return in_end;
   }
-
-  static int
-  noverseq(const IntSequence &s)
-  {
-    IntSequence seq(s);
-    return noverseq_ip((IntSequence &) s);
-  }
-private:
-  static int noverseq_ip(IntSequence &s);
 };
 
 /* Here is an abstraction for unfolded tensor. We provide a pure