// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
// SPDX-FileCopyrightText: The Eigen Authors
// SPDX-License-Identifier: MPL-2.0

// References:
//  [1] P. J. Davis, "Circulant Matrices", Wiley, 1979. Diagonalization of a
//      circulant matrix by the DFT and the closed-form eigenstructure used by
//      eigenvalues()/eigenvectors(); the SVD and pseudo-inverse below follow
//      from it by taking moduli/phases of the eigenvalues.
//  [2] R. M. Gray, "Toeplitz and Circulant Matrices: A Review", Foundations and
//      Trends in Communications and Information Theory, 2(3), 2006.
//  [3] G. H. Golub and C. F. Van Loan, "Matrix Computations", 4th ed., Johns
//      Hopkins University Press, 2013, chapter 4.8 (circulant systems and
//      FFT-based products) and chapter 5.4 (numerical rank conventions).
//  [4] J. J. Dongarra, J. R. Bunch, C. B. Moler and G. W. Stewart, "LINPACK
//      Users' Guide", SIAM, 1979. determinant()'s balanced accumulation follows
//      the convention of its xGEDI routines, which return determinants as a
//      (fraction, exponent) pair to avoid spurious overflow/underflow.
//  [5] P. H. Sterbenz, "Floating-Point Computation", Prentice-Hall, 1974.
//      Scaling by a power of two is exact, the property the balanced
//      accumulation relies on.

#ifndef EIGEN_STRUCTURED_CIRCULANT_H
#define EIGEN_STRUCTURED_CIRCULANT_H

// IWYU pragma: private
#include "./InternalHeaderCheck.h"

namespace Eigen {

template <typename Scalar_, int Size_ = Dynamic>
class Circulant;

namespace internal {

template <typename Scalar_, int Size_>
struct traits<Circulant<Scalar_, Size_>> {
  using Scalar = Scalar_;
  using StorageKind = Dense;
  using XprKind = MatrixXpr;
  using StorageIndex = int;
  static constexpr int RowsAtCompileTime = Size_;
  static constexpr int ColsAtCompileTime = Size_;
  static constexpr int MaxRowsAtCompileTime = Size_;
  static constexpr int MaxColsAtCompileTime = Size_;
  // Deliberately no NestByRefBit: transpose(), conjugate() and adjoint() return
  // owning temporaries, so Product must nest the operator by value for a
  // delayed-evaluated product expression to keep its left factor alive. The copy
  // is O(n), negligible against the O(n log n) product evaluation.
  static constexpr int Flags = 0;
};

template <typename Scalar_, int Size_>
struct evaluator_traits<Circulant<Scalar_, Size_>> {
  using Kind = IndexBased;
  using Shape = StructuredShape;
};

}  // namespace internal

/** \ingroup StructuredMatrices_Module
 * \class Circulant
 * \brief An \c n x \c n circulant matrix represented by its first column.
 *
 * For first column \f$c\f$ and unitary DFT matrix \f$F\f$,
 * \f[ C_{ij}=c_{(i-j)\bmod n}, \qquad C=F^*\operatorname{diag}(\sqrt{n}\,Fc)F. \f]
 * Thus its eigenvalues -- the operator's \em symbol -- are \f$\sqrt{n}\,Fc\f$,
 * the DFT of \f$c\f$, computed once at construction and reused by every
 * product. This yields an O(n log n) matrix-vector product (\c operator*), an
 * O(n log n) direct (pseudo-inverse) solve (\ref solve), and closed-form
 * factorizations: the eigendecomposition (\ref eigenvalues, \ref eigenvectors)
 * and the SVD (\ref singularValues, \ref matrixU, \ref matrixV) in the Fourier
 * basis, plus \ref rank, \ref inverse and \ref determinant. The class is closed
 * under \ref transpose, \ref conjugate and \ref adjoint, which reuse the cached
 * symbol instead of recomputing FFTs.
 *
 * The operator stores its own copy of the generating column and derives from
 * \c EigenBase. Because \c operator* returns an Eigen product expression, a
 * \c Circulant also drops into the matrix-free iterative solvers, and it can be
 * assigned to a dense matrix when an explicit representation is needed. As with
 * any matrix-free operator, the iterative solvers must be instantiated with
 * \c IdentityPreconditioner (e.g.
 * \c ConjugateGradient<Circulant<double>,Lower|Upper,IdentityPreconditioner>):
 * the default preconditioners read individual coefficients through \c col() or
 * \c InnerIterator, which the structured operators do not expose.
 *
 * \tparam Scalar_ the scalar type, real or complex.
 * \tparam Size_ the dimension at compile time, or \c Dynamic (the default).
 *
 * \par Examples
 * \ref StructuredExampleLoveEquation "Love's integral equation and a preconditioned Toeplitz system",
 * \ref StructuredExampleTorusHelmholtz "the screened Poisson equation on a torus",
 * \ref StructuredExampleSpectralDifferentiation "Fourier spectral differentiation".
 *
 * \sa class Toeplitz, makeCirculant()
 */
template <typename Scalar_, int Size_>
class Circulant : public EigenBase<Circulant<Scalar_, Size_>> {
 public:
  using Scalar = Scalar_;
  using RealScalar = typename NumTraits<Scalar>::Real;
  using StorageIndex = int;
  using Complex = std::complex<RealScalar>;
  using GeneratorType = Matrix<Scalar, Size_, 1>;
  using ComplexVector = Matrix<Complex, Dynamic, 1>;
  using RealVector = Matrix<RealScalar, Size_, 1>;
  using ComplexMatrix = Matrix<Complex, Size_, Size_>;

