t_polynomial.hh

/*
 * Copyright © 2004 Ondra Kamenik
 * Copyright © 2019 Dynare Team
 *
 * This file is part of Dynare.
 *
 * Dynare is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * Dynare is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 */

// Tensor polynomial evaluation.

/* We need to evaluate a tensor polynomial of the form:
                                                                   ₙ
   [g_x]_α₁ [x]^α₁ + [g_x²]_α₁α₂ [x]^α₁ [x]^α₂ + … + [g_xⁿ]_α₁…αₙ  ∏  [x]^αᵢ
                                                                  ⁱ⁼¹
   where x is a column vector.

   We have basically two options. The first is to use the formula above,
   the second is to use a Horner-like formula:

    ⎡ ⎡                                                   ⎤ ⎤
    ⎢ ⎢⎡                                  ⎤               ⎥ ⎥
    ⎢…⎢⎢[g_xⁿ⁻¹] + [g_xⁿ]_α₁…αₙ₋₁αₙ [x]^αₙ⎥       [x]^αₙ₋₁⎥…⎥  [x]^α₁
    ⎢ ⎢⎣                                  ⎦α₁…αₙ₋₁        ⎥ ⎥
    ⎣ ⎣                                                   ⎦ ⎦α₁

   Alternatively, we can put the the polynomial into a more compact form

               ₙ ⎡1⎤αᵢ
    [G]_α₁…αₙ  ∏ ⎢ ⎥
              ⁱ⁼¹⎣x⎦

   Then the polynomial evaluation becomes just a matrix multiplication of the
   vector power.

   Here we define the tensor polynomial as a container of full symmetry
   tensors and add an evaluation methods. We have two sorts of
   containers, folded and unfolded. For each type we declare two methods
   implementing the above formulas. We define classes for the
   compactification of the polynomial. The class derives from the tensor
   and has a eval method. */

#include "t_container.hh"
#include "fs_tensor.hh"
#include "rfs_tensor.hh"
#include "tl_static.hh"
#include "pascal_triangle.hh"

#include <memory>

/* Just to make the code nicer, we implement a Kronecker power of a
   vector encapsulated in the following class. It has getNext() method
   which returns either folded or unfolded row-oriented single column
   Kronecker power of the vector according to the type of a dummy
   argument. This allows us to use the type dependent code in templates
   below.

   The implementation of the Kronecker power is that we maintain the last
   unfolded power. If unfolded getNext() is called, we Kronecker multiply