diff --git a/failtest/CMakeLists.txt b/failtest/CMakeLists.txt
index a9f9467..4ddae13 100644
--- a/failtest/CMakeLists.txt
+++ b/failtest/CMakeLists.txt
@@ -71,6 +71,9 @@
 
 ei_add_failtest("erf_no_scalar_overload")
 ei_add_failtest("erfc_no_scalar_overload")
+ei_add_failtest("vandermonde_int")
+ei_add_failtest("vandermonde_rectangular_square")
+ei_add_failtest("bjorckpereyra_rectangular")
 
 ei_add_failtest("structured_bindings_dynamic_matrix")
 ei_add_failtest("structured_bindings_dynamic_array")
diff --git a/failtest/bjorckpereyra_rectangular.cpp b/failtest/bjorckpereyra_rectangular.cpp
new file mode 100644
index 0000000..f49ca3b
--- /dev/null
+++ b/failtest/bjorckpereyra_rectangular.cpp
@@ -0,0 +1,18 @@
+// SPDX-FileCopyrightText: The Eigen Authors
+// SPDX-License-Identifier: MPL-2.0
+
+#include "../unsupported/Eigen/StructuredMatrices"
+
+#ifdef EIGEN_SHOULD_FAIL_TO_BUILD
+constexpr int Cols = 3;
+#else
+constexpr int Cols = 2;
+#endif
+
+int main() {
+  Eigen::Vector2d nodes;
+  nodes << 0.0, 1.0;
+  Eigen::Vandermonde<double, 2, Cols> vandermonde(nodes, Cols);
+  Eigen::BjorckPereyra<double> solver(vandermonde);
+  return static_cast<int>(solver.rows());
+}
diff --git a/failtest/vandermonde_int.cpp b/failtest/vandermonde_int.cpp
new file mode 100644
index 0000000..9f8b2af
--- /dev/null
+++ b/failtest/vandermonde_int.cpp
@@ -0,0 +1,17 @@
+// SPDX-FileCopyrightText: The Eigen Authors
+// SPDX-License-Identifier: MPL-2.0
+
+#include "../unsupported/Eigen/StructuredMatrices"
+
+#ifdef EIGEN_SHOULD_FAIL_TO_BUILD
+using Scalar = int;
+#else
+using Scalar = double;
+#endif
+
+int main() {
+  Eigen::Matrix<Scalar, 2, 1> nodes;
+  nodes << Scalar(1), Scalar(2);
+  Eigen::Vandermonde<Scalar, 2, 2> vandermonde(nodes);
+  return static_cast<int>(vandermonde.rows());
+}
diff --git a/failtest/vandermonde_rectangular_square.cpp b/failtest/vandermonde_rectangular_square.cpp
new file mode 100644
index 0000000..6a682dc
--- /dev/null
+++ b/failtest/vandermonde_rectangular_square.cpp
@@ -0,0 +1,17 @@
+// SPDX-FileCopyrightText: The Eigen Authors
+// SPDX-License-Identifier: MPL-2.0
+
+#include "../unsupported/Eigen/StructuredMatrices"
+
+#ifdef EIGEN_SHOULD_FAIL_TO_BUILD
+constexpr int Cols = 3;
+#else
+constexpr int Cols = 2;
+#endif
+
+int main() {
+  Eigen::Vector2d nodes;
+  nodes << 0.0, 1.0;
+  Eigen::Vandermonde<double, 2, Cols> vandermonde(nodes);
+  return static_cast<int>(vandermonde.cols());
+}
diff --git a/unsupported/Eigen/StructuredMatrices b/unsupported/Eigen/StructuredMatrices
index 4e9d160..310beca 100644
--- a/unsupported/Eigen/StructuredMatrices
+++ b/unsupported/Eigen/StructuredMatrices
@@ -36,11 +36,17 @@
  *     operands, and diagonal (in particular identity) factors are stored and
  *     applied in diagonal form;
  *   - \c DiagonalPlusLowRank : a diagonal matrix plus a rank-k correction, with
- *     O(nk) products and O(nk^2) Woodbury solves, closed under inversion.
+ *     O(nk) products and O(nk^2) Woodbury solves, closed under inversion;
+ *   - \c Vandermonde : a Vandermonde matrix stored as its nodes, with Horner
+ *     products and O(n^2) Björck-Pereyra primal/dual solves
+ *     (\c BjorckPereyra).
  *
  * The operator types derive from \c EigenBase and store only compact generators
  * or factors. The FFT-backed operators (\c Circulant, \c Toeplitz and \c Hankel)
  * also keep a precomputed DFT symbol that every product reuses.
+ * The operator types derive from \c EigenBase and store only their generating
+ * vectors; the FFT-backed operators (\c Circulant, \c Toeplitz, \c Bccb) also
+ * keep a precomputed DFT symbol that every product reuses.
  * Because they expose \c operator* returning an Eigen product expression, they
  * also plug directly into the matrix-free iterative solvers
  * (\c ConjugateGradient, \c GMRES, \c MINRES, ...) without forming the dense
@@ -55,6 +61,11 @@
  * ones reusing the cached symbol), which in particular feeds the least-squares
  * solvers \c LSMR and \c LeastSquaresConjugateGradient (again with
  * \c IdentityPreconditioner).
+ * The \c Circulant, \c Toeplitz and \c Hankel operators are closed under
+ * transposition (\c transpose(),
+ * \c conjugate(), \c adjoint() return operators of the same kind, reusing the
+ * cached symbol), which in particular feeds the least-squares solvers \c LSMR
+ * and \c LeastSquaresConjugateGradient (again with \c IdentityPreconditioner).
  * \c Circulant additionally exposes its closed-form eigendecomposition and SVD
  * in the Fourier basis, a pseudo-inverse (minimum-norm least-squares) solve,
  * \c rank(), \c inverse() and \c determinant().
@@ -76,6 +87,7 @@
 #include "src/StructuredMatrices/Hankel.h"
 #include "src/StructuredMatrices/KroneckerOperator.h"
 #include "src/StructuredMatrices/DiagonalPlusLowRank.h"
+#include "src/StructuredMatrices/Vandermonde.h"
 // IWYU pragma: end_exports
 
 #include "../../Eigen/src/Core/util/ReenableStupidWarnings.h"
diff --git a/unsupported/Eigen/src/StructuredMatrices/StructuredMatrixUtils.h b/unsupported/Eigen/src/StructuredMatrices/StructuredMatrixUtils.h
index 66bafba..510ba6a 100644
--- a/unsupported/Eigen/src/StructuredMatrices/StructuredMatrixUtils.h
+++ b/unsupported/Eigen/src/StructuredMatrices/StructuredMatrixUtils.h
@@ -39,6 +39,11 @@
 // fewer than a couple of packets (measured crossover on AVX2 hardware).
 constexpr Index structured_scalar_threshold() { return 16; }
 
+// Numerical scale exponents are independent of Eigen's configurable dimension
+// index. A 32-bit Index can overflow while accumulating O(n^2) factor
+// exponents for dimensions that are otherwise practical.
+using structured_exponent_type = numext::int64_t;
+
 /** \internal Balanced mantissa*2^e arithmetic shared by the structured
  * operators' determinant-style accumulations (the split fraction/exponent
  * convention of LINPACK's xGEDI; see the per-class references).
@@ -51,7 +56,8 @@
 template <typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex>
 struct structured_balance_impl {
   using RealScalar = typename NumTraits<Scalar>::Real;
-  static Scalar run(const Scalar& z, Index& exponent) {
+  template <typename Exponent>
+  static Scalar run(const Scalar& z, Exponent& exponent) {
     const RealScalar mag = numext::maxi(numext::abs(numext::real(z)), numext::abs(numext::imag(z)));
     if (!(mag > RealScalar(0)) || !(numext::isfinite)(mag)) return z;
     int e;
@@ -67,7 +73,8 @@
 
 template <typename Scalar>
 struct structured_balance_impl<Scalar, false> {
-  static Scalar run(const Scalar& x, Index& exponent) {
+  template <typename Exponent>
+  static Scalar run(const Scalar& x, Exponent& exponent) {
     const Scalar mag = numext::abs(x);
     if (!(mag > Scalar(0)) || !(numext::isfinite)(mag)) return x;
     int e;
@@ -79,8 +86,8 @@
   static Scalar apply_exponent(const Scalar& x, int e) { return numext::ldexp(x, e); }
 };
 
-template <typename Scalar>
-Scalar structured_balance(const Scalar& z, Index& exponent) {
+template <typename Scalar, typename Exponent>
+Scalar structured_balance(const Scalar& z, Exponent& exponent) {
   return structured_balance_impl<Scalar>::run(z, exponent);
 }
 
@@ -88,9 +95,9 @@
  * component-wise for complex scalars. ldexp saturates cleanly to zero /
  * infinity (preserving signs) once the exponent leaves the representable
  * range; the clamp only guards the narrowing to int. */
-template <typename Scalar>
-Scalar structured_ldexp_clamped(const Scalar& z, Index exponent) {
-  constexpr Index kMaxExponent = Index(1) << 24;
+template <typename Scalar, typename Exponent>
+Scalar structured_ldexp_clamped(const Scalar& z, Exponent exponent) {
+  constexpr Exponent kMaxExponent = Exponent(1) << 24;
   const int e = static_cast<int>(numext::mini(numext::maxi(exponent, -kMaxExponent), kMaxExponent));
   return structured_balance_impl<Scalar>::apply_exponent(z, e);
 }
@@ -115,6 +122,27 @@
   return structured_ldexp_clamped(Scalar(1) / zs, -e);
 }
 