  static constexpr int RowsAtCompileTime = Size_;
  static constexpr int ColsAtCompileTime = Size_;
  static constexpr int MaxRowsAtCompileTime = Size_;
  static constexpr int MaxColsAtCompileTime = Size_;
  static constexpr int SizeAtCompileTime = internal::size_at_compile_time(Size_, Size_);
  static constexpr int MaxSizeAtCompileTime = SizeAtCompileTime;
  static constexpr bool IsRowMajor = false;
  // Deliberately no IsVectorAtCompileTime: Ref<const Circulant>'s default StrideType
  // argument reads it, so its absence makes internal::is_ref_compatible SFINAE to
  // false and keeps the iterative solvers on their matrix-free path.

  /** Builds a circulant matrix from its first column \a col.
   *
   * When the matrix is large enough for products to take the FFT path, the DFT of
   * \a col -- the eigenvalues of the matrix -- is computed here, once, and reused
   * by every subsequent product and solve. */
  template <typename Derived>
  explicit Circulant(const MatrixBase<Derived>& col) : m_col(col) {
    EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
    EIGEN_INTERNAL_ASSERT_NONEMPTY_MSG(m_col, "Circulant generator must be non-empty");
    if (m_col.size() > internal::structured_direct_threshold()) {
      m_symbol = computeSymbol();
      m_prodSymbol = computeProdSymbol(m_col);
    }
  }

  EIGEN_DEVICE_FUNC Index rows() const { return m_col.size(); }
  EIGEN_DEVICE_FUNC Index cols() const { return m_col.size(); }

  /** \returns the generating first column. */
  const GeneratorType& column() const { return m_col; }

  /** \returns the symbol of the operator, i.e. the DFT of the generating column.
   * Its entries are the eigenvalues of the matrix. Cached when the operator is
   * large enough for products to take the FFT path, computed on the fly for small
   * operators. */
  ComplexVector symbol() const { return m_symbol.size() > 0 ? m_symbol : computeSymbol(); }

  /** \returns the coefficient at row \a row and column \a col. */
  Scalar coeff(Index row, Index col) const {
    Index k = row - col;
    if (k < 0) k += rows();
    return m_col.coeff(k);
  }

  /** \returns the transpose of \c *this, itself a \c Circulant operator: the one
   * generated by the index-reversed column. The cached symbols, when present, are
   * reused -- the symbols of the transpose are the index reversals of the symbols
   * -- so no FFT is recomputed. */
  Circulant transpose() const {
    const Index n = rows();
    GeneratorType col(n);
    col[0] = m_col[0];
    if (n > 1) col.tail(n - 1) = m_col.tail(n - 1).reverse();
    return Circulant(col, internal::structured_reverse_symbol(m_symbol),
                     internal::structured_reverse_symbol(m_prodSymbol));
  }

