StructuredMatrices: Add Vandermonde operator and Bjorck-Pereyra solver libeigen/eigen!2691 Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com> Co-authored-by: Rasmus Munk Larsen <rlarsen@nvidia.com>
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()); +}