+/** \internal \returns \c a - b guarded against spurious overflow: when the
+ * plain difference of two finite values overflows to infinity, it is
+ * recomputed from the halved operands with \a e set to 1 so the caller can
+ * carry the factor of two in its running exponent (the balanced accumulations
+ * fold it into \c exponent at the call site). The recomputation is exact where
+ * it matters: a difference of finite values only overflows when both operands
+ * are huge, normal values, whose halves and halved difference are exactly
+ * representable. For complex scalars the finiteness tests are component-wise;
+ * a subnormal component riding along with a huge one loses its last bit to the
+ * halving, an error far below one ulp of the factor's magnitude. Non-finite
+ * operands (genuine Inf/NaN nodes) propagate untouched, with \a e = 0. */
+template <typename Scalar>
+Scalar structured_guarded_diff(const Scalar& a, const Scalar& b, int& e) {
+  using RealScalar = typename NumTraits<Scalar>::Real;
+  e = 0;
+  const Scalar t = a - b;
+  if ((numext::isfinite)(t) || !(numext::isfinite)(a) || !(numext::isfinite)(b)) return t;
+  e = 1;
+  return a * RealScalar(0.5) - b * RealScalar(0.5);
+}
+
 /** \internal \returns the indices sorted by decreasing precomputed modulus
  * \a mods (each modulus is computed once, not on every comparison); the shared
  * ordering of the operators' singularValues()/matrixU()/matrixV(). The sort is
@@ -337,7 +365,7 @@
     // The length-one DFT is the identity and is unsupported by kissfft. The
     // pointwise step still runs, scaled around just as the transforms are on the
     // general path.
-    const int budget = std::numeric_limits<RealScalar>::max_exponent - 2;
+    const int budget = NumTraits<RealScalar>::max_exponent() - 2;
     ComplexVector xf(1);
     for (Index k = 0; k < rhs.cols(); ++k) {
       int colExp;
@@ -357,9 +385,10 @@
     return;
   }
 
-  int log2p = 0;
-  for (Index t = p; t > 0; t /= 2) ++log2p;
-  const int budget = std::numeric_limits<RealScalar>::max_exponent - 2 * log2p - 2;
+  // The bit width of p (>= 1 here), an upper bound for the ceil(log2 p) of the
+  // magnitude bound above -- one bit looser at power-of-two transform lengths.
+  const int log2p = log2_floor(static_cast<std::make_unsigned_t<Index>>(p)) + 1;
+  const int budget = NumTraits<RealScalar>::max_exponent() - 2 * log2p - 2;
 
   auto&& fft = structured_fft_engine<RealScalar>();
   ComplexVector xt = ComplexVector::Zero(p);
@@ -450,12 +479,14 @@
   template <typename Dest>
   static void evalTo(Dest& dst, const Op& lhs, const Rhs& rhs) {
     dst.setZero();
-    lhs.addProduct(dst, rhs, Scalar(1));
+    scaleAndAddTo(dst, lhs, rhs, Scalar(1));
   }
 
   template <typename Dest>
   static void scaleAndAddTo(Dest& dst, const Op& lhs, const Rhs& rhs, const Scalar& alpha) {
-    lhs.addProduct(dst, rhs, alpha);
+    using RhsNested = typename nested_eval<Rhs, Op::RowsAtCompileTime>::type;
+    RhsNested actualRhs(rhs);
+    lhs.addProduct(dst, actualRhs, alpha);
   }
 };
 
diff --git a/unsupported/Eigen/src/StructuredMatrices/Vandermonde.h b/unsupported/Eigen/src/StructuredMatrices/Vandermonde.h
new file mode 100644
index 0000000..f218f09
--- /dev/null
+++ b/unsupported/Eigen/src/StructuredMatrices/Vandermonde.h
@@ -0,0 +1,770 @@
+// 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] N. J. Higham, "Accuracy and Stability of Numerical Algorithms", 2nd ed.,
+//      SIAM, 2002, chapter 27. Avoiding spurious overflow by rescaling with
+//      powers of two, the technique behind determinant()'s balanced
+//      accumulation and the scaled Horner recurrence of addProduct().
+//  [2] P. H. Sterbenz, "Floating-Point Computation", Prentice-Hall, 1974.
+//      Scaling by a power of two is exact, the property the balanced
+//      accumulation and the scaled Horner recurrence rely on.
+//  [3] 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.
+//  [4] L. Reichel, "Newton Interpolation at Leja Points", BIT 30 (1990),
+//      332--346. Leja ordering controls growth in the Newton representation;
+//      BjorckPereyra uses it for genuinely complex node sets.
+
+#ifndef EIGEN_STRUCTURED_VANDERMONDE_H
+#define EIGEN_STRUCTURED_VANDERMONDE_H
+
+// IWYU pragma: private
+#include "./InternalHeaderCheck.h"
+
+namespace Eigen {
+
+template <typename Scalar_, int Rows_ = Dynamic, int Cols_ = Dynamic>
+class Vandermonde;
+
+template <typename Scalar_>
+class BjorckPereyra;
+
+namespace internal {
+
+template <typename Scalar_, int Rows_, int Cols_>
+struct traits<Vandermonde<Scalar_, Rows_, Cols_>> {
+  using Scalar = Scalar_;
+  using StorageKind = Dense;
+  using XprKind = MatrixXpr;
+  using StorageIndex = int;
+  static constexpr int RowsAtCompileTime = Rows_;
+  static constexpr int ColsAtCompileTime = Cols_;
+  static constexpr int MaxRowsAtCompileTime = Rows_;
+  static constexpr int MaxColsAtCompileTime = Cols_;
+  // Deliberately no NestByRefBit: the makeVandermonde() factories (and any
+  // function returning the operator by value) produce 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(m), negligible
+  // against the O(mn) product evaluation.
+  static constexpr int Flags = Rows_ == 1 && Cols_ != 1 ? RowMajorBit : 0;
+};
+
+template <typename Scalar_, int Rows_, int Cols_>
+struct evaluator_traits<Vandermonde<Scalar_, Rows_, Cols_>> {
+  using Kind = IndexBased;
+  using Shape = StructuredShape;
+};
+
+// Core rewrites alpha * (lhs * rhs) as (alpha * lhs) * rhs. The scaled
+// Vandermonde wrapper needs coefficient and BLAS metadata to participate in
+// that general dense-expression machinery.
+template <typename Scalar_, int Rows_, int Cols_>
+struct evaluator<Vandermonde<Scalar_, Rows_, Cols_>> : evaluator_base<Vandermonde<Scalar_, Rows_, Cols_>> {
+  using XprType = Vandermonde<Scalar_, Rows_, Cols_>;
+  using Scalar = Scalar_;
+  static constexpr int CoeffReadCost = HugeCost;
+  static constexpr int Flags = traits<XprType>::Flags;
+  static constexpr int Alignment = 0;
+
+  EIGEN_DEVICE_FUNC constexpr EIGEN_STRONG_INLINE explicit evaluator(const XprType& xpr) : m_xpr(xpr) {}
+
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar coeff(Index row, Index col) const { return m_xpr.coeff(row, col); }
+
+ private:
+  const XprType& m_xpr;
+};
+
+template <typename Scalar_, int Rows_, int Cols_>
+struct blas_traits<Vandermonde<Scalar_, Rows_, Cols_>> {
+  using XprType = Vandermonde<Scalar_, Rows_, Cols_>;
+  using Scalar = Scalar_;
+  using ExtractType = const XprType&;
+  using ExtractType_ = XprType;
+  using DirectLinearAccessType = XprType;
+  static constexpr bool IsComplex = NumTraits<Scalar>::IsComplex;
+  static constexpr bool IsTransposed = false;
+  static constexpr bool NeedToConjugate = false;
+  static constexpr bool HasUsableDirectAccess = false;
+  static constexpr bool HasScalarFactor = false;
+  static EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE ExtractType extract(const XprType& x) { return x; }
+  static EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE Scalar extractScalarFactor(const XprType&) { return Scalar(1); }
+};
+
+template <typename Scalar_>
+struct traits<BjorckPereyra<Scalar_>> : traits<Matrix<Scalar_, Dynamic, Dynamic>> {
+  using XprKind = MatrixXpr;
+  using StorageKind = SolverStorage;
+  using StorageIndex = int;
+  using BaseTraits = traits<Matrix<Scalar_, Dynamic, Dynamic>>;
+  static constexpr int Flags = BaseTraits::Flags & RowMajorBit;
+  static constexpr int CoeffReadCost = Dynamic;
+};
+
+}  // namespace internal
+
+/** \ingroup StructuredMatrices_Module
+ * \class Vandermonde
+ * \brief An \c m x \c n Vandermonde matrix represented by its node vector.
+ *
+ * A Vandermonde matrix has entry \c (i,j) equal to \f$ x_i^j \f$, where \c x is
+ * the node vector. Thus
+ * \f[ (Va)_i = \sum_{j=0}^{n-1} a_j x_i^j, \qquad
+ *     p_{n-1}=a_{n-1},\quad p_j=a_j+x_i p_{j+1}. \f]
+ * The class stores only the \c m nodes; products evaluate this Horner recurrence
+ * rule at O(mn) operations -- the same cost as a dense product, but with O(m)
+ * storage and without ever forming the matrix.
+ *
+ * Square systems are solved in O(n^2) by the Björck-Pereyra algorithm (class
+ * \ref BjorckPereyra), whose \c transpose().solve() form covers the dual
+ * (moment) system. There is no fast transposed \em product, so the class is not
+ * closed under transposition and rectangular least-squares problems are best
+ * handled by a dense QR of the materialized matrix.
+ *
+ * Because \c operator* returns an Eigen product expression, a \c Vandermonde
+ * 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 BiCGSTAB<Vandermonde<double>,IdentityPreconditioner>): the default
+ * preconditioners read individual coefficients through \c col() or
+ * \c InnerIterator, which the structured operators do not expose.
+ *
+ * \warning Vandermonde matrices with real nodes are exponentially
+ * ill-conditioned: the condition number grows at least like \f$ 2^n \f$ for any
+ * real node configuration (Beckermann, 2000). Solves remain surprisingly
+ * accurate for monotone node sets and sign-alternating right-hand sides
+ * (Björck-Pereyra's celebrated property, see Higham, ASNA ch. 22), but forward
+ * errors necessarily scale with the conditioning in general. Complex nodes on
+ * the unit circle are the well-conditioned case: for the n-th roots of unity,
+ * \f$ V/\sqrt{n} \f$ is unitary.
+ *
+ * \tparam Scalar_ a floating-point-like real or complex scalar supporting
+ * Eigen's scalar math hooks, including \c isfinite, \c frexp and \c ldexp (and
+ * \c log for complex solver ordering). Integer types are rejected.
+ * \tparam Rows_ the number of rows (nodes) at compile time, or \c Dynamic.
+ * \tparam Cols_ the number of columns (powers) at compile time, or \c Dynamic.
+ *
+ * \sa class BjorckPereyra, makeVandermonde()
+ */
+template <typename Scalar_, int Rows_, int Cols_>
+class Vandermonde : public EigenBase<Vandermonde<Scalar_, Rows_, Cols_>> {
+ public:
+  using Derived = Vandermonde;
+  using StorageBaseType = Vandermonde;
+  using Scalar = Scalar_;
+  using RealScalar = typename NumTraits<Scalar>::Real;
+  using StorageIndex = int;
+  static constexpr int NodeOptions = Rows_ == Dynamic ? AutoAlign : DontAlign;
+  using NodeVector = Matrix<Scalar, Rows_, 1, NodeOptions>;
+  using Nested = Vandermonde;
+
+  static constexpr int RowsAtCompileTime = Rows_;
+  static constexpr int ColsAtCompileTime = Cols_;
+  static constexpr int MaxRowsAtCompileTime = Rows_;
+  static constexpr int MaxColsAtCompileTime = Cols_;
+  static constexpr int SizeAtCompileTime = internal::size_at_compile_time(Rows_, Cols_);
+  static constexpr int MaxSizeAtCompileTime = SizeAtCompileTime;
+  static constexpr int Flags = internal::traits<Vandermonde>::Flags;
+  static constexpr bool IsRowMajor = (Flags & RowMajorBit) != 0;
+  // Deliberately no IsVectorAtCompileTime: Ref<const Vandermonde>'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.
+
+  EIGEN_STATIC_ASSERT_NON_INTEGER(RealScalar)
+  EIGEN_MAKE_SCALAR_BINARY_OP_ONTHELEFT(operator*, internal::scalar_product_op)
+
+  /** Builds an \c m x \a cols Vandermonde matrix from the \c m nodes \a nodes. */
+  template <typename Derived>
+  Vandermonde(const MatrixBase<Derived>& nodes, Index cols) : m_x(nodes), m_cols(cols) {
+    EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+    eigen_assert(m_x.size() > 0 && m_cols > 0 && "Vandermonde must be non-empty");
+    eigen_assert((Cols_ == Dynamic || Cols_ == cols) && "cols does not match the compile-time column count");
+  }
+
+  /** Builds the square Vandermonde matrix of the nodes \a nodes. */
+  template <typename Derived>
+  explicit Vandermonde(const MatrixBase<Derived>& nodes) : Vandermonde(nodes, nodes.size()) {
+    EIGEN_STATIC_ASSERT(Rows_ == Dynamic || Cols_ == Dynamic || Rows_ == Cols_, YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES)
+    EIGEN_STATIC_ASSERT(
+        Cols_ == Dynamic || Derived::SizeAtCompileTime == Dynamic || Cols_ == Derived::SizeAtCompileTime,
+        YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES)
+  }
+
+  EIGEN_DEVICE_FUNC Index rows() const { return m_x.size(); }
+  EIGEN_DEVICE_FUNC Index cols() const { return m_cols; }
+
+  /** \returns the node vector. */
+  const NodeVector& nodes() const { return m_x; }
+
+  /** \returns the coefficient at row \a row and column \a col, \f$ x_i^j \f$. */
+  EIGEN_DEVICE_FUNC Scalar coeff(Index row, Index col) const {
+    Scalar p(1);
+    const Scalar xi = m_x.coeff(row);
+    for (Index t = 0; t < col; ++t) p *= xi;
+    return p;
+  }
+
+  /** \returns the determinant of a \b square Vandermonde matrix through the
+   * closed form \f$ \prod_{i<j} (x_j - x_i) \f$, in O(n^2) operations. The
+   * product is accumulated in the balanced form \c m * 2^e (the split
+   * fraction/exponent determinant convention of LINPACK's xGEDI [3]) -- every
+   * factor and the running product are renormalized to unit magnitude with the
+   * power of two tracked separately, an exact rescaling [2] -- so the partial
+   * products can neither overflow nor underflow when the determinant itself is
+   * representable, whatever the spread of the nodes. Zero factors (repeated
+   * nodes, giving an exactly singular matrix) and non-finite factors propagate
+   * exactly. */
+  Scalar determinant() const {
+    eigen_assert(rows() == cols() && "Vandermonde::determinant requires a square matrix");
+    const Index n = rows();
+    Scalar det(1);
+    internal::structured_exponent_type exponent = 0;
+    for (Index j = 1; j < n; ++j)
+      for (Index i = 0; i < j; ++i) {
+        // A node difference can overflow even though the determinant is
+        // representable (e.g. nodes near +-max); the guarded difference then
+        // enters at half scale with the factor of two carried by the running
+        // exponent.
+        int shift;
+        const Scalar diff = internal::structured_guarded_diff(m_x.coeff(j), m_x.coeff(i), shift);
+        exponent += shift;
+        det = internal::structured_balance(det * internal::structured_balance(diff, exponent), exponent);
+      }
+    // ldexp saturates cleanly to zero / infinity once the accumulated exponent
+    // leaves the representable range; the clamp only guards the narrowing to int.
+    return internal::structured_ldexp_clamped(det, exponent);
+  }
+
+  /** \internal Writes the dense representation into \a dst: column \c j is the
+   * elementwise product of column \c j-1 with the nodes, so only cumulative
+   * columnwise products are involved. Invoked through \c dense = vandermonde; */
+  template <typename Dest>
+  void evalTo(Dest& dst) const {
+    dst.col(0).setOnes();
+    for (Index j = 1; j < m_cols; ++j) dst.col(j) = dst.col(j - 1).cwiseProduct(m_x);
+  }
+
+  /** \internal Computes \c dst += (*this), see evalTo(). */
+  template <typename Dest>
+  void addTo(Dest& dst) const {
+    NodeVector p = NodeVector::Ones(rows());
+    dst.col(0) += p;
+    for (Index j = 1; j < m_cols; ++j) {
+      p = p.cwiseProduct(m_x);
+      dst.col(j) += p;
+    }
+  }
+
+  /** \internal Computes \c dst -= (*this), see evalTo(). */
+  template <typename Dest>
+  void subTo(Dest& dst) const {
+    NodeVector p = NodeVector::Ones(rows());
+    dst.col(0) -= p;
+    for (Index j = 1; j < m_cols; ++j) {
+      p = p.cwiseProduct(m_x);
+      dst.col(j) -= p;
+    }
+  }
+
+  /** \returns the product expression \c (*this) * \a a: the polynomial with
+   * ascending coefficients \c a (per column) evaluated at every node by Horner's
+   * rule, at O(mn) operations and O(1) extra storage. The expression carries the
+   * default product tag, so assigning it behaves like any dense product: a
+   * temporary resolves aliasing between the destination and \a a, and
+   * \c .noalias() skips it. */
+  template <typename Rhs>
+  Product<Vandermonde, Rhs> operator*(const MatrixBase<Rhs>& a) const {
+    EIGEN_STATIC_ASSERT(ColsAtCompileTime == Dynamic || Rhs::RowsAtCompileTime == Dynamic ||
+                            int(ColsAtCompileTime) == int(Rhs::RowsAtCompileTime),
+                        INVALID_MATRIX_PRODUCT)
+    eigen_assert(a.rows() == cols() && "invalid product: dimensions do not match");
+    return Product<Vandermonde, Rhs>(*this, a.derived());
+  }
+
+  /** \internal Computes \c dst += alpha * (*this) * rhs by Horner's rule.
+   * \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.
+   *
+   * Horner intermediates can overflow even when the polynomial value itself is
+   * representable (e.g. coefficients near the overflow threshold evaluated at a
+   * node of magnitude 1/2). Each (node, column) pair is therefore screened with
+   * a conservative exponent bound: when no intermediate can overflow -- every
+   * input of moderate magnitude -- the plain Horner loop runs, matching the naive
+   * evaluation step for step, though not necessarily bit for bit: the compiler
+   * may contract one loop's multiply-add into an FMA and not the other's.
+   * Otherwise scaledHorner() keeps the running value in
+   * the balanced form m * 2^e of determinant(). Non-finite nodes or coefficients
+   * also take the plain loop, which propagates Inf/NaN entrywise like a dense
+   * product; and 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.
+   *
+   * The column's finiteness and its exponent bound come from the single
+   * fast-max pass of internal::structured_exponent_bound_finite(), which is not
+   * guaranteed to propagate NaN. That is sufficient here for the same reason as
+   * in the FFT products: an Inf in NaN-free data always surfaces in a fast
+   * maximum, and a column containing NaN yields the NaN results dense-product
+   * semantics require through the plain loop and the scaled recurrence alike
+   * (every Horner step folds the NaN coefficient in, and the balancing helpers
+   * pass non-finite values through), so missing a NaN cannot change the
+   * result. */
+  template <typename Dest, typename Rhs, typename ProductScalar>
+  void addProduct(Dest& dst, const Rhs& rhs, const ProductScalar& alpha) const {
+    const Index m = rows(), n = m_cols;
+    using Exponent = internal::structured_exponent_type;
+    eigen_assert(rhs.rows() == n && "invalid product: dimensions do not match");
+    const bool unitAlpha = alpha == ProductScalar(1);
+    int log2n = 0;  // n < 2^log2n: bounds the number of addends of the Horner sum
+    for (Index t = n; t > 0; t /= 2) ++log2n;
+    for (Index k = 0; k < rhs.cols(); ++k) {
+      int colExp;  // max modulus < 2^colExp; 0 for a zero or non-finite column
+      const bool colFinite = internal::structured_exponent_bound_finite(rhs.col(k), colExp);
+      for (Index i = 0; i < m; ++i) {
+        const Scalar xi = m_x.coeff(i);
+        // |p_j| < 2^(colExp+log2n+(n-1) max(xiExp,0)+1).
+        const Exponent intermediateBound =
+            Exponent(colExp) + Exponent(log2n) + (Exponent(n) - 1) * Exponent(numext::maxi(exponentBound(xi), 0)) + 2;
+        const bool plain = !colFinite || !(numext::isfinite)(xi) ||
+                           intermediateBound <= Exponent(NumTraits<RealScalar>::max_exponent());
+        ProductScalar acc;
+        if (plain) {
+          acc = rhs.coeff(n - 1, k);
+          for (Index j = n - 2; j >= 0; --j) acc = acc * xi + rhs.coeff(j, k);
+        } else {
+          acc = scaledHorner<ProductScalar>(xi, rhs, k);
+        }
+        dst.coeffRef(i, k) += unitAlpha ? acc : ProductScalar(alpha * acc);
+      }
+    }
+  }
+
+ private:
+  /** \internal \returns an exponent bound \c e with \c |z| < 2^e (the modulus
+   * for a complex \a z), or 0 for a zero or non-finite \a z; the
+   * single-coefficient analogue of internal::structured_exponent_bound(). */
+  template <typename T>
+  static int exponentBound(const T& z) {
+    return exponentBoundImpl(z, internal::bool_constant<NumTraits<T>::IsComplex>());
+  }
+
+  template <typename T>
+  static int exponentBoundImpl(const T& z, std::false_type) {
+    if (!(numext::abs(z) > T(0)) || !(numext::isfinite)(z)) return 0;
+    int e;
+    EIGEN_USING_STD(frexp);
+    frexp(z, &e);
+    return e;
+  }
+
+  template <typename T>
+  static int exponentBoundImpl(const T& z, std::true_type) {
+    using Real = typename NumTraits<T>::Real;
+    const Real mag = numext::maxi(numext::abs(numext::real(z)), numext::abs(numext::imag(z)));
+    if (!(mag > Real(0)) || !(numext::isfinite)(mag)) return 0;
+    int e;
+    EIGEN_USING_STD(frexp);
+    frexp(mag, &e);
+    return e + 1;  // the modulus is at most sqrt(2) times the largest component
+  }
+
+  /** \internal \returns \a z * 2^e computed through two exact half-factors, so
+   * the factors themselves stay representable for the exponent swings the scaled
+   * Horner recurrence produces (|e| up to about twice the scalar's exponent
+   * range on the negative side, at most one exponent range on the positive
+   * side); a factor past the underflow threshold flushes to zero together with
+   * the then-negligible contribution it scales. */
+  template <typename T>
+  static T twoHalfScale(const T& z, internal::structured_exponent_type e) {
+    using Real = typename NumTraits<T>::Real;
+    constexpr internal::structured_exponent_type kMaxExponent = internal::structured_exponent_type(1) << 24;
+    const int ec = static_cast<int>(numext::mini(numext::maxi(e, -kMaxExponent), kMaxExponent));
+    const Real h1 = numext::ldexp(Real(1), ec / 2);
+    const Real h2 = numext::ldexp(Real(1), ec - ec / 2);
+    return (z * h1) * h2;
+  }
+
+  /** \internal Evaluates the polynomial with ascending coefficients
+   * \c rhs.col(k) at the node \a xi, keeping the running value in the balanced
+   * form \c acc * 2^exponent of determinant() (overflow-avoiding power-of-two
+   * rescaling [1], exact by [2]): the node enters through its unit mantissa
+   * with its exponent folded into the running one, the mantissa is
+   * renormalized after every step, and each coefficient is folded into the
+   * running frame scaled by an exact power of two split into two half-factors
+   * (when the coefficient dominates the frame, the frame is rebased onto the
+   * coefficient's exponent instead). Intermediates can therefore neither
+   * overflow nor underflow, and the final ldexp saturates to +-Inf / +-0 exactly
+   * where the true value leaves the representable range.
+   *
+   * An exactly zero mantissa carries no scale, so the frame is reset before
+   * every fold: after an exact cancellation (or a zero node annihilating the
+   * running value) a stale huge frame would otherwise underflow the next small
+   * coefficient to zero. A cancellation that leaves a tiny nonzero mantissa
+   * needs no such care -- the frexp renormalization rebases the frame to the
+   * surviving magnitude (which is a multiple of the operands' unit roundoff,
+   * hence never subnormal for real scalars, and frexp is exact on subnormal
+   * component values regardless).
+   * \pre the node is finite and the column passed the fast-max routing
+   * predicate of addProduct(): it holds no Inf without an accompanying NaN. A
+   * NaN-bearing column can reach the recurrence when the fast maximum misses
+   * the NaN; every helper passes non-finite values through, so it produces the
+   * NaN result dense-product semantics require. */
+  template <typename ProductScalar, typename Rhs>
+  ProductScalar scaledHorner(const Scalar& xi, const Rhs& rhs, Index k) const {
+    using Exponent = internal::structured_exponent_type;
+    Exponent xiE = 0;
+    const Scalar xiMant = internal::structured_balance(xi, xiE);  // xi = xiMant * 2^xiE, exactly
+    ProductScalar acc(0);
+    Exponent exponent = 0;  // running value = acc * 2^exponent
+    for (Index j = m_cols - 1; j >= 0; --j) {
+      if (j < m_cols - 1) {
+        exponent += xiE;
+        acc = internal::structured_balance(acc * xiMant, exponent);
+      }
+      // A zero value has no scale: reset the frame so the next coefficient
+      // enters at its own magnitude instead of underflowing in a stale one.
+      if (acc == ProductScalar(0)) exponent = 0;
+      const ProductScalar aj(rhs.coeff(j, k));
+      if (aj == ProductScalar(0)) continue;
+      const Exponent ajExp = exponentBound(aj);
+      if (exponent < ajExp) {
+        // The coefficient dominates the running frame: rebase onto the
+        // coefficient's exponent. The running value rescales exactly, or
+        // underflows harmlessly once it is negligible against the coefficient.
+        acc = twoHalfScale(acc, exponent - ajExp);
+        exponent = ajExp;
+      }
+      acc = internal::structured_balance(acc + twoHalfScale(aj, -exponent), exponent);
+    }
+    return internal::structured_ldexp_clamped(acc, exponent);
+  }
+
+  NodeVector m_x;
+  Index m_cols;
+};
+
+/** \ingroup StructuredMatrices_Module
+ * \returns an \c m x \a cols \ref Vandermonde operator with node vector \a nodes;
+ * the compile-time row count is deduced from \a nodes. */
+template <typename Derived>
+Vandermonde<typename Derived::Scalar, Derived::SizeAtCompileTime, Dynamic> makeVandermonde(
+    const MatrixBase<Derived>& nodes, Index cols) {
+  return Vandermonde<typename Derived::Scalar, Derived::SizeAtCompileTime, Dynamic>(nodes, cols);
+}
+
+/** \ingroup StructuredMatrices_Module
+ * \returns the square \ref Vandermonde operator of the nodes \a nodes. */
+template <typename Derived>
+Vandermonde<typename Derived::Scalar, Derived::SizeAtCompileTime, Derived::SizeAtCompileTime> makeVandermonde(
+    const MatrixBase<Derived>& nodes) {
+  return Vandermonde<typename Derived::Scalar, Derived::SizeAtCompileTime, Derived::SizeAtCompileTime>(nodes);
+}
+
+/** \ingroup StructuredMatrices_Module
+ * \class BjorckPereyra
+ * \brief Björck-Pereyra O(n^2) solver for square Vandermonde systems.
+ *
+ * Solves \c V*a = f -- polynomial interpolation: find the coefficients of the
+ * polynomial taking values \c f at the nodes -- in O(n^2) operations and O(n)
+ * storage, via divided differences in the Newton basis followed by the basis
+ * change to monomials (Björck & Pereyra, 1970; Golub & Van Loan, Alg. 4.6.2).
+ * The transposed (dual, or moment) system \f$ V^T w = b \f$ is solved by the
+ * companion dual recurrences through the standard \c SolverBase idiom:
+ * \code
+ *   BjorckPereyra<double> bp(V);            // or bp.compute(V);
+ *   VectorXd a = bp.solve(f);                // solve V   * a = f
+ *   VectorXd w = bp.transpose().solve(b);    // solve V^T * w = b
+ *   VectorXd u = bp.adjoint().solve(b);      // solve V^H * u = b
+ * \endcode
+ *
+ * There is no factorization: \c compute() stores the nodes (and flags exactly
+ * repeated or non-finite nodes through \c info()), and each solve runs the
+ * O(n^2) recurrences directly. Genuinely complex node sets are put in a
+ * deterministic Leja order to control growth in the Newton representation;
+ * real nodes, including complex scalars with zero imaginary parts, retain their
+ * input order and its useful monotonicity properties.
+ *
+ * Despite the exponential conditioning of real-node Vandermonde matrices, the
+ * computed solution is often far more accurate than the conditioning suggests:
+ * for monotonically ordered nodes and a right-hand side with alternating signs
+ * the forward error is governed by a small relative-perturbation bound
+ * independent of the condition number (Higham, ASNA ch. 22).
+ *
+ * \tparam Scalar_ a floating-point-like real or complex scalar supporting
+ * Eigen's scalar math hooks, including \c isfinite (and \c abs and \c log for
+ * complex node ordering). Integer types are rejected.
+ *
+ * \sa class Vandermonde
+ */
+template <typename Scalar_>
+class BjorckPereyra : public SolverBase<BjorckPereyra<Scalar_>> {
+ public:
+  using Base = SolverBase<BjorckPereyra>;
+  friend class SolverBase<BjorckPereyra>;
+  EIGEN_GENERIC_PUBLIC_INTERFACE(BjorckPereyra)
+  EIGEN_STATIC_ASSERT_NON_INTEGER(RealScalar)
+  using NodeVector = Matrix<Scalar, Dynamic, 1>;
+
+  /** Default constructor; call \ref compute before \ref solve. */
+  BjorckPereyra() : m_isInitialized(false), m_info(InvalidInput) {}
+
+  /** Constructs the solver for the square Vandermonde matrix \a V. */
+  template <int Rows_, int Cols_>
+  explicit BjorckPereyra(const Vandermonde<Scalar, Rows_, Cols_>& V) : m_isInitialized(false), m_info(InvalidInput) {
+    compute(V);
+  }
+
+  /** Stores the nodes of the square Vandermonde matrix \a V, checks that they
+   * are finite and distinct, and computes a Leja order for genuinely complex
+   * nodes. */
+  template <int Rows_, int Cols_>
+  BjorckPereyra& compute(const Vandermonde<Scalar, Rows_, Cols_>& V) {
+    EIGEN_STATIC_ASSERT(Rows_ == Dynamic || Cols_ == Dynamic || Rows_ == Cols_, YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES)
+    eigen_assert(V.rows() == V.cols() && "BjorckPereyra requires a square Vandermonde matrix");
+    m_x = V.nodes();
+    m_order.clear();
+    m_info = Success;
+    const Index n = m_x.size();
+    if (!m_x.allFinite()) m_info = InvalidInput;
+    for (Index j = 1; j < n && m_info == Success; ++j)
+      for (Index i = 0; i < j; ++i) {
+        if (m_x[i] == m_x[j]) {
+          m_info = NumericalIssue;
+          break;
+        }
+      }
+    if (m_info == Success) initializeNodeOrder(m_x, internal::bool_constant<NumTraits<Scalar>::IsComplex>());
+    m_isInitialized = true;
+    return *this;
+  }
+
+  Index rows() const noexcept { return m_x.size(); }
+  Index cols() const noexcept { return m_x.size(); }
+
+  /** \returns \c Success, \c NumericalIssue when the nodes contain an exact
+   * duplicate (the matrix is singular), or \c InvalidInput for a non-finite
+   * node. */
+  ComputationInfo info() const {
+    eigen_assert(m_isInitialized && "BjorckPereyra is not initialized.");
+    return m_info;
+  }
+
+#ifdef EIGEN_PARSED_BY_DOXYGEN
+  /** \returns the solution \c a of \c V*a = \a f, as a lazily evaluated
+   * expression. Supports multiple right-hand sides. The transposed and adjoint
+   * systems are solved through \c transpose().solve(b) and \c adjoint().solve(b).
+   * \pre \ref compute has been called. */
+  template <typename Rhs>
+  inline const Solve<BjorckPereyra, Rhs> solve(const MatrixBase<Rhs>& f) const;
+#endif
+
+#ifndef EIGEN_PARSED_BY_DOXYGEN
+  /** \internal Primal solve V*a = rhs: divided differences (Newton coefficients),
+   * then the Newton-to-monomial basis change. */
+  template <typename RhsType, typename DstType>
+  void _solve_impl(const RhsType& rhs, DstType& dst) const {
+    using RhsScalar = typename RhsType::Scalar;
+    using WorkScalar = typename DstType::Scalar;
+    using ProductOp = internal::scalar_product_op<Scalar, RhsScalar>;
+    EIGEN_CHECK_BINARY_COMPATIBILITY(ProductOp, Scalar, RhsScalar)
+
+    const Index n = m_x.size();
+    dst = rhs;
+    Matrix<WorkScalar, Dynamic, 1> permuted;
+    if (!m_order.empty()) permuted.resize(n);
+    for (Index k = 0; k < rhs.cols(); ++k) {
+      auto a = dst.col(k);
+      if (!m_order.empty()) {
+        for (Index i = 0; i < n; ++i) permuted[i] = a[m_order[static_cast<std::size_t>(i)]];
+        a = permuted;
+      }
+      for (Index j = 0; j < n - 1; ++j)
+        for (Index i = n - 1; i > j; --i) a[i] = (a[i] - a[i - 1]) / (m_x[i] - m_x[i - j - 1]);
+      for (Index j = n - 2; j >= 0; --j)
+        for (Index i = j; i < n - 1; ++i) a[i] -= m_x[j] * a[i + 1];
+    }
+  }
+
+  /** \internal Transposed (dual) solve V^T*w = rhs: the transposes of the primal
+   * elementary steps, applied in reverse order; conjugated on the way in and out
+   * for the adjoint. */
+  template <bool Conjugate, typename RhsType, typename DstType>
+  void _solve_impl_transposed(const RhsType& rhs, DstType& dst) const {
+    using RhsScalar = typename RhsType::Scalar;
+    using WorkScalar = typename DstType::Scalar;
+    using ProductOp = internal::scalar_product_op<Scalar, RhsScalar>;
+    EIGEN_CHECK_BINARY_COMPATIBILITY(ProductOp, Scalar, RhsScalar)
+
+    const Index n = m_x.size();
+    dst = rhs.template conjugateIf<Conjugate>();
+    Matrix<WorkScalar, Dynamic, 1> permuted;
+    if (!m_order.empty()) permuted.resize(n);
+    for (Index k = 0; k < rhs.cols(); ++k) {
+      auto w = dst.col(k);
+      for (Index j = 0; j < n - 1; ++j)
+        for (Index i = n - 1; i > j; --i) w[i] -= m_x[j] * w[i - 1];
+      for (Index j = n - 2; j >= 0; --j) {
+        w.tail(n - j - 1).array() /= (m_x.tail(n - j - 1) - m_x.head(n - j - 1)).array();
+        for (Index i = j; i < n - 1; ++i) w[i] -= w[i + 1];
+      }
+
+      if (!m_order.empty()) {
+        permuted = w;
+        for (Index i = 0; i < n; ++i) w[m_order[static_cast<std::size_t>(i)]] = permuted[i];
+      }
+    }
+    if (Conjugate) dst = dst.conjugate().eval();
+  }
+#endif
+
+ private:
+  static RealScalar lejaLogAbs(const Scalar& z) {
+    const RealScalar re = numext::abs(numext::real(z));
+    const RealScalar im = numext::abs(numext::imag(z));
+    const RealScalar scale = numext::maxi(re, im);
+    if (scale == RealScalar(0)) return -NumTraits<RealScalar>::infinity();
+    const RealScalar scaledRe = re / scale;
+    const RealScalar scaledIm = im / scale;
+    return numext::log(scale) + RealScalar(0.5) * numext::log(scaledRe * scaledRe + scaledIm * scaledIm);
+  }
+
+  static RealScalar lejaLogDistance(const Scalar& a, const Scalar& b) {
+    int exponent;
+    const Scalar difference = internal::structured_guarded_diff(a, b, exponent);
+    return lejaLogAbs(difference) + RealScalar(exponent) * numext::log(RealScalar(2));
+  }
+
+  void initializeNodeOrder(const NodeVector&, std::false_type) {}
+
+  void initializeNodeOrder(const NodeVector& nodes, std::true_type) {
+    using Real = typename NumTraits<Scalar>::Real;
+    const Index n = nodes.size();
+    bool genuinelyComplex = false;
+    for (Index i = 0; i < n; ++i) genuinelyComplex = genuinelyComplex || numext::imag(nodes[i]) != Real(0);
+    if (!genuinelyComplex || n < 2) return;
+
+    const NodeVector original = nodes;
+
+    m_order.resize(static_cast<std::size_t>(n));
+    std::vector<char> selected(static_cast<std::size_t>(n), 0);
+    std::vector<RealScalar> scores(static_cast<std::size_t>(n), RealScalar(0));
+    Index next;  // the Leja sequence starts at a node of largest modulus
+    original.unaryExpr(&lejaLogAbs).maxCoeff(&next);
+
+    for (Index position = 0; position < n; ++position) {
+      m_order[static_cast<std::size_t>(position)] = next;
+      selected[static_cast<std::size_t>(next)] = 1;
+      if (position + 1 == n) break;
+
+      Index candidate = -1;
+      RealScalar candidateScore = -NumTraits<RealScalar>::infinity();
+      for (Index i = 0; i < n; ++i) {
+        if (selected[static_cast<std::size_t>(i)]) continue;
+        scores[static_cast<std::size_t>(i)] += lejaLogDistance(original[i], original[next]);
+        if (candidate < 0 || scores[static_cast<std::size_t>(i)] > candidateScore) {
+          candidate = i;
+          candidateScore = scores[static_cast<std::size_t>(i)];
+        }
+      }
+      next = candidate;
+    }
+
+    for (Index i = 0; i < n; ++i) m_x[i] = original[m_order[static_cast<std::size_t>(i)]];
+  }
+
+  NodeVector m_x;
+  std::vector<Index> m_order;
+  bool m_isInitialized;
+  ComputationInfo m_info;
+};
+
+namespace internal {
+
+/** \internal Solve results use the scalar promoted from the solver and RHS.
+ * Core's generic solve traits retain the RHS scalar, which would discard the
+ * imaginary part when a complex Vandermonde is applied to a real RHS. */
+template <typename SolverScalar, typename RhsScalar,
+          bool Compatible = has_ReturnType<ScalarBinaryOpTraits<SolverScalar, RhsScalar>>::value>
+struct bjorck_pereyra_result_scalar {
+  using type = SolverScalar;
+};
+
+template <typename SolverScalar, typename RhsScalar>
+struct bjorck_pereyra_result_scalar<SolverScalar, RhsScalar, true> {
+  using type = typename ScalarBinaryOpTraits<SolverScalar, RhsScalar>::ReturnType;
+};
+
+template <typename SolverScalar, typename RhsType>
+struct bjorck_pereyra_solve_traits {
+  using ResultScalar = typename bjorck_pereyra_result_scalar<SolverScalar, typename RhsType::Scalar>::type;
+  using PlainObject =
+      typename make_proper_matrix_type<ResultScalar, Dynamic, RhsType::ColsAtCompileTime, RhsType::PlainObject::Options,
+                                       Dynamic, RhsType::MaxColsAtCompileTime>::type;
+};
+
+template <typename Scalar_, typename RhsType>
+struct solve_traits<BjorckPereyra<Scalar_>, RhsType, Dense> : bjorck_pereyra_solve_traits<Scalar_, RhsType> {};
+
+template <typename Scalar_, typename RhsType>
+struct solve_traits<Transpose<const BjorckPereyra<Scalar_>>, RhsType, Dense>
+    : bjorck_pereyra_solve_traits<Scalar_, RhsType> {};
+
+template <typename Scalar_, typename RhsType>
+struct solve_traits<CwiseUnaryOp<scalar_conjugate_op<Scalar_>, const Transpose<const BjorckPereyra<Scalar_>>>, RhsType,
+                    Dense> : bjorck_pereyra_solve_traits<Scalar_, RhsType> {};
+
+template <typename Factor, typename Scalar_, int Rows_, int Cols_, typename Plain, typename Rhs>
+struct scaled_vandermonde_product_impl
+    : generic_product_impl_base<
+          CwiseBinaryOp<scalar_product_op<Factor, Scalar_>, const CwiseNullaryOp<scalar_constant_op<Factor>, Plain>,
+                        const Vandermonde<Scalar_, Rows_, Cols_>>,
+          Rhs, scaled_vandermonde_product_impl<Factor, Scalar_, Rows_, Cols_, Plain, Rhs>> {
+  using Op = Vandermonde<Scalar_, Rows_, Cols_>;
+  using ScaledOp = CwiseBinaryOp<scalar_product_op<Factor, Scalar_>,
+                                 const CwiseNullaryOp<scalar_constant_op<Factor>, Plain>, const Op>;
+  using Scalar = typename Product<ScaledOp, Rhs>::Scalar;
+
+  template <typename Dest>
+  static void scaleAndAddTo(Dest& dst, const ScaledOp& lhs, const Rhs& rhs, const Scalar& alpha) {
+    using RhsNested = typename nested_eval<Rhs, Rows_>::type;
+    RhsNested actualRhs(rhs);
+    lhs.rhs().addProduct(dst, actualRhs, Scalar(alpha * lhs.lhs().functor().m_other));
+  }
+};
+
+// Preserve the Horner kernel after Core introduces the scaled wrapper above;
+// otherwise the wrapper has DenseShape and falls back to a coefficient product.
+#define EIGEN_SCALED_VANDERMONDE_PRODUCT_IMPL(ProductTag)                                                        \
+  template <typename Factor, typename Scalar_, int Rows_, int Cols_, typename Plain, typename Rhs>               \
+  struct generic_product_impl<                                                                                   \
+      CwiseBinaryOp<scalar_product_op<Factor, Scalar_>, const CwiseNullaryOp<scalar_constant_op<Factor>, Plain>, \
+                    const Vandermonde<Scalar_, Rows_, Cols_>>,                                                   \
+      Rhs, DenseShape, DenseShape, ProductTag>                                                                   \
+      : scaled_vandermonde_product_impl<Factor, Scalar_, Rows_, Cols_, Plain, Rhs> {};
+
+EIGEN_SCALED_VANDERMONDE_PRODUCT_IMPL(CoeffBasedProductMode)
+EIGEN_SCALED_VANDERMONDE_PRODUCT_IMPL(LazyCoeffBasedProductMode)
+EIGEN_SCALED_VANDERMONDE_PRODUCT_IMPL(OuterProduct)
+EIGEN_SCALED_VANDERMONDE_PRODUCT_IMPL(InnerProduct)
+EIGEN_SCALED_VANDERMONDE_PRODUCT_IMPL(GemvProduct)
+EIGEN_SCALED_VANDERMONDE_PRODUCT_IMPL(GemmProduct)
+
+#undef EIGEN_SCALED_VANDERMONDE_PRODUCT_IMPL
+
+template <typename Scalar_, int Rows_, int Cols_, typename Rhs, int ProductTag>
+struct generic_product_impl<Vandermonde<Scalar_, Rows_, Cols_>, Rhs, StructuredShape, DenseShape, ProductTag>
+    : structured_product_impl<Vandermonde<Scalar_, Rows_, Cols_>, Rhs> {};
+
+}  // namespace internal
+
+}  // namespace Eigen
+
+#endif  // EIGEN_STRUCTURED_VANDERMONDE_H
diff --git a/unsupported/benchmarks/StructuredMatrices/CMakeLists.txt b/unsupported/benchmarks/StructuredMatrices/CMakeLists.txt
index f87f4ea..26c2428 100644
--- a/unsupported/benchmarks/StructuredMatrices/CMakeLists.txt
+++ b/unsupported/benchmarks/StructuredMatrices/CMakeLists.txt
@@ -8,3 +8,4 @@
 eigen_add_benchmark(bench_structured_transpose bench_structured_transpose.cpp)
 eigen_add_benchmark(bench_structured_kronecker bench_structured_kronecker.cpp)
 eigen_add_benchmark(bench_structured_dplr bench_structured_dplr.cpp)