  /** \returns the complex conjugate of \c *this, itself a \c Circulant operator.
   * The cached symbols, when present, are reused: the symbols of the conjugate are
   * the conjugated index reversals of the symbols. */
  Circulant conjugate() const {
    return Circulant(m_col.conjugate(), internal::structured_reverse_symbol(m_symbol).conjugate(),
                     internal::structured_reverse_symbol(m_prodSymbol).conjugate());
  }

  /** \returns the adjoint of \c *this, itself a \c Circulant operator. The cached
   * symbols, when present, are reused: the symbols of the adjoint are the
   * elementwise conjugates of the symbols (the eigenvalues conjugate while the
   * Fourier eigenbasis stays fixed). */
  Circulant adjoint() const {
    const Index n = rows();
    GeneratorType col(n);
    col[0] = numext::conj(m_col[0]);
    if (n > 1) col.tail(n - 1) = m_col.tail(n - 1).reverse().conjugate();
    return Circulant(col, m_symbol.conjugate(), m_prodSymbol.conjugate());
  }

  /** \internal Writes the dense representation into \a dst; column \c j is the
   * generator rotated downwards by \c j. Row-major destinations use reversed
   * generator segments along each row. Invoked through \c dense = circulant; */
  template <typename Dest>
  void evalTo(Dest& dst) const {
    const Index n = rows();
    EIGEN_IF_CONSTEXPR (Dest::IsRowMajor) {
      for (Index i = 0; i < n; ++i) {
        dst.row(i).head(i + 1) = m_col.head(i + 1).reverse().transpose();
        dst.row(i).tail(n - i - 1) = m_col.tail(n - i - 1).reverse().transpose();
      }
      return;
    }
    for (Index j = 0; j < n; ++j) {
      dst.col(j).head(j) = m_col.tail(j);
      dst.col(j).tail(n - j) = m_col.head(n - j);
    }
  }

  /** \internal Computes \c dst += (*this), see evalTo(). */
  template <typename Dest>
  void addTo(Dest& dst) const {
    const Index n = rows();
    EIGEN_IF_CONSTEXPR (Dest::IsRowMajor) {
      for (Index i = 0; i < n; ++i) {
        dst.row(i).head(i + 1) += m_col.head(i + 1).reverse().transpose();
        dst.row(i).tail(n - i - 1) += m_col.tail(n - i - 1).reverse().transpose();
      }
      return;
    }
    for (Index j = 0; j < n; ++j) {
      dst.col(j).head(j) += m_col.tail(j);
      dst.col(j).tail(n - j) += m_col.head(n - j);
    }
  }

  /** \internal Computes \c dst -= (*this), see evalTo(). */
  template <typename Dest>
  void subTo(Dest& dst) const {
    const Index n = rows();
    EIGEN_IF_CONSTEXPR (Dest::IsRowMajor) {
      for (Index i = 0; i < n; ++i) {
        dst.row(i).head(i + 1) -= m_col.head(i + 1).reverse().transpose();
        dst.row(i).tail(n - i - 1) -= m_col.tail(n - i - 1).reverse().transpose();
      }
      return;
    }
    for (Index j = 0; j < n; ++j) {
      dst.col(j).head(j) -= m_col.tail(j);
      dst.col(j).tail(n - j) -= m_col.head(n - j);
    }
  }

  /** \returns the product expression \c (*this) * \a x, evaluated through a fast
   * FFT-based matrix-vector product. The expression carries the default product
   * tag, so assigning it behaves like any dense product: a temporary resolves
   * aliasing between the destination and \a x, and \c .noalias() skips it. */
  template <typename Rhs>
  Product<Circulant, Rhs> operator*(const MatrixBase<Rhs>& x) const {
    return Product<Circulant, Rhs>(*this, x.derived());
  }

  /** \returns the minimum-norm least-squares solution of \c (*this) * x = b,
   * computed directly in the Fourier domain. Symbol entries whose modulus reaches
   * the rank threshold (see \ref rank) are inverted; the remaining ones are
   * treated as exact zeros, so the result is the pseudo-inverse applied to \a b.
   * For a non-singular operator this is the exact solution
   * \c x = ifft( fft(b) ./ fft(c) ). Supports multiple right-hand sides. */
  template <typename Rhs>
  Matrix<Scalar, Size_, Rhs::ColsAtCompileTime> solve(const MatrixBase<Rhs>& b) const {
    EIGEN_INTERNAL_ASSERT_SOLVER_RHS_MSG(*this, b, false, "right-hand side has the wrong number of rows");
    const Index n = rows();
    const ComplexVector s = symbol();
    const RealVector mods = s.cwiseAbs();
    const ComplexVector sinv = internal::structured_pinv_symbol(s, mods, internal::structured_rank_threshold(s, mods));
    Matrix<Scalar, Size_, Rhs::ColsAtCompileTime> x(n, b.cols());
    x.setZero();
    internal::structured_fft_apply(x, sinv, n, b.derived(), Scalar(1));
    return x;
  }

  /** \returns the numerical rank: the number of symbol entries whose modulus is
   * no smaller than the threshold \c n * epsilon * max_k|symbol[k]|, clamped from
   * below by the smallest normal number. This is the same threshold \ref solve
   * uses to decide which Fourier components to invert, and the comparison is
   * strict like SVDBase's, so an entry sitting exactly on the threshold still
   * counts as non-zero. */
  Index rank() const {
    const ComplexVector s = symbol();
    const RealVector mods = s.cwiseAbs();
    return (!(mods.array() < internal::structured_rank_threshold(s, mods))).count();  // NaN entries count as non-zero
  }

  /** \returns the inverse of \c *this, itself a \c Circulant operator: the one
   * generated by \c ifft(1 ./ symbol), the first column of the inverse matrix.
   * \warning The operator must be non-singular; use \ref solve for a
   * pseudo-inverse solve of a rank-deficient operator. */
  Circulant inverse() const {
    const Index n = rows();
    const ComplexVector sinv = symbol().cwiseInverse();
    GeneratorType col(n);
    if (n == 1) {
      col = internal::structured_scalar_part_impl<Scalar>::run(sinv);
    } else {
      auto&& fft = internal::structured_fft_engine<RealScalar>();
      ComplexVector ct(n);
      fft.inv(ct, sinv, n);
      col = internal::structured_scalar_part_impl<Scalar>::run(ct);
    }
    return Circulant(col, n > internal::structured_direct_threshold() ? sinv : ComplexVector(),
                     n > internal::structured_direct_threshold() ? computeProdSymbol(col) : ComplexVector());
  }

  /** \returns the determinant, i.e. the product of the eigenvalues (the symbol
   * entries). The product is accumulated in the balanced form \c m * 2^e (the
   * split fraction/exponent determinant convention of LINPACK's xGEDI [4]) --
   * every factor and the running product are renormalized to unit magnitude with
   * the power of two tracked separately -- so the partial products can neither
   * overflow nor underflow when the determinant itself is representable, whatever
   * the ordering of large and small eigenvalues. For a real operator the product
   * is real up to roundoff, and its real part is returned. */
  Scalar determinant() const {
    const ComplexVector s = symbol();
    Complex det(1);
    Index exponent = 0;
    for (Index k = 0; k < s.size(); ++k)
      det = internal::structured_balance(Complex(det * internal::structured_balance(s[k], exponent)), exponent);
    return internal::structured_scalar_part_impl<Scalar>::run_scalar(internal::structured_ldexp_clamped(det, exponent));
  }

  /** \returns the eigenvalues: eigenvalue \c k is \c symbol()[k], and its
   * (unit-norm) eigenvector is the Fourier vector \c f_k with
   * \f$ (f_k)_j = e^{2\pi i j k / n} / \sqrt{n} \f$, see \ref eigenvectors.
   * Every circulant matrix is diagonalized by this same Fourier basis. */
  ComplexVector eigenvalues() const { return symbol(); }