+eigen_add_benchmark(bench_structured_vandermonde bench_structured_vandermonde.cpp)
diff --git a/unsupported/benchmarks/StructuredMatrices/bench_structured_vandermonde.cpp b/unsupported/benchmarks/StructuredMatrices/bench_structured_vandermonde.cpp
new file mode 100644
index 0000000..05c7b9d
--- /dev/null
+++ b/unsupported/benchmarks/StructuredMatrices/bench_structured_vandermonde.cpp
@@ -0,0 +1,148 @@
+// Benchmarks for the Vandermonde operator: the O(mn)-flops / O(m)-storage
+// Horner product against the equivalent dense GEMV on the materialized matrix
+// (comparable O(mn) arithmetic, O(mn) storage), and the O(n^2)
+// Björck-Pereyra square solve against a dense PartialPivLU factor-and-solve
+// (O(n^3)). These cases cover nominal real-input performance. Both solve
+// variants time the full pipeline from the stored representation to the
+// solution; the dense one is additionally handed the materialized matrix for
+// free (built outside the loop).
+// SPDX-FileCopyrightText: The Eigen Authors
+// SPDX-License-Identifier: MPL-2.0
+
+#include <benchmark/benchmark.h>
+#include <Eigen/Core>
+#include <Eigen/LU>
+#include <unsupported/Eigen/StructuredMatrices>
+
+using namespace Eigen;
+
+typedef Matrix<double, Dynamic, 1> Vec;
+typedef Matrix<double, Dynamic, Dynamic> Mat;
+
+// Chebyshev nodes in (-1, 1): deterministic, bounded, and distinct (no
+// Björck-Pereyra division by zero or exact-zero LU pivot).
+static Vec chebyshevNodes(Index n) {
+  Vec x(n);
+  for (Index i = 0; i < n; ++i) x[i] = std::cos(EIGEN_PI * double(2 * i + 1) / double(2 * n));
+  return x;
+}
+
+static Vec polynomialCoefficients(Index n) {
+  Vec a(n);
+  for (Index i = 0; i < n; ++i) a[i] = std::cos(0.25 * double(i + 1)) / double(i + 1);
+  return a;
+}
+
+// Long-double Horner reference used only before timing. The dimension-scaled
+// bound accounts for at most n multiply-adds in each polynomial evaluation.
+static bool validateProduct(benchmark::State& state, const Vec& x, const Vec& a, const Vec& actual) {
+  Vec expected(x.size());
+  for (Index i = 0; i < x.size(); ++i) {
+    long double acc = static_cast<long double>(a[a.size() - 1]);
+    for (Index j = a.size() - 2; j >= 0; --j)
+      acc = acc * static_cast<long double>(x[i]) + static_cast<long double>(a[j]);
+    expected[i] = static_cast<double>(acc);
+  }
+  const double tol = 64.0 * double(a.size()) * NumTraits<double>::epsilon();
+  if (!actual.allFinite() || (actual - expected).norm() > tol * numext::maxi(1.0, expected.norm())) {
+    state.SkipWithError("Vandermonde product failed validation");
+    return false;
+  }
+  return true;
+}
+
+// The large Chebyshev-node systems are extremely ill-conditioned, so a random
+// right-hand side has no meaningful double-precision monomial solution. The
+// exact constant polynomial f_i=1 gives both algorithms the finite,
+// independently known solution a=e_0 and keeps this a reproducible nominal-path
+// benchmark. Validate the solution once before entering either timed loop.
+static bool validateConstantSolve(benchmark::State& state, const Vec& actual) {
+  Vec expected = Vec::Zero(actual.size());
+  expected[0] = 1.0;
+  const double tol = 64.0 * double(actual.size()) * NumTraits<double>::epsilon();
+  if (!actual.allFinite() || (actual - expected).norm() > tol) {
+    state.SkipWithError("Vandermonde solve failed validation");
+    return false;
+  }
+  return true;
+}
+
+// --- Product y = V * a: Horner rule on the nodes ---
+static void BM_VandermondeProduct(benchmark::State& state) {
+  const Index m = state.range(0), n = state.range(1);
+  Vec x = chebyshevNodes(m), a = polynomialCoefficients(n), y(m);
+  Vandermonde<double> V(x, n);
+  Vec check = V * a;
+  if (!validateProduct(state, x, a, check)) return;
+  for (auto _ : state) {
+    y.noalias() = V * a;
+    benchmark::DoNotOptimize(y.data());
+    benchmark::ClobberMemory();
+  }
+}
+BENCHMARK(BM_VandermondeProduct)
+    ->Args({64, 64})
+    ->Args({256, 256})
+    ->Args({1024, 1024})
+    ->Args({4096, 4096})
+    ->Args({4096, 256})
+    ->Args({256, 4096})
+    ->ArgNames({"m", "n"});
+
+// --- Product y = dense * a: GEMV on the pre-materialized matrix ---
+static void BM_VandermondeProductDense(benchmark::State& state) {
+  const Index m = state.range(0), n = state.range(1);
+  Vec x = chebyshevNodes(m), a = polynomialCoefficients(n), y(m);
+  Vandermonde<double> V(x, n);
+  Mat dense = V;  // materialized once, outside the timed loop
+  Vec check = dense * a;
+  if (!validateProduct(state, x, a, check)) return;
+  for (auto _ : state) {
+    y.noalias() = dense * a;
+    benchmark::DoNotOptimize(y.data());
+    benchmark::ClobberMemory();
+  }
+}
+BENCHMARK(BM_VandermondeProductDense)
+    ->Args({64, 64})
+    ->Args({256, 256})
+    ->Args({1024, 1024})
+    ->Args({4096, 4096})
+    ->Args({4096, 256})
+    ->Args({256, 4096})
+    ->ArgNames({"m", "n"});
+
+// --- Square solve V * a = f: Björck-Pereyra divided-difference recurrences ---
+static void BM_BjorckPereyraSolve(benchmark::State& state) {
+  const Index n = state.range(0);
+  Vec x = chebyshevNodes(n), f = Vec::Ones(n), a(n);
+  Vandermonde<double> V(x);
+  BjorckPereyra<double> bpCheck(V);
+  Vec check = bpCheck.solve(f);
+  if (!validateConstantSolve(state, check)) return;
+  for (auto _ : state) {
+    BjorckPereyra<double> bp(V);  // compute() includes the O(n^2) duplicate scan
+    a = bp.solve(f);
+    benchmark::DoNotOptimize(a.data());
+    benchmark::ClobberMemory();
+  }
+}
+BENCHMARK(BM_BjorckPereyraSolve)->Arg(16)->Arg(64)->Arg(256)->Arg(1024);
+
+// --- Square solve dense * a = f: PartialPivLU on the pre-materialized matrix ---
+static void BM_VandermondeSolveDense(benchmark::State& state) {
+  const Index n = state.range(0);
+  Vec x = chebyshevNodes(n), f = Vec::Ones(n), a(n);
+  Vandermonde<double> V(x);
+  Mat dense = V;  // materialized once, outside the timed loop
+  PartialPivLU<Mat> luCheck(dense);
+  Vec check = luCheck.solve(f);
+  if (!validateConstantSolve(state, check)) return;
+  for (auto _ : state) {
+    PartialPivLU<Mat> lu(dense);  // the O(n^3) elimination
+    a = lu.solve(f);
+    benchmark::DoNotOptimize(a.data());
+    benchmark::ClobberMemory();
+  }
+}
+BENCHMARK(BM_VandermondeSolveDense)->Arg(16)->Arg(64)->Arg(256)->Arg(1024);
diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt
index a84ab5c..17078a0 100644
--- a/unsupported/test/CMakeLists.txt
+++ b/unsupported/test/CMakeLists.txt
@@ -125,6 +125,8 @@
     ei_add_property(EIGEN_MISSING_BACKENDS "OpenGL, ")
 endif()
 