  /** \returns the unitary matrix of eigenvectors: column \c k is the Fourier
   * vector \c f_k matching \c eigenvalues()[k].
   * \note The eigenvector matrix is materialized as a dense \c n x \c n matrix;
   * unlike the other methods of this class this costs O(n^2) storage. */
  ComplexMatrix eigenvectors() const {
    const Index n = rows();
    const ComplexVector roots = fourierRoots();
    ComplexMatrix F(n, n);
    for (Index k = 0; k < n; ++k) fourierColumn(F, roots, k, k);
    return F;
  }

  /** \returns the singular values, sorted in decreasing order: the moduli of the
   * symbol entries. The ordering is shared with \ref matrixU and \ref matrixV, so
   * together they form the SVD \c *this = U * singularValues().asDiagonal() * V^H.
   * A real operator has a conjugate-symmetric symbol, so the two modes of a
   * conjugate pair carry exactly equal singular values. */
  RealVector singularValues() const {
    ComplexVector s;
    RealVector mods;
    const std::vector<Index> perm = svdOrdering(s, mods);
    RealVector sv(s.size());
    for (Index t = 0; t < s.size(); ++t) sv[t] = mods[perm[t]];
    return sv;
  }

  /** \returns the matrix of left singular vectors \c U: column \c t is the Fourier
   * vector of the t-th largest symbol entry, scaled by its phase (phase 1 for a
   * zero entry). Dense \c n x \c n, see the note in \ref eigenvectors. */
  ComplexMatrix matrixU() const {
    const Index n = rows();
    ComplexVector s;
    RealVector mods;
    const std::vector<Index> perm = svdOrdering(s, mods);
    const ComplexVector roots = fourierRoots();
    ComplexMatrix U(n, n);
    for (Index t = 0; t < n; ++t) {
      fourierColumn(U, roots, perm[t], t);
      const RealScalar a = mods[perm[t]];
      if (a > RealScalar(0)) U.col(t) *= s[perm[t]] / a;
    }
    return U;
  }

  /** \returns the matrix of right singular vectors \c V: column \c t is the
   * Fourier vector of the t-th largest symbol entry. Dense \c n x \c n, see the
   * note in \ref eigenvectors. */
  ComplexMatrix matrixV() const {
    const Index n = rows();
    ComplexVector s;
    RealVector mods;
    const std::vector<Index> perm = svdOrdering(s, mods);
    const ComplexVector roots = fourierRoots();
    ComplexMatrix V(n, n);
    for (Index t = 0; t < n; ++t) fourierColumn(V, roots, perm[t], t);
    return V;
  }

  /** \internal Computes \c dst += alpha * (*this) * rhs. \c ProductScalar is the
   * promoted scalar of the product (complex when a real operator is applied to a
   * complex right-hand side); the accumulation runs in the promoted type. */
  template <typename Dest, typename Rhs, typename ProductScalar>
  void addProduct(Dest& dst, const Rhs& rhs, const ProductScalar& alpha) const {
    const Index n = rows();
    EIGEN_INTERNAL_ASSERT_PRODUCT_SHAPE_MSG(*this, rhs, "invalid product: dimensions do not match");
    if (n <= internal::structured_direct_threshold()) {
      directProduct(dst, rhs, alpha);
      return;
    }
    // Products use the padded embedding symbol when the operator size is not
    // 5-smooth (see computeProdSymbol()); the exact-size transform of such a
    // size runs through kissfft's quadratic generic butterfly.
    const ComplexVector& s = m_prodSymbol.size() > 0 ? m_prodSymbol : m_symbol;
    internal::structured_fft_apply(dst, s, n, rhs, alpha);
  }

 private:
  /** \internal Builds an operator from a generator and already-known symbols
   * (empty for small operators; \a prodSymbol also empty for 5-smooth sizes,
   * which need no padded product embedding), skipping the FFTs of the public
   * constructor. Used by transpose(), conjugate(), adjoint() and inverse(),
   * whose symbols are cheap transformations of the existing ones. */
  Circulant(const GeneratorType& col, const ComplexVector& symbol, const ComplexVector& prodSymbol)
      : m_col(col), m_symbol(symbol), m_prodSymbol(prodSymbol) {}