+ei_add_test(structured_vandermonde)
+ei_add_test(structured_vandermonde_int_index)
 ei_add_test(polynomialsolver)
 ei_add_test(polynomialutils)
 ei_add_test(splines)
diff --git a/unsupported/test/structured_vandermonde.cpp b/unsupported/test/structured_vandermonde.cpp
new file mode 100644
index 0000000..43c3648
--- /dev/null
+++ b/unsupported/test/structured_vandermonde.cpp
@@ -0,0 +1,816 @@
+// 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
+
+#include "main.h"
+
+#include <unsupported/Eigen/StructuredMatrices>
+
+using namespace Eigen;
+
+// Reference dense Vandermonde built entry-wise from the nodes.
+template <typename Scalar>
+Matrix<Scalar, Dynamic, Dynamic> reference_vandermonde(const Matrix<Scalar, Dynamic, 1>& x, Index n) {
+  const Index m = x.size();
+  Matrix<Scalar, Dynamic, Dynamic> dense(m, n);
+  for (Index i = 0; i < m; ++i) {
+    Scalar p(1);
+    for (Index j = 0; j < n; ++j) {
+      dense(i, j) = p;
+      p *= x[i];
+    }
+  }
+  return dense;
+}
+
+template <typename Scalar>
+void test_vandermonde_product(Index m, Index n) {
+  typedef Matrix<Scalar, Dynamic, 1> Vec;
+  typedef Matrix<Scalar, Dynamic, Dynamic> Mat;
+
+  Vec x = Vec::Random(m);
+  Vandermonde<Scalar> V(x, n);
+  Mat dense = reference_vandermonde<Scalar>(x, n);
+
+  Mat Vd = V;
+  VERIFY_IS_APPROX(Vd, dense);
+  for (Index t = 0; t < 5; ++t) {
+    Index i = internal::random<Index>(0, m - 1), j = internal::random<Index>(0, n - 1);
+    VERIFY_IS_APPROX(V.coeff(i, j), dense(i, j));
+  }
+
+  Vec a = Vec::Random(n);
+  VERIFY_IS_APPROX((V * a).eval(), (dense * a).eval());
+
+  Mat A = Mat::Random(n, 3);
+  VERIFY_IS_APPROX((V * A).eval(), (dense * A).eval());
+
+  // Accumulation form exercised by the iterative solvers.
+  Vec y = Vec::Random(m);
+  Vec y0 = y;
+  y.noalias() += V * a;
+  VERIFY_IS_APPROX(y, (y0 + dense * a).eval());
+}
+
+// Well-conditioned interpolation: Chebyshev-like nodes in [-1,1] keep the
+// conditioning around (1+sqrt(2))^n, so for moderate n a modest tolerance wide
+// of roundoff verifies the primal and dual solves against the exact solution.
+template <typename Scalar>
+void test_bjorck_pereyra(Index n) {
+  typedef typename NumTraits<Scalar>::Real RealScalar;
+  typedef Matrix<Scalar, Dynamic, 1> Vec;
+  typedef Matrix<Scalar, Dynamic, Dynamic> Mat;
+
+  Vec x(n);
+  for (Index i = 0; i < n; ++i)
+    x[i] = Scalar(RealScalar(std::cos(double(EIGEN_PI) * double(2 * i + 1) / double(2 * n))));  // Chebyshev nodes
+  Vandermonde<Scalar> V(x);
+  Mat dense = reference_vandermonde<Scalar>(x, n);
+
+  BjorckPereyra<Scalar> bp(V);
+  VERIFY(bp.info() == Success);
+
+  // Forward-error bound for the (moderately conditioned) Chebyshev-node system.
+  const RealScalar tol = RealScalar(5e7) * NumTraits<RealScalar>::epsilon();  // ~1e-8 in double
+
+  // Primal: interpolate values of a known polynomial and recover its coefficients.
+  Vec aTrue = Vec::Random(n);
+  Vec f = dense * aTrue;
+  Vec a = bp.solve(f);
+  VERIFY((a - aTrue).norm() <= tol * aTrue.norm());
+
+  // Dual (moment) system through the SolverBase transpose idiom.
+  Vec wTrue = Vec::Random(n);
+  Vec b = dense.transpose() * wTrue;
+  Vec w = bp.transpose().solve(b);
+  VERIFY((w - wTrue).norm() <= tol * wTrue.norm());
+
+  // Adjoint solve.
+  Vec c = dense.adjoint() * wTrue;
+  Vec u = bp.adjoint().solve(c);
+  VERIFY((u - wTrue).norm() <= tol * wTrue.norm());
+
+  // Multiple right-hand sides.
+  Mat F = Mat::Random(n, 3);
+  Mat A = bp.solve(F);
+  VERIFY_IS_APPROX((dense * A).eval(), F);
+}
+
+// The n-th roots of unity make V/sqrt(n) exactly unitary. Verify against the
+// analytic inverse a = V^H f / n to catch catastrophic Newton-basis growth.
+template <typename RealScalar>
+void test_bjorck_pereyra_roots_of_unity(Index n) {
+  typedef std::complex<RealScalar> Scalar;
+  typedef Matrix<Scalar, Dynamic, 1> Vec;
+
+  Vec x(n);
+  for (Index i = 0; i < n; ++i)
+    x[i] = std::polar(RealScalar(1), RealScalar(2 * EIGEN_PI) * RealScalar(i) / RealScalar(n));
+  Vandermonde<Scalar> V(x);
+
+  BjorckPereyra<Scalar> bp(V);
+  VERIFY(bp.info() == Success);
+
+  // The first cardinal polynomial has every coefficient equal to 1/n. This
+  // deterministic case exposes the severe Newton-basis growth of cyclic node
+  // order while remaining analytically checkable.
+  Vec f = Vec::Zero(n);
+  f[0] = Scalar(1);
+  Vec a = bp.solve(f);
+  Vec aRef = Vec::Constant(n, Scalar(RealScalar(1) / RealScalar(n)));
+  // The recurrence performs O(n^2) complex operations; the fixed factor also
+  // covers growth in its intermediate Newton coefficients.
+  const RealScalar tol = RealScalar(128) * RealScalar(n) * RealScalar(n) * NumTraits<RealScalar>::epsilon();
+  VERIFY((a - aRef).norm() <= tol * aRef.norm());
+}
+
+// Eigen::half has explicit conversion from the wider type returned by the
+// standard scaling functions. Instantiate both product and solver paths so the
+// range-protection helpers retain support for narrow floating-point scalars.
+void test_vandermonde_half() {
+  typedef Matrix<half, 2, 1> Vec;
+  Vec nodes, coefficients;
+  nodes << half(0), half(2);
+  coefficients << half(3), half(2);
+
+  Vandermonde<half, 2, 2> V(nodes);
+  Vec values = V * coefficients;
+  VERIFY_IS_EQUAL(float(values[0]), 3.0f);
+  VERIFY_IS_EQUAL(float(values[1]), 7.0f);
+
+  BjorckPereyra<half> bp(V);
+  Vec recovered = bp.solve(values);
+  VERIFY_IS_EQUAL(float(recovered[0]), 3.0f);
+  VERIFY_IS_EQUAL(float(recovered[1]), 2.0f);
+}
+
+// Björck-Pereyra's celebrated accuracy property (Higham, ASNA ch. 22): for
+// monotone nodes and an alternating-sign right-hand side the forward error is
+// tiny even though the matrix conditioning is astronomical. Deterministic.
+void test_bjorck_pereyra_higham() {
+  typedef Matrix<double, Dynamic, 1> Vec;
+
+  const Index n = 20;
+  Vec x(n);
+  for (Index i = 0; i < n; ++i) x[i] = double(i + 1) / 32.0;  // exactly represented, monotone in (0,1)
+  Vandermonde<double> V(x);
+
+  // f_i = (-1)^i: alternating signs.
+  Vec f(n);
+  for (Index i = 0; i < n; ++i) f[i] = (i % 2 == 0) ? 1.0 : -1.0;
+  BjorckPereyra<double> bp(V);
+  Vec a = bp.solve(f);
+
+  // Independent long-double coefficient reference. For equally spaced nodes
+  // with h=1/32, the Newton weights are Delta^k f_0 / (k! h^k) = (-64)^k/k!.
+  // Accumulate weight_k * product_{j<k}(z-x_j) directly in the monomial basis.
+  Matrix<long double, Dynamic, 1> aRef = Matrix<long double, Dynamic, 1>::Zero(n);
+  Matrix<long double, Dynamic, 1> basis = Matrix<long double, Dynamic, 1>::Zero(n);
+  basis[0] = 1.0L;
+  long double weight = 1.0L;
+  for (Index k = 0; k < n; ++k) {
+    aRef += weight * basis;
+    if (k + 1 == n) break;
+    const long double node = static_cast<long double>(x[k]);
+    basis[k + 1] = basis[k];
+    for (Index j = k; j > 0; --j) basis[j] = basis[j - 1] - node * basis[j];
+    basis[0] *= -node;
+    weight *= -64.0L / static_cast<long double>(k + 1);
+  }
+
+  const Vec roundedRef = aRef.cast<double>();
+  const double maxRelativeError = ((a - roundedRef).cwiseAbs().array() / roundedRef.cwiseAbs().array()).maxCoeff();
+  const double coefficientTol = 8.0 * double(n) * double(n) * NumTraits<double>::epsilon();
+  VERIFY(maxRelativeError <= coefficientTol);
+}
+
+// The transpose and adjoint recurrences must be distinguished for genuinely
+// complex nodes. A rotated root-of-unity grid keeps both systems well
+// conditioned while avoiding the unrotated grid's extra symmetries.
+void test_bjorck_pereyra_complex_transpose_adjoint(Index n) {
+  typedef std::complex<double> Scalar;
+  typedef Matrix<double, Dynamic, 1> RealVec;
+  typedef Matrix<Scalar, Dynamic, 1> Vec;
+  typedef Matrix<Scalar, Dynamic, Dynamic> Mat;
+
+  Vec x(n);
+  for (Index i = 0; i < n; ++i) x[i] = std::polar(1.0, double(2 * EIGEN_PI) * (double(i) + 0.25) / double(n));
+  Vandermonde<Scalar> V(x);
+  Mat dense = reference_vandermonde<Scalar>(x, n);
+  BjorckPereyra<Scalar> bp(V);
+
+  Vec expected = Vec::Random(n);
+  Vec transposed = bp.transpose().solve(dense.transpose() * expected);
+  Vec adjoint = bp.adjoint().solve(dense.adjoint() * expected);
+  const double tol = 256.0 * double(n) * double(n) * NumTraits<double>::epsilon();
+  const double scale = numext::maxi(1.0, expected.norm());
+  VERIFY((transposed - expected).norm() <= tol * scale);
+  VERIFY((adjoint - expected).norm() <= tol * scale);
+
+  // A real right-hand side is compatible with complex nodes and promotes the
+  // solve result to complex. The rotated Fourier matrix has V^H V = n I, so
+  // all three analytic inverse forms are available without another solver.
+  RealVec realRhs = RealVec::Random(n);
+  Vec primalMixed = bp.solve(realRhs);
+  Vec transposedMixed = bp.transpose().solve(realRhs);
+  Vec adjointMixed = bp.adjoint().solve(realRhs);
+  Vec primalRef = dense.adjoint() * realRhs / Scalar(double(n));
+  Vec transposedRef = dense.conjugate() * realRhs / Scalar(double(n));
+  Vec adjointRef = dense * realRhs / Scalar(double(n));
+  const double mixedScale = numext::maxi(1.0, realRhs.norm());
+  VERIFY((primalMixed - primalRef).norm() <= tol * mixedScale);
+  VERIFY((transposedMixed - transposedRef).norm() <= tol * mixedScale);
+  VERIFY((adjointMixed - adjointRef).norm() <= tol * mixedScale);
+
+  // The reverse promotion direction remains valid as well: real nodes with a
+  // complex right-hand side produce complex work values and output.
+  RealVec realNodes(2);
+  realNodes << -1.0, 1.0;
+  BjorckPereyra<double> realBp{Vandermonde<double>(realNodes)};
+  Vec complexRhs(2), complexExpected(2);
+  complexRhs << Scalar(1.0, 1.0), Scalar(3.0, -2.0);
+  complexExpected << Scalar(2.0, -0.5), Scalar(1.0, -1.5);
+  Vec realNodeMixed = realBp.solve(complexRhs);
+  VERIFY((realNodeMixed - complexExpected).norm() <= 4.0 * NumTraits<double>::epsilon() * complexExpected.norm());
+
+  // Materializing an aliased RHS must precede the internal Leja permutation.
+  Vec aliased = Vec::Random(n);
+  const Vec aliasInput = (aliased + Vec::Ones(n)).eval();
+  const Vec aliasReference = bp.solve(aliasInput);
+  aliased = bp.solve(aliased + Vec::Ones(n));
+  VERIFY((aliased - aliasReference).norm() <= tol * numext::maxi(1.0, aliasReference.norm()));
+}
+
+// Repeated nodes make the matrix exactly singular, while non-finite nodes make
+// the solver input invalid. Both must be reported before solve evaluation.
+void test_bjorck_pereyra_singular() {
+  typedef Matrix<double, Dynamic, 1> Vec;
+  Vec x(5);
+  x << 0.1, 0.7, 0.3, 0.7, 0.9;
+  Vandermonde<double> V(x);
+  BjorckPereyra<double> bp(V);
+  VERIFY(bp.info() == NumericalIssue);
+
+  x[3] = std::numeric_limits<double>::quiet_NaN();
+  bp.compute(Vandermonde<double>(x));
+  VERIFY(bp.info() == InvalidInput);
+}
+
+// Aliased products: the products carry the default product tag, so assignment
+// materializes a temporary exactly like a dense product and x = V * x,
+// x += V * x must come out as if the right-hand side had been copied first.
+template <typename Scalar>
+void test_vandermonde_aliased_product(Index n) {
+  typedef Matrix<Scalar, Dynamic, 1> Vec;
+  typedef Matrix<Scalar, Dynamic, Dynamic> Mat;
+
+  Vec x = Vec::Random(n);
+  Vandermonde<Scalar> V(x);  // square, so x = V * x type-checks
+  Mat dense = reference_vandermonde<Scalar>(x, n);
+
+  Vec b = Vec::Random(n);
+  Vec y = b;
+  y = V * y;
+  VERIFY_IS_APPROX(y, (dense * b).eval());
+
+  y = b;
+  y += V * y;
+  VERIFY_IS_APPROX(y, (b + dense * b).eval());
+
+  y = b;
+  y -= V * y;
+  VERIFY_IS_APPROX(y, (b - dense * b).eval());
+
+  Mat B = Mat::Random(n, 3);
+  Mat Y = B;
+  Y = V * Y;
+  VERIFY_IS_APPROX(Y, (dense * B).eval());
+}
+
+// Aliasing beyond the same-object case: the default-product temporary must also
+// resolve right-hand-side expressions that reference the destination, overlapping
+// views of one buffer, and rectangular self-assignments where the destination is
+// resized by the assignment.
+template <typename Scalar>
+void test_vandermonde_aliased_expression(Index n) {
+  typedef Matrix<Scalar, Dynamic, 1> Vec;
+  typedef Matrix<Scalar, Dynamic, Dynamic> Mat;
+
+  Vec nodes = Vec::Random(n);
+  Vandermonde<Scalar> V(nodes);  // square
+  Mat dense = reference_vandermonde<Scalar>(nodes, n);
+
+  // Right-hand-side expression referencing the destination.
+  Vec x = Vec::Random(n), x0 = x;
+  x = V * (x + Vec::Ones(n));
+  VERIFY_IS_APPROX(x, (dense * (x0 + Vec::Ones(n))).eval());
+
+  // Overlapping (shifted) segments of one buffer.
+  Vec buf = Vec::Random(n + 1);
+  Vec expected = dense * buf.tail(n);
+  buf.head(n) = V * buf.tail(n);
+  VERIFY_IS_APPROX(buf.head(n).eval(), expected);
+
+  // Rectangular self-assignment (the reviewer's 5x4 case): x = V * x resizes
+  // the destination from 4 to 5, so the product must be captured before the
+  // destination storage is touched.
+  Vec tallNodes = Vec::Random(5);
+  Vandermonde<Scalar> Vt(tallNodes, 4);
+  Mat denseT = reference_vandermonde<Scalar>(tallNodes, 4);
+  Vec z = Vec::Random(4), z0 = z;
+  z = Vt * z;
+  VERIFY_IS_EQUAL(z.size(), 5);
+  VERIFY_IS_APPROX(z, (denseT * z0).eval());
+}
+
+// A delayed-evaluated product expression must keep its structured factor alive:
+// the makeVandermonde() factory returns an owning temporary, so Product has to
+// nest the operator by value (no NestByRefBit) or `expr` dangles.
+template <typename Scalar>
+void test_vandermonde_delayed_product(Index m, Index n) {
+  typedef Matrix<Scalar, Dynamic, 1> Vec;
+  typedef Matrix<Scalar, Dynamic, Dynamic> Mat;
+
+  STATIC_CHECK(!std::is_reference<typename internal::ref_selector<Vandermonde<Scalar>>::type>::value);
+
+  Vec x = Vec::Random(m), a = Vec::Random(n);
+  Mat dense = reference_vandermonde<Scalar>(x, n);
+
+  auto expr = makeVandermonde(x, n) * a;  // the factory temporary dies with the full expression
+  Vec scribble = Vec::Random(2 * m);      // reuses the temporary's freed heap storage
+  Vec y = expr;
+  VERIFY_IS_APPROX(y, (dense * a).eval());
+  VERIFY_IS_EQUAL(scribble.size(), 2 * m);  // keep the scribble alive across the evaluation
+}
+
+// Mixed-scalar products: a real operator applied to a complex right-hand side
+// (and a complex operator applied to a real one) promotes to the complex
+// product scalar, so alpha and the Horner accumulation must run in the promoted
+// type rather than the operator scalar.
+template <typename RealScalar>
+void test_vandermonde_mixed_scalar(Index m, Index n) {
+  typedef std::complex<RealScalar> Complex;
+  typedef Matrix<RealScalar, Dynamic, 1> RVec;
+  typedef Matrix<Complex, Dynamic, 1> CVec;
+  typedef Matrix<Complex, Dynamic, Dynamic> CMat;
+
+  RVec x = RVec::Random(m);
+  Vandermonde<RealScalar> V(x, n);
+  CMat dense = reference_vandermonde<RealScalar>(x, n).template cast<Complex>();
+
+  CVec a = CVec::Random(n);
+  CVec y = V * a;
+  VERIFY_IS_APPROX(y, (dense * a).eval());
+
+  CVec y0 = CVec::Random(m);
+  y = y0;
+  y.noalias() += V * a;
+  VERIFY_IS_APPROX(y, (y0 + dense * a).eval());
+
+  CVec xc = CVec::Random(m);
+  Vandermonde<Complex> Vc(xc, n);
+  CMat denseC = reference_vandermonde<Complex>(xc, n);
+  RVec ar = RVec::Random(n);
+  CVec z = Vc * ar;
+  VERIFY_IS_APPROX(z, (denseC * ar).eval());
+}
+
+template <typename Scalar>
+void test_vandermonde_determinant(Index n) {
+  typedef Matrix<Scalar, Dynamic, 1> Vec;
+  typedef Matrix<Scalar, Dynamic, Dynamic> Mat;
+
+  Vec x = Vec::Random(n);
+  Vandermonde<Scalar> V(x);
+  Mat dense = reference_vandermonde<Scalar>(x, n);
+  VERIFY_IS_APPROX(V.determinant(), dense.determinant());
+}
+
+// Reference for the wide-dynamic-range determinant tests: the same factor sequence with explicit
+// exponent tracking, halving any factor whose node difference overflows (exact, since a difference
+// only overflows for huge normal operands). Power-of-two rescaling is exact, so this reproduces the
+// exact product up to one rounding per multiplication, representable where naive partial products
+// leave the double range. Finite nodes only.
+double reference_vandermonde_det(const Matrix<double, Dynamic, 1>& x) {
+  double mantissa = 1.0;
+  Index exponent = 0;
+  for (Index j = 1; j < x.size(); ++j)
+    for (Index i = 0; i < j; ++i) {
+      double diff = x[j] - x[i];
+      if (!(numext::isfinite)(diff)) {
+        diff = 0.5 * x[j] - 0.5 * x[i];
+        ++exponent;
+      }
+      int e;
+      mantissa *= std::frexp(diff, &e);
+      exponent += e;
+      mantissa = std::frexp(mantissa, &e);
+      exponent += e;
+    }
+  return std::ldexp(mantissa, static_cast<int>(exponent));
+}
+
+// Balanced determinant accumulation (reviewer repro on MR 2691): the running
+// product is kept as mantissa * 2^e with exact frexp/ldexp renormalization, so
+// partial products that leave the representable range cannot destroy a
+// representable determinant, while zeros and genuine overflow still propagate.
+void test_vandermonde_determinant_scaled() {
+  typedef Matrix<double, Dynamic, 1> Vec;
+  typedef std::complex<double> Complex;
+  // The implementation accumulates the identical factor sequence as the
+  // reference above, so the two sides differ only by the exactness of the
+  // power-of-two renormalization.
+  const double kBalancedTol = 16 * NumTraits<double>::epsilon();
+  // The analytic values below additionally absorb the representation error of
+  // the decimal node literals and of pow(10, 27.5), amplified by the 15-factor
+  // product.
+  const double kAnalyticTol = 100 * NumTraits<double>::epsilon();
+
+  const double a = 1e-110, M = std::pow(10.0, 27.5);
+  {
+    // Underflow side (the reviewer's exact node set): the a-spaced factors alone
+    // reach 2e-330, below the smallest subnormal, so the naive running product
+    // flushes to an exact zero -- yet the determinant is 864 * a^3 * M^12 ~ 864.
+    Vec x(6);
+    x << 0.0, a, 2 * a, M, 2 * M, 3 * M;
+    Vandermonde<double> V(x);
+    const double det = V.determinant();
+    const double ref = reference_vandermonde_det(x);
+    VERIFY((numext::isfinite)(det));
+    VERIFY(numext::abs(det / ref - 1.0) <= kBalancedTol);
+    VERIFY(numext::abs(det / 864.0 - 1.0) <= kAnalyticTol);
+  }
+  {
+    // Overflow-side analogue: with the large nodes first the naive running
+    // product tops out at ~4.3e327 (infinity) three factors before the end, yet
+    // the determinant is a representable -864 * B^12 * s^3 = -8.64e305.
+    const double B = 1e28, s = 1e-11;
+    Vec x(6);
+    x << B, 2 * B, 3 * B, 0.0, s, 2 * s;
+    Vandermonde<double> V(x);
+    const double det = V.determinant();
+    const double ref = reference_vandermonde_det(x);
+    VERIFY((numext::isfinite)(det));
+    VERIFY(numext::abs(det / ref - 1.0) <= kBalancedTol);
+    VERIFY(numext::abs(det / -8.64e305 - 1.0) <= kAnalyticTol);
+  }
+  {
+    // Genuinely overflowing determinant, ~3.5e424: must saturate to +infinity.
+    const double T = 1e28;
+    Vec x(6);
+    x << 0.0, T, 2 * T, 3 * T, 4 * T, 5 * T;
+    Vandermonde<double> V(x);
+    const double det = V.determinant();
+    VERIFY((numext::isinf)(det) && det > 0.0);
+  }
+  {
+    // Genuinely underflowing determinant: subnormal-spaced nodes give 15 factors
+    // below 2^-1071, a product near 2^-16000, which must saturate to zero with
+    // the sign of the factored form: +0 for ascending nodes (all factors
+    // positive), -0 for the descending permutation (15 negative factors).
+    const double d = std::numeric_limits<double>::denorm_min();
+    Vec x(6);
+    x << 0.0, d, 2 * d, 3 * d, 4 * d, 5 * d;
+    const double det = Vandermonde<double>(x).determinant();
+    VERIFY(det == 0.0 && !std::signbit(det));
+    const double detNeg = Vandermonde<double>(Vec(x.reverse())).determinant();
+    VERIFY(detNeg == 0.0 && std::signbit(detNeg));
+  }
+  {
+    // A repeated node makes the matrix exactly singular: the zero factor must
+    // propagate to an exact 0 even though the naive running product overflows
+    // beforehand (which would end in 0 * inf = NaN).
+    Vec x(5);
+    x << 0.0, 1e300, 1e-300, 1e300, 5.0;
+    Vandermonde<double> V(x);
+    VERIFY(V.determinant() == 0.0);
+  }
+  {
+    // Complex path of the balanced accumulation: purely imaginary nodes i*x
+    // turn every factor of the underflow-side case into i*(x_j - x_i), so
+    // det = i^15 * 864 = -864i.
+    Vec xr(6);
+    xr << 0.0, a, 2 * a, M, 2 * M, 3 * M;
+    Matrix<Complex, Dynamic, 1> x = Complex(0, 1) * xr.cast<Complex>();
+    Vandermonde<Complex> V(x);
+    const Complex det = V.determinant();
+    const Complex ref(0.0, -reference_vandermonde_det(xr));
+    VERIFY((numext::isfinite)(numext::abs(det)));
+    VERIFY(numext::abs(det - ref) <= kBalancedTol * numext::abs(ref));
+  }
+}
+
+// Reviewer reproducer on MR 2691: a node difference can overflow while the determinant stays
+// representable. For nodes [-DBL_MAX, DBL_MAX, d, 2d, ..., 5d] with d the smallest subnormal,
+// DBL_MAX - (-DBL_MAX) = 2^1025 overflows, yet det = -576 DBL_MAX^11 d^10 ~ -3.1633e160. The
+// overflowing difference must be halved (exact for these huge normal operands) and the factor of
+// two carried in the running exponent.
+void test_vandermonde_determinant_overflowing_differences() {
+  typedef Matrix<double, Dynamic, 1> Vec;
+  typedef std::complex<double> Complex;
+  const double M = (std::numeric_limits<double>::max)();
+  const double d = std::numeric_limits<double>::denorm_min();
+  const double kBalancedTol = 16 * NumTraits<double>::epsilon();
+  const double kAnalyticTol = 100 * NumTraits<double>::epsilon();
+
+  Vec x(7);
+  x << -M, M, d, 2 * d, 3 * d, 4 * d, 5 * d;
+  const double det = Vandermonde<double>(x).determinant();
+  VERIFY((numext::isfinite)(det));
+
+  // Scaled reference computation with the identical factor sequence.
+  const double ref = reference_vandermonde_det(x);
+  VERIFY(numext::abs(det / ref - 1.0) <= kBalancedTol);
+
+  // Analytic factored form, accumulated with the same exponent bookkeeping:
+  // every huge factor rounds to +-DBL_MAX (and the halved leading factor is
+  // exactly DBL_MAX * 2), the d-spaced block contributes exactly 288 * 2^-10740,
+  // so det = -2 * 288 * DBL_MAX^11 * 2^-10740.
+  double mantissa = -576.0;
+  Index exponent = -10740;
+  for (int t = 0; t < 11; ++t) {
+    int e;
+    mantissa *= std::frexp(M, &e);
+    exponent += e;
+    mantissa = std::frexp(mantissa, &e);
+    exponent += e;
+  }
+  const double expected = std::ldexp(mantissa, static_cast<int>(exponent));
+  VERIFY(numext::abs(det / expected - 1.0) <= kAnalyticTol);
+  VERIFY(numext::abs(det / -3.1633e160 - 1.0) <= 1e-3);  // the reviewer's quoted value
+
+  // Complex nodes i*x: every factor becomes i*(x_j - x_i) -- the overflow now
+  // sits in the imaginary components, which the component-wise finiteness check
+  // must catch -- and det picks up i^21 = i.
+  Matrix<Complex, Dynamic, 1> xc = Complex(0, 1) * x.cast<Complex>();
+  const Complex detc = Vandermonde<Complex>(xc).determinant();
+  const Complex refc = Complex(0, 1) * Complex(ref);
+  VERIFY((numext::isfinite)(detc));
+  VERIFY(numext::abs(detc - refc) <= kBalancedTol * numext::abs(refc));
+}
+
+// Reviewer reproducer on MR 2691: Horner intermediates can overflow while the value stays
+// representable. At node 1/2 with coefficients [0, DBL_MAX, DBL_MAX] the intermediate 1.5 DBL_MAX
+// is Inf, but the value is exactly fl(0.75 DBL_MAX). The scaled path must return it bit-exactly,
+// the plain path must match the naive loop for moderate data, and genuinely unrepresentable values
+// must still saturate.
+void test_vandermonde_scaled_horner() {
+  typedef Matrix<double, Dynamic, 1> Vec;
+  typedef std::complex<double> Complex;
+  typedef Matrix<Complex, Dynamic, 1> CVec;
+  const double M = (std::numeric_limits<double>::max)();
+
+  {
+    // The reviewer's exact reproducer, in each accumulation form.
+    Vec x(1);
+    x << 0.5;
+    Vandermonde<double> V(x, 3);
+    Vec a(3);
+    a << 0.0, M, M;
+    Vec y = V * a;
+    VERIFY_IS_EQUAL(y[0], 0.75 * M);
+    y.setZero();
+    y.noalias() += V * a;
+    VERIFY_IS_EQUAL(y[0], 0.75 * M);
+    // Two identical columns exercise the per-column screening.
+    Matrix<double, Dynamic, Dynamic> A(3, 2);
+    A.col(0) = a;
+    A.col(1) = a;
+    Matrix<double, Dynamic, Dynamic> Y = V * A;
+    VERIFY_IS_EQUAL(Y(0, 0), 0.75 * M);
+    VERIFY_IS_EQUAL(Y(0, 1), 0.75 * M);
+  }
+  {
+    // Complex path: the same coefficients rotated by i keep the overflow in a
+    // single component; the value is i * fl(0.75 * DBL_MAX).
+    CVec xc(1);
+    xc << Complex(0.5, 0.0);
+    Vandermonde<Complex> Vc(xc, 3);
+    CVec ac(3);
+    ac << Complex(0), Complex(0, M), Complex(0, M);
+    CVec yc = Vc * ac;
+    VERIFY_IS_EQUAL(numext::real(yc[0]), 0.0);
+    VERIFY_IS_EQUAL(numext::imag(yc[0]), 0.75 * M);
+  }
+  {
+    // Genuine overflow must still saturate: value 3 * DBL_MAX at the node 2.
+    Vec x(1);
+    x << 2.0;
+    Vandermonde<double> V(x, 2);
+    Vec a(2);
+    a << M, M;
+    Vec y = V * a;
+    VERIFY((numext::isinf)(y[0]) && y[0] > 0.0);
+    a << M, -M;  // value -DBL_MAX: representable again, sign preserved
+    y = V * a;
+    VERIFY_IS_EQUAL(y[0], -M);
+  }
+  {
+    // A vanished accumulator has no scale (reviewer repro): after the running
+    // value becomes exactly zero -- a zero node annihilating it, or an exact
+    // cancellation -- the frame exponent must reset, or the trailing small
+    // coefficient underflows in the stale huge frame and the value 0 is
+    // returned instead of the coefficient itself.
+    const double mn = (std::numeric_limits<double>::min)();
+    const double dm = std::numeric_limits<double>::denorm_min();
+    Vec x(1);
+    x << 0.0;  // zero node: the value is exactly a[0]
+    Vandermonde<double> V0(x, 3);
+    Vec a(3);
+    a << mn, M, M;
+    Vec y = V0 * a;
+    VERIFY_IS_EQUAL(y[0], mn);
+
+    x << 1.0;  // exact cancellation: M - M + dm = dm
+    Vandermonde<double> V1(x, 3);
+    a << dm, -M, M;
+    y = V1 * a;
+    VERIFY_IS_EQUAL(y[0], dm);
+  }
+  {
+    // Complex variant: the real components cancel exactly while the tiny
+    // imaginary coefficient must survive in a fresh frame.
+    const double dm = std::numeric_limits<double>::denorm_min();
+    CVec xc(1);
+    xc << Complex(1.0, 0.0);
+    Vandermonde<Complex> Vc(xc, 3);
+    CVec ac(3);
+    ac << Complex(0.0, dm), Complex(-M, 0.0), Complex(M, 0.0);
+    CVec yc = Vc * ac;
+    VERIFY_IS_EQUAL(numext::real(yc[0]), 0.0);
+    VERIFY_IS_EQUAL(numext::imag(yc[0]), dm);
+  }
+  {
+    // Near-cancellation: M - (M - 2^971) leaves the tiny nonzero intermediate
+    // 2^971 (mantissa 2^-53 in the huge frame), which the frexp renormalization
+    // must rebase; the trailing coefficient 1 then rounds away exactly as in
+    // real arithmetic: fl(2^971 + 1) = 2^971.
+    const double big = std::ldexp(1.0, 971);  // ulp(DBL_MAX), and M - big is exact
+    Vec x(1);
+    x << 1.0;
+    Vandermonde<double> V(x, 3);
+    Vec a(3);
+    a << 1.0, -(M - big), M;
+    Vec y = V * a;
+    VERIFY_IS_EQUAL(y[0], big);
+  }
+  {
+    // Moderate data keeps the plain path, so this matches a naive Horner loop to Horner's own forward
+    // error bound (Higham, ASNA 2nd ed., section 5.1): |p(x) - fl(p(x))| <= gamma_{2n} sum_j |a_j| |x|^j
+    // <= 2n eps sum_j |a_j| for |x| <= 1. Not bit-identical, since the compiler may contract
+    // acc * x[i] + a[j] into an FMA in one of the two loops and not the other.
+    const Index m = 7, n = 6;
+    Vec x = Vec::Random(m), a = Vec::Random(n);
+    Vandermonde<double> V(x, n);
+    Vec y = V * a;
+    const double bound = double(2 * n) * NumTraits<double>::epsilon() * a.cwiseAbs().sum();
+    for (Index i = 0; i < m; ++i) {
+      double acc = a[n - 1];
+      for (Index j = n - 2; j >= 0; --j) acc = acc * x[i] + a[j];
+      VERIFY(numext::abs(y[i] - acc) <= bound);
+    }
+  }
+}
+
+template <typename Scalar, int M, int N>
+void test_vandermonde_fixed() {
+  typedef Matrix<Scalar, M, 1> NodeVec;
+  typedef Matrix<Scalar, Dynamic, 1> Vec;
+  typedef Matrix<Scalar, M, N> MatMN;
+
+  NodeVec x = NodeVec::Random();
+  Vandermonde<Scalar, M, N> V(x, N);
+  STATIC_CHECK((Vandermonde<Scalar, M, N>::RowsAtCompileTime == M));
+  STATIC_CHECK((Vandermonde<Scalar, M, N>::ColsAtCompileTime == N));
+  STATIC_CHECK((internal::remove_all_t<decltype(makeVandermonde(x))>::ColsAtCompileTime == M));
+
+  MatMN dense = V;
+  VERIFY_IS_APPROX(dense, MatMN(reference_vandermonde<Scalar>(Vec(x), N)));
+
+  Matrix<Scalar, N, 1> a = Matrix<Scalar, N, 1>::Random();
+  Matrix<Scalar, M, 1> y = V * a;
+  VERIFY_IS_APPROX(y, (dense * a).eval());
+}
+
+// Core rewrites alpha * (V * a) as (alpha * V) * a. Keep the two scaled
+// expression forms together because they share the same structured wrapper.
+void test_vandermonde_expression_regressions() {
+  typedef Matrix<double, 1, 1> Vec1;
+  typedef Matrix<double, 3, 1> Vec3;
+  typedef Matrix<double, 1, 3> Row3;
+  typedef Matrix<double, 3, 3> Mat3;
+  typedef Matrix<double, 3, 2> Mat32;
+
+  Vec1 x;
+  x << 0.5;
+  Vandermonde<double, 1, 3> V(x, 3);
+  Row3 dense = V;
+  Row3 scaledDense = 2.0 * V;
+  VERIFY_IS_APPROX(scaledDense, (2.0 * dense).eval());
+
+  Vec3 a;
+  a << 1.0, 2.0, 3.0;
+  Vec1 scaledProduct = 0.5 * (V * a);
+  VERIFY_IS_APPROX(scaledProduct, (0.5 * (dense * a)).eval());
+
+  // A product on the right must be evaluated before the structured Horner
+  // evaluator asks it for coefficients.
+  Mat3 B = Mat3::Random();
+  Mat32 C = Mat32::Random();
+  VERIFY_IS_APPROX((V * (B * C)).eval(), (dense * (B * C)).eval());
+}
+
+void test_vandermonde_dimension_checks() {
+  VectorXd x(2);
+  x << 0.0, 1.0;
+  VERIFY_RAISES_ASSERT((Vandermonde<double, Dynamic, 3>(x)));
+
+  Vandermonde<double> rectangular(x, 3);
+  VERIFY_RAISES_ASSERT((BjorckPereyra<double>(rectangular)));
+}
+
+// Fixed-size Vandermonde expressions own the operator by value so temporaries
+// returned from factories remain alive. Their node storage must therefore avoid
+// over-alignment: generic Product and cwise wrapper types do not supply aligned
+// operator new under C++14.
+void test_vandermonde_fixed_heap_expressions() {
+  typedef Matrix<double, 4, 1> Vec4;
+  typedef Matrix<double, 4, 4> Mat4;
+
+  Vec4 x;
+  x << -1.0, -0.25, 0.5, 1.0;
+  Vec4 a;
+  a << 1.0, 2.0, 3.0, 4.0;
+
+  typedef decltype(makeVandermonde(x) * a) ProductExpression;
+  typedef decltype(2.0 * makeVandermonde(x)) ScaledExpression;
+  STATIC_CHECK(alignof(ProductExpression) <= alignof(std::max_align_t));
+  STATIC_CHECK(alignof(ScaledExpression) <= alignof(std::max_align_t));
+
+  ProductExpression* delayedProduct = new ProductExpression(makeVandermonde(x) * a);
+  ScaledExpression* delayedScaled = new ScaledExpression(2.0 * makeVandermonde(x));
+
+  const Mat4 dense = reference_vandermonde<double>(Matrix<double, Dynamic, 1>(x), 4);
+  Vec4 y = *delayedProduct;
+  Mat4 scaled = *delayedScaled;
+  VERIFY_IS_APPROX(y, (dense * a).eval());
+  VERIFY_IS_APPROX(scaled, (2.0 * dense).eval());
+  delete delayedProduct;
+  delete delayedScaled;
+}
+
+EIGEN_DECLARE_TEST(structured_vandermonde) {
+  for (int i = 0; i < g_repeat; ++i) {
+    // Horner products, dense assignment, coefficient access.
+    CALL_SUBTEST_1((test_vandermonde_product<double>(1, 1)));
+    CALL_SUBTEST_1((test_vandermonde_product<double>(8, 8)));
+    CALL_SUBTEST_1((test_vandermonde_product<double>(20, 12)));  // tall
+    CALL_SUBTEST_1((test_vandermonde_product<double>(12, 20)));  // wide
+    CALL_SUBTEST_1((test_vandermonde_product<float>(10, 10)));
+    CALL_SUBTEST_1((test_vandermonde_product<std::complex<double>>(9, 7)));
+    CALL_SUBTEST_1((test_vandermonde_product<std::complex<float>>(7, 9)));
+
+    // Björck-Pereyra primal/dual/adjoint solves.
+    CALL_SUBTEST_2((test_bjorck_pereyra<double>(1)));
+    CALL_SUBTEST_2((test_bjorck_pereyra<double>(2)));
+    CALL_SUBTEST_2((test_bjorck_pereyra<double>(10)));
+    CALL_SUBTEST_2((test_bjorck_pereyra<double>(14)));
+    CALL_SUBTEST_2((test_bjorck_pereyra<std::complex<double>>(10)));
+    CALL_SUBTEST_2((test_bjorck_pereyra_roots_of_unity<double>(64)));
+    CALL_SUBTEST_2((test_bjorck_pereyra_roots_of_unity<float>(32)));
+    CALL_SUBTEST_2(test_vandermonde_half());
+    CALL_SUBTEST_2(test_bjorck_pereyra_higham());
+    CALL_SUBTEST_2(test_bjorck_pereyra_complex_transpose_adjoint(12));
+    CALL_SUBTEST_2(test_bjorck_pereyra_singular());
+
+    // Closed-form determinant and fixed sizes.
+    CALL_SUBTEST_3((test_vandermonde_determinant<double>(8)));
+    CALL_SUBTEST_3((test_vandermonde_determinant<std::complex<double>>(7)));
+    CALL_SUBTEST_3(test_vandermonde_determinant_scaled());
+    CALL_SUBTEST_3(test_vandermonde_determinant_overflowing_differences());
+    CALL_SUBTEST_3(test_vandermonde_scaled_horner());
+    CALL_SUBTEST_3((test_vandermonde_fixed<double, 6, 4>()));
+    CALL_SUBTEST_3((test_vandermonde_fixed<double, 5, 5>()));
+    CALL_SUBTEST_3((test_vandermonde_fixed<std::complex<float>, 4, 4>()));
+    CALL_SUBTEST_3(test_vandermonde_expression_regressions());
+    CALL_SUBTEST_3(test_vandermonde_fixed_heap_expressions());
+
+    // Product regressions: aliasing, operand lifetime, mixed scalars.
+    CALL_SUBTEST_4((test_vandermonde_aliased_product<double>(8)));
+    CALL_SUBTEST_4((test_vandermonde_aliased_product<double>(17)));
+    CALL_SUBTEST_4((test_vandermonde_aliased_product<std::complex<double>>(9)));
+    CALL_SUBTEST_4((test_vandermonde_aliased_expression<double>(8)));
+    CALL_SUBTEST_4((test_vandermonde_aliased_expression<std::complex<double>>(9)));
+    CALL_SUBTEST_4((test_vandermonde_delayed_product<double>(12, 9)));
+    CALL_SUBTEST_4((test_vandermonde_delayed_product<std::complex<double>>(8, 8)));
+    CALL_SUBTEST_4((test_vandermonde_mixed_scalar<double>(10, 8)));
+    CALL_SUBTEST_4((test_vandermonde_mixed_scalar<float>(7, 9)));
+  }
+
+  CALL_SUBTEST_3(test_vandermonde_dimension_checks());
+}
diff --git a/unsupported/test/structured_vandermonde_int_index.cpp b/unsupported/test/structured_vandermonde_int_index.cpp
new file mode 100644
index 0000000..7897cd6
--- /dev/null
+++ b/unsupported/test/structured_vandermonde_int_index.cpp
@@ -0,0 +1,51 @@
+// SPDX-FileCopyrightText: The Eigen Authors
+// SPDX-License-Identifier: MPL-2.0
+
+#ifdef EIGEN_DEFAULT_DENSE_INDEX_TYPE
+#undef EIGEN_DEFAULT_DENSE_INDEX_TYPE
+#endif
+#define EIGEN_DEFAULT_DENSE_INDEX_TYPE int
+#include "main.h"
+
+#include <unsupported/Eigen/StructuredMatrices>
+
+using namespace Eigen;
+
+// Roughly 2.2 million factors near the top of the double range exceed a
+// 32-bit Index exponent even though the determinant has the unambiguous result
+// +Inf. The internal exponent type must not wrap.
+void test_vandermonde_determinant_int_index() {
+  typedef Matrix<double, Dynamic, 1> Vec;
+  const Index n = 2100;
+  const double low = -9e307, high = 9e307;
+  Vec x(n);
+  for (Index i = 0; i < n; ++i) {
+    const double t = double(i) / double(n - 1);
+    x[i] = (1.0 - t) * low + t * high;
+  }
+
+  const double det = Vandermonde<double>(x).determinant();
+  VERIFY((numext::isinf)(det));
+  VERIFY(det > 0.0);
+}
+
+// The scaled-Horner frame can likewise exceed a 32-bit Index exponent. A
+// degree-2,099,999 monomial at DBL_MAX must overflow to +Inf rather than wrap
+// its accumulated exponent and underflow to zero.
+void test_vandermonde_scaled_horner_int_index() {
+  typedef Matrix<double, Dynamic, 1> Vec;
+  const Index n = 2100000;
+  Vec x(1), a = Vec::Zero(n);
+  x[0] = (std::numeric_limits<double>::max)();
+  a[n - 1] = 1.0;
+
+  Vec y = Vandermonde<double>(x, n) * a;
+  VERIFY((numext::isinf)(y[0]));
+  VERIFY(y[0] > 0.0);
+}
+
+EIGEN_DECLARE_TEST(structured_vandermonde_int_index) {
+  STATIC_CHECK(sizeof(Index) == sizeof(int));
+  CALL_SUBTEST_1(test_vandermonde_determinant_int_index());
+  CALL_SUBTEST_2(test_vandermonde_scaled_horner_int_index());
+}