  /** \internal \returns the symbol products should use for a generator \a col of
   * size n: empty when n is 5-smooth (the exact-size transform is already fast),
   * otherwise the DFT of the circulant embedding of size
   * \c p = fft_next_good_size(2n - 1) -- the length-n cyclic convolution
   * evaluated as a padded linear convolution, exactly the Toeplitz embedding of
   * the circulant. kissfft falls back to a quadratic generic butterfly for prime
   * factors other than 2, 3 and 5, so the exact-size transform of e.g. a prime n
   * costs orders of magnitude more than the padded one. The spectral operations
   * (eigenvalues, solve, rank, determinant, SVD) keep the exact-size symbol,
   * whose entries are the eigenvalues. */
  static ComplexVector computeProdSymbol(const GeneratorType& col) {
    const Index n = col.size();
    if (internal::fft_next_good_size(n) == n) return ComplexVector();
    const Index p = internal::fft_next_good_size(2 * n - 1);
    ComplexVector embedding = ComplexVector::Zero(p);
    embedding.head(n) = col.template cast<Complex>();
    embedding.tail(n - 1) = col.tail(n - 1).template cast<Complex>();
    ComplexVector symbol(p);
    auto&& fft = internal::structured_fft_engine<RealScalar>();
    fft.fwd(symbol, embedding, p);
    return symbol;
  }

  /** \internal Direct O(n^2) kernel for column \a k of the right-hand side:
   * computes \c dst.col(k) += alpha * (*this) * rhs.col(k) without transforms,
   * for operators below the FFT threshold. */
  template <typename Dest, typename Rhs, typename ProductScalar>
  void directProductColumn(Dest& dst, const Rhs& rhs, Index k, const ProductScalar& alpha) const {
    const Index n = rows();
    // A unit alpha must not multiply: even the identity complex scalar (1,0)
    // pollutes an (Inf,0) value with NaN through the 0*Inf cross term.
    const bool unitAlpha = alpha == ProductScalar(1);
    if (n <= internal::structured_scalar_threshold()) {
      // Tiny sizes: a plain scalar loop beats the segment-based path below, whose
      // per-segment setup dominates when segments hold only a few entries.
      for (Index i = 0; i < n; ++i) {
        ProductScalar acc(0);
        for (Index j = 0; j < n; ++j) acc += coeff(i, j) * rhs.coeff(j, k);
        dst.coeffRef(i, k) += unitAlpha ? acc : ProductScalar(alpha * acc);
      }
      return;
    }

    // Column j is c rotated downward by j; split it into two contiguous segments.
    auto dstCol = dst.col(k);
    for (Index j = 0; j < n; ++j) {
      const ProductScalar xj = unitAlpha ? ProductScalar(rhs.coeff(j, k)) : ProductScalar(alpha * rhs.coeff(j, k));
      dstCol.head(j) += xj * m_col.tail(j);
      dstCol.tail(n - j) += xj * m_col.head(n - j);
    }
  }

  /** \internal Direct O(n^2) product kernel over every column, see
   * directProductColumn(). */
  template <typename Dest, typename Rhs, typename ProductScalar>
  void directProduct(Dest& dst, const Rhs& rhs, const ProductScalar& alpha) const {
    for (Index k = 0; k < rhs.cols(); ++k) directProductColumn(dst, rhs, k, alpha);
  }

  /** \internal Computes the symbol \a s and the moduli \a mods that order the
   * SVD, and \returns the permutation sorting those moduli in decreasing order --
   * the single ordering singularValues(), matrixU() and matrixV() agree on.
   *
   * For a real generator the symbol is conjugate-symmetric, so modes \c k and
   * \c n-k have equal moduli; taken through the complex transform they agree only
   * to roundoff. Ranking such a pair by the computed moduli leaves its order
   * decided by those last bits, and the three accessors each derive the ordering
   * from their own symbol computation, which nothing requires to agree down to
   * the last bit: one of them ranking a pair the other way pairs a left singular
   * vector with the wrong right one and breaks the reconstruction by an O(1)
   * amount. The moduli of a conjugate pair are therefore tied exactly, which the
   * stable sort resolves by index in every caller. Sharing one value between the
   * pair costs nothing: it cancels out of \c U * S, where matrixU() divides the
   * mode's phase by it and singularValues() returns it.
   */
  std::vector<Index> svdOrdering(ComplexVector& s, RealVector& mods) const {
    const Index n = rows();
    s = symbol();
    EIGEN_IF_CONSTEXPR (!NumTraits<Scalar>::IsComplex) {
      const Index pairs = (n - 1) / 2;
      mods.resize(n);
      mods.head(n - pairs) = s.head(n - pairs).cwiseAbs();
      mods.tail(pairs) = mods.segment(1, pairs).reverse();
    } else {
      mods = s.cwiseAbs();
    }
    return internal::structured_svd_permutation(mods);
  }

  // Cache the n distinct roots once per Fourier matrix, preserving polar()'s
  // bounded arguments and rounding while avoiding n^2 trigonometric evaluations.
  ComplexVector fourierRoots() const {
    const Index n = rows();
    const RealScalar scale = RealScalar(1) / numext::sqrt(RealScalar(n));
    ComplexVector roots(n);
    for (Index j = 0; j < n; ++j)
      roots[j] = std::polar(scale, RealScalar(2 * EIGEN_PI) * RealScalar(j) / RealScalar(n));
    return roots;
  }

  // Column k uses roots[j*k mod n]; incremental modular indexing avoids j*k overflow.
  void fourierColumn(ComplexMatrix& F, const ComplexVector& roots, Index k, Index dstCol) const {
    const Index n = rows();
    Index jk = 0;
    for (Index j = 0; j < n; ++j) {
      F(j, dstCol) = roots[jk];
      jk += k;
      if (jk >= n) jk -= n;
    }
  }

  /** \internal \returns the DFT of the generating column. */
  ComplexVector computeSymbol() const {
    const Index n = m_col.size();
    const ComplexVector cc = m_col.template cast<Complex>();
    if (n == 1) return cc;  // the DFT of a single sample is the identity
    ComplexVector symbol(n);
    auto&& fft = internal::structured_fft_engine<RealScalar>();
    fft.fwd(symbol, cc, n);
    return symbol;
  }

  GeneratorType m_col;
  ComplexVector m_symbol;
  // The padded embedding symbol products use when the size is not 5-smooth;
  // empty otherwise. See computeProdSymbol().
  ComplexVector m_prodSymbol;
};

/** \ingroup StructuredMatrices_Module
 * \returns a \ref Circulant operator with first column \a col. The compile-time
 * size of the operator is deduced from \a col. */
template <typename Derived>
Circulant<typename Derived::Scalar, Derived::SizeAtCompileTime> makeCirculant(const MatrixBase<Derived>& col) {
  return Circulant<typename Derived::Scalar, Derived::SizeAtCompileTime>(col);
}

namespace internal {

// StructuredShape avoids ambiguity with the stock dense product specializations.
template <typename Scalar_, int Size_, typename Rhs, int ProductTag>
struct generic_product_impl<Circulant<Scalar_, Size_>, Rhs, StructuredShape, DenseShape, ProductTag>
    : structured_product_impl<Circulant<Scalar_, Size_>, Rhs> {};

}  // namespace internal

#if EIGEN_COMP_CXXVER < 17
template <typename Scalar_, int Size_>
constexpr int Circulant<Scalar_, Size_>::RowsAtCompileTime;
template <typename Scalar_, int Size_>
constexpr int Circulant<Scalar_, Size_>::ColsAtCompileTime;
template <typename Scalar_, int Size_>
constexpr int Circulant<Scalar_, Size_>::MaxRowsAtCompileTime;
template <typename Scalar_, int Size_>
constexpr int Circulant<Scalar_, Size_>::MaxColsAtCompileTime;
template <typename Scalar_, int Size_>
constexpr int Circulant<Scalar_, Size_>::SizeAtCompileTime;
template <typename Scalar_, int Size_>
constexpr int Circulant<Scalar_, Size_>::MaxSizeAtCompileTime;
template <typename Scalar_, int Size_>
constexpr bool Circulant<Scalar_, Size_>::IsRowMajor;
#endif

}  // namespace Eigen

#endif  // EIGEN_STRUCTURED_CIRCULANT_H