// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2026 Rasmus Munk Larsen <rmlarsen@gmail.com>
//
// 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-License-Identifier: MPL-2.0

#ifndef EIGEN_TRIDIAGONAL_BISECTION_H
#define EIGEN_TRIDIAGONAL_BISECTION_H

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

namespace Eigen {

/** \eigenvalues_module \ingroup Eigenvalues_Module
 *
 * \brief Selects which eigenvalues to compute in a spectral-bisection solve.
 *
 * Passed to TridiagonalEigenSolver::computeEigenvalues(). The Sturm-sequence
 * bisection algorithm can compute an arbitrary contiguous subset of the
 * spectrum, which the implicit-QR algorithm cannot. Three modes are provided:
 *  - all(): every eigenvalue (the default),
 *  - indices(il, iu): the eigenvalues with 0-based ascending indices in the
 *    half-open range [il, iu) (i.e. the (iu - il) smallest starting at il),
 *  - values(vl, vu): the eigenvalues lying in the half-open interval [vl, vu)
 *    (lower-closed, upper-open, like indices(); an eigenvalue exactly equal to
 *    vl is selected, one equal to vu is not). Endpoint membership is resolved
 *    against the computed eigenvalue approximations, so an eigenvalue within
 *    roughly \c n*eps*||T|| of an endpoint -- but not exactly on it -- may
 *    resolve to either side of that endpoint.
 *
 * In every mode the computed eigenvalues are returned in non-decreasing order.
 */
struct EigenvalueRange {
  enum Kind { All, ByIndex, ByValue };

  Kind kind;
  Index il;
  Index iu;
  // Endpoints are held as long double so values() does not silently collapse two nearby endpoints
  // to the same number when the solver's RealScalar is long double (they are narrowed on use).
  long double vl;
  long double vu;

  /** \returns a selector for the entire spectrum. */
  static EigenvalueRange all() { return EigenvalueRange{All, 0, 0, 0.0L, 0.0L}; }

  /** \returns a selector for the eigenvalues with 0-based indices in [il, iu). */
  static EigenvalueRange indices(Index il, Index iu) { return EigenvalueRange{ByIndex, il, iu, 0.0L, 0.0L}; }

  /** \returns a selector for the eigenvalues in the half-open interval [vl, vu). */
  static EigenvalueRange values(long double vl, long double vu) { return EigenvalueRange{ByValue, 0, 0, vl, vu}; }
};

namespace internal {

/** \internal
 *
 * Evaluates Sturm sequences for a real symmetric tridiagonal matrix T at a
 * batch of shift points, using packet math to process many points in parallel.
 *
 * For a shift \c x the routine computes the number of eigenvalues of T that are
 * less than \c x via the recurrence
 *   \f$ q_0 = \alpha_0 - x,\quad q_i = (\alpha_i - \beta_{i-1}^2 / q_{i-1}) - x \f$,
 * counting how many pivots \f$ q_i \f$ are non-positive. To avoid division by
 * zero and overflow, every pivot that drops to or below \c pivmin is treated as
 * a (small) negative pivot, counted, and clamped to \c -pivmin. The same guard
 * is applied to every pivot including the first, so the scalar tail produces
 * counts identical to the vectorized body.
 *
 * \param[in] alpha     diagonal of T, length \c n.
 * \param[in] beta_sq   squared off-diagonal of T, \f$ \beta_i^2 \f$, length \c n-1.
 * \param[in] n         dimension of T.
 * \param[in] pivmin    smallest allowed pivot magnitude, strictly positive,
 *                      e.g. \c safemin * max_i beta_i^2 (floored to \c safemin).
 * \param[in] eval_points  the \c num_points shift values at which to evaluate.
 * \param[out] count    on output \c count[j] is the number of eigenvalues of T
 *                      less than \c eval_points[j], as an exact integer. Must
 *                      have room for \c num_points entries.
 */
// Evaluates the Sturm count for \c kUnroll packets of shift points at once, fully unrolled so the
// per-lane state (x, q, running count) stays in vector registers. The in-register count accumulates
// in RealScalar (one padd per row, no cross-lane traffic) and is flushed into the exact int64 output
// every kFlushPeriod rows -- comfortably before it could reach 2^digits, where RealScalar stops
// representing consecutive integers (float: 2^24) -- so the stored counts are exact for any n. The
// row loop is chunked on the flush period rather than testing per row, keeping the hot inner loop
// unchanged; matrices with n below the period (the usual case) run it exactly once.
// See tridiagonal_sturm_counts().
template <int kUnroll, typename RealScalar>
EIGEN_STRONG_INLINE void tridiagonal_sturm_block(const RealScalar* alpha, const RealScalar* beta_sq, Index n,
                                                 RealScalar pivmin, const RealScalar* eval_points,
                                                 numext::int64_t* count, Index start) {
  using Packet = typename packet_traits<RealScalar>::type;
  constexpr int kPacketSize = unpacket_traits<Packet>::size;
  constexpr int kFlushShift = NumTraits<RealScalar>::digits() - 2 < 30 ? NumTraits<RealScalar>::digits() - 2 : 30;
  constexpr Index kFlushPeriod = Index(1) << kFlushShift;
  const Packet pivmin_p = pset1<Packet>(pivmin);
  const Packet neg_pivmin_p = pset1<Packet>(-pivmin);
  const Packet one_p = pset1<Packet>(RealScalar(1));

  numext::int64_t* cnt = count + start;
  Packet x[kUnroll], q[kUnroll], c[kUnroll];
  const Packet alpha_0 = pset1<Packet>(alpha[0]);
  EIGEN_UNROLL_LOOP
  for (int k = 0; k < kUnroll; ++k) {
    x[k] = ploadu<Packet>(eval_points + start + k * kPacketSize);
    q[k] = psub(alpha_0, x[k]);
    const Packet mask = pcmp_le(q[k], pivmin_p);
    c[k] = pand(one_p, mask);
    q[k] = pselect(mask, pmin(q[k], neg_pivmin_p), q[k]);
  }
  // Convert the in-register lane counts to the exact int64 output: overwrite on the first flush,
  // add into the previously flushed partials thereafter. The scratch buffer lives outside the lambda:
  // MSVC captures kPacketSize instead of reading it as a constant, rejecting the bound inside.
  RealScalar buf[kPacketSize];
  auto store_counts = [&](bool accumulate) {
    EIGEN_UNROLL_LOOP
    for (int k = 0; k < kUnroll; ++k) {
      pstoreu(buf, c[k]);
      EIGEN_UNROLL_LOOP
      for (int l = 0; l < kPacketSize; ++l) {
        if (accumulate)
          cnt[k * kPacketSize + l] += numext::int64_t(buf[l]);
        else
          cnt[k * kPacketSize + l] = numext::int64_t(buf[l]);
      }
    }
  };
  bool flushed = false;
  Index i = 1;
  while (true) {
    const Index i_end = numext::mini(n, i + kFlushPeriod);
    for (; i < i_end; ++i) {
      const Packet alpha_i = pset1<Packet>(alpha[i]);
      const Packet beta_sq_im1 = pset1<Packet>(beta_sq[i - 1]);
      EIGEN_UNROLL_LOOP
      for (int k = 0; k < kUnroll; ++k) {
        // q = (alpha_i - beta_{i-1}^2 / q) - x.
        q[k] = psub(psub(alpha_i, pdiv(beta_sq_im1, q[k])), x[k]);
        const Packet mask = pcmp_le(q[k], pivmin_p);
        c[k] = padd(c[k], pand(one_p, mask));
        q[k] = pselect(mask, pmin(q[k], neg_pivmin_p), q[k]);
      }
    }
    if (i >= n) break;
    // Flush the in-register counts (exact: at most kFlushPeriod + 1 increments so far) into the
    // integer output and restart the accumulators. Only reached for n above the flush period.
    store_counts(flushed);
    EIGEN_UNROLL_LOOP
    for (int k = 0; k < kUnroll; ++k) c[k] = pzero(c[k]);
    flushed = true;
  }
  // Final conversion to the exact integer counts (adding any flushed partials).
  store_counts(flushed);
}

template <typename RealScalar>
void tridiagonal_sturm_counts(const RealScalar* alpha, const RealScalar* beta_sq, Index n, RealScalar pivmin,
                              const RealScalar* eval_points, numext::int64_t* count, Index num_points) {
  using Packet = typename packet_traits<RealScalar>::type;
  constexpr Index kPacketSize = Index(unpacket_traits<Packet>::size);
  // num_points is a point count and so never negative; clearing the sign bit makes that explicit to
  // the compiler, which lets the division by the power-of-two kPacketSize lower to a shift (and tidies
  // the trailing cascade loops) instead of emitting signed-division sign-correction code.
  num_points &= (std::numeric_limits<Index>::max)();
  const Index full = num_points / kPacketSize;  // number of whole packets

  // Process whole packets in register-resident blocks, cascading 8 -> 4 -> 2 -> 1 so that even a
  // small batch gets enough independent recurrences to hide the pdiv latency. Eight packets is the
  // measured sweet spot on AVX2 (more would spill the per-lane state out of vector registers).
  Index p = 0;
  for (; p + 8 <= full; p += 8)
    tridiagonal_sturm_block<8>(alpha, beta_sq, n, pivmin, eval_points, count, p * kPacketSize);
  for (; p + 4 <= full; p += 4)
    tridiagonal_sturm_block<4>(alpha, beta_sq, n, pivmin, eval_points, count, p * kPacketSize);
  for (; p + 2 <= full; p += 2)
    tridiagonal_sturm_block<2>(alpha, beta_sq, n, pivmin, eval_points, count, p * kPacketSize);
  for (; p + 1 <= full; p += 1)
    tridiagonal_sturm_block<1>(alpha, beta_sq, n, pivmin, eval_points, count, p * kPacketSize);

  // Scalar tail for the remaining (< kPacketSize) points. The pivot recurrence is bit-identical to
  // the packet path above; the count accumulates directly in the exact integer type.
  for (Index j = full * kPacketSize; j < num_points; ++j) {
    const RealScalar xj = eval_points[j];
    RealScalar qj = alpha[0] - xj;
    numext::int64_t cj = 0;
    if (qj <= pivmin) {
      ++cj;
      qj = numext::mini(qj, -pivmin);
    }
    for (Index i = 1; i < n; ++i) {
      qj = (alpha[i] - beta_sq[i - 1] / qj) - xj;
      if (qj <= pivmin) {
        ++cj;
        qj = numext::mini(qj, -pivmin);
      }
    }
    count[j] = cj;
  }
}

/** \internal
 *
 * Runs the (\c t_hi - \c t_lo) independent bisections for a contiguous block of
 * target eigenvalue indices in lock-step, writing the converged shift midpoints
 * (in the internally normalized scale) to \c out[0 .. t_hi-t_lo).
 *
 * This is the unit of parallel work in tridiagonal_bisection(): the bisections of
 * two disjoint index blocks share no state, so each thread owns one block and the
 * blocks communicate only at the join. The vectorized batch evaluator keeps every
 * lane of every active bracket busy within a block.
 *
 * \param[in] alpha       normalized diagonal of T (length \c n).
 * \param[in] beta_sq     normalized squared off-diagonal (length \c n-1).
 * \param[in] n           dimension of T.
 * \param[in] pivmin      smallest allowed pivot magnitude (see tridiagonal_sturm_counts()).
 * \param[in] bracket_lo  lower end of the initial search bracket (count == 0 below it).
 * \param[in] bracket_hi  upper end of the initial search bracket (count == n above it).
 * \param[in] t_lo, t_hi  0-based target indices for this block, half-open [t_lo, t_hi).
 * \param[in] max_iters   maximum number of bisection steps.
 * \param[in] abs_tol     absolute width below which a bracket is considered converged.
 * \param[out] out        converged midpoints for indices [t_lo, t_hi); length t_hi-t_lo.
 */
template <typename RealScalar>
void tridiagonal_bisection_block(const RealScalar* alpha, const RealScalar* beta_sq, Index n, RealScalar pivmin,
                                 RealScalar bracket_lo, RealScalar bracket_hi, Index t_lo, Index t_hi, int max_iters,
                                 RealScalar abs_tol, RealScalar* out) {
  using ArrayType = Array<RealScalar, Dynamic, 1>;
  using CountArrayType = Array<numext::int64_t, Dynamic, 1>;
  const Index m = t_hi - t_lo;
  if (m <= 0) return;

  // The eigenvalue with 0-based index i is the value where count(x) crosses from <= i to > i.
  // Counts and targets are exact 64-bit integers (see tridiagonal_sturm_counts()).
  ArrayType lower = ArrayType::Constant(m, bracket_lo);
  ArrayType upper = ArrayType::Constant(m, bracket_hi);
  const CountArrayType targets = CountArrayType::LinSpaced(m, t_lo, t_hi - 1);
  CountArrayType counts(m);
  ArrayType mid = RealScalar(0.5) * (lower + upper);

  // Each eigenvalue is recorded the iteration it first converges, using a per-element criterion that
  // depends only on that element's own bracket. This makes the result independent of how the index
  // range is grouped -- in particular bitwise identical for any number of threads -- because element
  // i's converged value never depends on when its neighbors finish. (A single global convergence test
  // would instead keep refining an already-converged eigenvalue until the slowest one in its group
  // caught up, so the last bits would shift with the thread count.)
  ArrayType result = mid;
  Array<bool, Dynamic, 1> done = Array<bool, Dynamic, 1>::Constant(m, false);

  // In the early iterations many eigenvalues still share a bracket, so their (sorted) midpoints are
  // identical and the Sturm count need only be evaluated at the distinct values and scattered back.
  // Once every midpoint is distinct the brackets only ever get finer, so we stop deduplicating and
  // evaluate all midpoints directly, avoiding the per-iteration bookkeeping (and its unpredictable
  // branch). Deduplication never changes the result: count(mid[i]) is unchanged. It only pays off
  // when there are many more points than packets (otherwise the per-point kernel cost is too small
  // to offset the bookkeeping), so the scratch is allocated only then.
  constexpr int kPacketSize = unpacket_traits<typename packet_traits<RealScalar>::type>::size;
  bool deduplicate = (m >= 64 * Index(kPacketSize));
  ArrayType distinct;
  CountArrayType counts_distinct;
  if (deduplicate) {
    distinct.resize(m);
    counts_distinct.resize(m);
  }
  for (int iter = 0; iter < max_iters; ++iter) {
    if (deduplicate) {
      const RealScalar* midp = mid.data();
      RealScalar* distp = distinct.data();
      Index nd = 0;
      distp[nd++] = midp[0];
      for (Index i = 1; i < m; ++i)
        if (midp[i] != midp[i - 1]) distp[nd++] = midp[i];
      tridiagonal_sturm_counts<RealScalar>(alpha, beta_sq, n, pivmin, distp, counts_distinct.data(), nd);
      const numext::int64_t* cdp = counts_distinct.data();
      numext::int64_t* cp = counts.data();
      Index g = 0;
      cp[0] = cdp[g];
      for (Index i = 1; i < m; ++i) {
        if (midp[i] != midp[i - 1]) ++g;
        cp[i] = cdp[g];
      }
      if (nd == m) deduplicate = false;  // all distinct: finer brackets stay distinct
    } else {
      tridiagonal_sturm_counts<RealScalar>(alpha, beta_sq, n, pivmin, mid.data(), counts.data(), m);
    }
    // count(mid) <= target  =>  eigenvalue is >= mid, raise the lower bound;
    // otherwise the eigenvalue is < mid, so lower the upper bound.
    const auto raise = (counts <= targets);
    lower = raise.select(mid, lower);
    upper = raise.select(upper, mid);
    const ArrayType new_mid = RealScalar(0.5) * (lower + upper);
    // Freeze each eigenvalue at the midpoint of its bracket the first time that bracket is tight
    // enough (width within abs_tol) or stops moving. Already-frozen entries keep their value.
    const auto converged = (new_mid == mid) || ((upper - lower) <= abs_tol);
    result = (!done && converged).select(new_mid, result);
    done = done || converged;
    mid = new_mid;
    if (done.all()) break;
  }
  // Any eigenvalue that did not converge within max_iters keeps its last midpoint, written straight
  // into the caller's output (no aliasing temporary: select is coefficient-wise and out is disjoint).
  Eigen::Map<ArrayType>(out, m) = done.select(result, mid);
}

/** \internal
 *
 * Returns the number of eigenvalues of the (normalized) symmetric tridiagonal matrix T that are
 * strictly less than \c x, evaluated by the same pivot recurrence as tridiagonal_sturm_counts() but
 * counting only strictly-negative pivots. This translates the endpoints of a half-open value range
 * [vl, vu) into a target index range. The count is only reliable when \c x is farther from every
 * eigenvalue than the recurrence's rounding displacement (~ \c n*eps*||T||): tridiagonal_bisection()
 * therefore evaluates it at endpoints widened outward by that bound -- yielding a guaranteed
 * superset of [vl, vu) -- and resolves exact endpoint membership afterwards, against the converged
 * eigenvalues themselves.
 *
 * As in the batch evaluator, a pivot whose magnitude reaches the floor \c pivmin is replaced by
 * +/- pivmin, keeping its sign, so the following division cannot overflow.
 *
 * \param[in] alpha    normalized diagonal of T (length \c n).
 * \param[in] beta_sq  normalized squared off-diagonal (length \c n-1).
 * \param[in] n        dimension of T.
 * \param[in] pivmin   smallest allowed pivot magnitude (see tridiagonal_sturm_counts()).
 * \param[in] x        the shift at which to count.
 */
template <typename RealScalar>
Index tridiagonal_sturm_count_below(const RealScalar* alpha, const RealScalar* beta_sq, Index n, RealScalar pivmin,
                                    RealScalar x) {
  RealScalar q = alpha[0] - x;
  Index count = (q < RealScalar(0)) ? 1 : 0;
  if (numext::abs(q) < pivmin) q = numext::copysign(pivmin, q);
  for (Index i = 1; i < n; ++i) {
    q = (alpha[i] - beta_sq[i - 1] / q) - x;
    if (q < RealScalar(0)) ++count;
    if (numext::abs(q) < pivmin) q = numext::copysign(pivmin, q);
  }
  return count;
}

/** \internal
 *
 * Computes a subset of the eigenvalues of a real symmetric tridiagonal matrix T
 * by Sturm-sequence spectral bisection (cf. LAPACK's xSTEBZ), accelerated by the
 * vectorized batch evaluator tridiagonal_sturm_counts().
 *
 * The eigenvalues are bracketed using the Gershgorin disc theorem, then refined
 * by running an independent bisection per requested eigenvalue, with all active
 * brackets evaluated together in each step so the packet evaluator stays full.
 *
 * \param[in] diag     diagonal of T (length \c n).
 * \param[in] subdiag  sub-diagonal of T (length \c n-1).
 * \param[in] range    which eigenvalues to compute (see EigenvalueRange).
 * \param[in] abs_tol  requested absolute accuracy; the effective tolerance is
 *                     max(abs_tol, eps * ||T||). Pass 0 for full precision.
 * \param[out] eivalues  filled with the selected eigenvalues in non-decreasing
 *                     order; resized to the number of selected eigenvalues.
 * \returns the number of eigenvalues written to \c eivalues.
 */
template <typename DiagType, typename SubdiagType, typename EivalType>
Index tridiagonal_bisection(const DiagType& diag, const SubdiagType& subdiag, const EigenvalueRange& range,
                            typename DiagType::Scalar abs_tol, EivalType& eivalues) {
  using RealScalar = typename DiagType::Scalar;
  using ArrayType = Array<RealScalar, Dynamic, 1>;
  EIGEN_STATIC_ASSERT(NumTraits<RealScalar>::IsInteger == 0 && NumTraits<RealScalar>::IsComplex == 0,
                      THIS_FUNCTION_IS_NOT_FOR_INTEGER_OR_COMPLEX_TYPES)

  const Index n = diag.size();
  if (n == 0) {
    eivalues.derived().resize(0);
    return 0;
  }

  // Normalize the matrix to O(1) to avoid overflow/underflow when squaring the
  // off-diagonal and during the Sturm recurrence; eigenvalues scale linearly, so the
  // scaling is undone at the very end. (The caller has already verified the input is
  // finite.) This mirrors the uniform scaling done in SelfAdjointEigenSolver::compute().
  // Divide each entry directly by the scale rather than multiplying by its reciprocal: when the scale
  // is subnormal, 1/scale overflows to infinity, the normalization silently disables itself, and the
  // Sturm recurrence underflows and returns wrong eigenvalues.
  RealScalar scale = diag.cwiseAbs().maxCoeff();
  if (n >= 2) scale = numext::maxi(scale, subdiag.cwiseAbs().maxCoeff());
  if (numext::is_exactly_zero(scale)) scale = RealScalar(1);

  // Local contiguous copies of the scaled matrix data, |off-diagonal|, and its square.
  const ArrayType alpha = diag.array() / scale;
  const ArrayType beta_abs = (n >= 2) ? ArrayType(subdiag.array().abs() / scale) : ArrayType(0);
  const ArrayType beta_sq = (n >= 2) ? ArrayType(beta_abs.square()) : ArrayType(0);

  // Smallest pivot allowed during the Sturm recurrence (positive, floored so
  // that a matrix with an all-zero off-diagonal still has pivmin > 0).
  const RealScalar eps = NumTraits<RealScalar>::epsilon();
  const RealScalar safemin = numext::maxi(RealScalar(1) / NumTraits<RealScalar>::highest(),
                                          (RealScalar(1) + eps) * (std::numeric_limits<RealScalar>::min)());
  const RealScalar max_beta_sq = (n >= 2) ? beta_sq.maxCoeff() : RealScalar(0);
  const RealScalar pivmin = safemin * numext::maxi(max_beta_sq, RealScalar(1));

  // Gershgorin bounds: row k has radius |beta_{k-1}| + |beta_k| (with the
  // missing boundary off-diagonals taken as zero).
  ArrayType radius(n);
  if (n == 1) {
    radius(0) = RealScalar(0);
  } else {
    radius(0) = beta_abs(0);
    radius(n - 1) = beta_abs(n - 2);
    if (n > 2) radius.segment(1, n - 2) = beta_abs.head(n - 2) + beta_abs.segment(1, n - 2);
  }
  RealScalar lambda_min = (alpha - radius).minCoeff();
  RealScalar lambda_max = (alpha + radius).maxCoeff();

  // Effective tolerance and outward expansion of the bracket so that
  // count(lambda_min) == 0 and count(lambda_max) == n (cf. LAPACK xSTEBZ).
  const RealScalar tnorm = numext::maxi(numext::abs(lambda_min), numext::abs(lambda_max));
  // Convergence tolerance: refine each bracket to ~1 ulp of ||T||. This is the role of xSTEBZ's
  // relative tolerance (it stops once b - a < RELFAC * ulp * max(|a|, |b|), RELFAC = 2), here
  // specialized to the matrix norm tnorm and folded together with any absolute tolerance the caller
  // requested.
  // The 2*pivmin floor mirrors xSTEBZ/xLAEBZ (which converge on max(abstol, pivmin, reltol*...)):
  // below the pivot floor the Sturm counts are meaningless anyway, and stopping there keeps the
  // bracket midpoints out of the subnormal range, where hardware with flush-to-zero packet
  // arithmetic (ARMv7 NEON) would evaluate them inconsistently between the packet and scalar paths.
  abs_tol = numext::maxi(abs_tol, numext::maxi(eps * tnorm, RealScalar(2) * pivmin));
  // Widen the Gershgorin bracket so that, despite rounding in the Sturm recurrence, count() really does
  // reach 0 at lambda_min and n at lambda_max. The n*eps*tnorm term bounds the worst-case count error
  // accumulated over the n recurrence steps; the 2*pivmin term covers the pivot floor. The 2.1 prefactor
  // is xSTEBZ's FUDGE factor: ideally 1 would suffice, but it is taken slightly larger to stay robust on
  // sloppy arithmetic, and (per xSTEBZ) widening the bracket this way only loosens the initial search
  // interval -- it has no effect on the accuracy of the converged eigenvalues.
  const RealScalar expand = RealScalar(2.1) * (RealScalar(n) * eps * tnorm + RealScalar(2) * pivmin);
  lambda_min -= expand;
  lambda_max += expand;

  // Determine the target indices [t_lo, t_hi) and the search bracket.
  Index t_lo = 0, t_hi = n;
  RealScalar bracket_lo = lambda_min, bracket_hi = lambda_max;
  bool value_filter = false;
  RealScalar vl_n = RealScalar(0), vu_n = RealScalar(0);
  RealScalar end_tol_lo = RealScalar(0), end_tol_hi = RealScalar(0);
  if (range.kind == EigenvalueRange::ByIndex) {
    eigen_assert(range.il >= 0 && range.il <= range.iu && range.iu <= n && "invalid eigenvalue index range");
    t_lo = range.il;
    t_hi = range.iu;
  } else if (range.kind == EigenvalueRange::ByValue) {
    eigen_assert(range.vl <= range.vu && "invalid eigenvalue value range");
    vl_n = RealScalar(range.vl) / scale;
    vu_n = RealScalar(range.vu) / scale;
    // A Sturm count taken exactly at an endpoint is unreliable when an eigenvalue sits there: the
    // rounded terminal pivot's sign is arbitrary, so an eigenvalue equal to vl could be silently
    // dropped, or one equal to vu kept, breaking the documented half-open [vl, vu) semantics.
    // Instead, count at endpoints widened outward by a bound on the count's rounding displacement
    // (the same form as the Gershgorin expansion above, per endpoint magnitude, plus the convergence
    // tolerance) -- a guaranteed superset of [vl, vu) -- converge that superset, and resolve
    // endpoint membership against the *converged* eigenvalues after the bisection (see below).
    // An infinite endpoint (given as such, or a finite long double that narrowed to infinity in
    // RealScalar) means unbounded on that side: the count is trivial there, and the tolerance
    // arithmetic must not run (inf - inf = NaN would empty the result).
    if ((numext::isfinite)(vl_n)) {
      end_tol_lo =
          RealScalar(2.1) * (RealScalar(n) * eps * numext::maxi(tnorm, numext::abs(vl_n)) + RealScalar(2) * pivmin) +
          abs_tol;
      t_lo = tridiagonal_sturm_count_below<RealScalar>(alpha.data(), beta_sq.data(), n, pivmin,
                                                       vl_n - RealScalar(2) * end_tol_lo);
      bracket_lo = numext::maxi(lambda_min, vl_n - RealScalar(2) * end_tol_lo);
    } else {
      t_lo = (vl_n < RealScalar(0)) ? 0 : n;
    }
    if ((numext::isfinite)(vu_n)) {
      end_tol_hi =
          RealScalar(2.1) * (RealScalar(n) * eps * numext::maxi(tnorm, numext::abs(vu_n)) + RealScalar(2) * pivmin) +
          abs_tol;
      t_hi = tridiagonal_sturm_count_below<RealScalar>(alpha.data(), beta_sq.data(), n, pivmin, vu_n + end_tol_hi);
      bracket_hi = numext::mini(lambda_max, vu_n + RealScalar(2) * end_tol_hi);
    } else {
      t_hi = (vu_n > RealScalar(0)) ? n : 0;
    }
    value_filter = true;
  }

  const Index m = t_hi - t_lo;
  eivalues.derived().resize(m);
  if (m <= 0) return 0;

  const int max_iters = NumTraits<RealScalar>::digits() + 2;
  ArrayType mid_all(m);
  RealScalar* out = mid_all.data();
  const RealScalar* alpha_p = alpha.data();
  const RealScalar* beta_sq_p = beta_sq.data();

  // The m eigenvalues are bisected independently, so the spectrum splits cleanly across threads with
  // a single fork/join: thread t owns the contiguous index block [t_lo + lo, t_lo + hi) and writes
  // its converged midpoints into the disjoint slice out[lo, hi). No communication until the join.
#if defined(EIGEN_HAS_OPENMP)
  int nthreads = 1;
  // Don't nest inside an existing parallel region, and only fork when there is enough work to
  // amortize the thread overhead while still leaving each thread a SIMD-friendly chunk of points.
  if (omp_get_num_threads() == 1) {
    constexpr Index kPacketSize = Index(unpacket_traits<typename packet_traits<RealScalar>::type>::size);
    // One work unit ~ one Sturm step (a packet division); kMinTaskSize is the minimum per thread.
    // Multiply in double: the product overflows a 32-bit Index for matrices well within reach.
    const double work = double(m) * double(n) * double(max_iters);
    const double kMinTaskSize = 131072.0;
    const Index work_threads = Index(work / kMinTaskSize);
    const Index point_threads = m / (8 * kPacketSize);
    const Index pb = numext::maxi(Index(1), numext::mini(work_threads, point_threads));
    nthreads = int(numext::mini(pb, Index(Eigen::nbThreads())));
  }
  if (nthreads > 1) {
#pragma omp parallel num_threads(nthreads)
    {
      const Index nt = omp_get_num_threads();
      const Index tid = omp_get_thread_num();
      // Balanced split: every thread gets floor(m/nt) or ceil(m/nt) consecutive indices, none empty.
      const Index lo = tid * m / nt;
      const Index hi = (tid + 1) * m / nt;
      tridiagonal_bisection_block<RealScalar>(alpha_p, beta_sq_p, n, pivmin, bracket_lo, bracket_hi, t_lo + lo,
                                              t_lo + hi, max_iters, abs_tol, out + lo);
    }
  } else
#endif
  {
    tridiagonal_bisection_block<RealScalar>(alpha_p, beta_sq_p, n, pivmin, bracket_lo, bracket_hi, t_lo, t_hi,
                                            max_iters, abs_tol, out);
  }

  // Enforce the documented non-decreasing order. Under uniform IEEE arithmetic the independent
  // brackets already yield sorted midpoints, but on hardware with non-uniform subnormal handling
  // (ARMv7, where NEON packet lanes flush to zero while the scalar VFP tail does not) neighbouring
  // targets refined through the two paths can disagree by ~pivmin at the subnormal boundary. The
  // cumulative max restores the invariant at O(m) cost and is a no-op for already-sorted output.
  for (Index j = 1; j < m; ++j) out[j] = numext::maxi(out[j], out[j - 1]);

  // ByValue: resolve endpoint membership against the converged eigenvalues. Comparing against both
  // endpoints shifted DOWN by the endpoint tolerance gives every eigenvalue within that tolerance of
  // an endpoint the documented closed/open treatment of the endpoint itself: an eigenvalue equal to
  // vl is deterministically kept and one equal to vu deterministically dropped, however the rounded
  // Sturm counts at the endpoints came out. Eigenvalues near (but not on) an endpoint may resolve to
  // either side, as with LAPACK xSTEBZ, whose endpoint counts carry the same rounding ambiguity.
  Index m_out = m;
  if (value_filter) {
    // Infinite endpoints pass through unshifted: every finite eigenvalue satisfies >= -inf / < +inf.
    const RealScalar keep_lo = (numext::isfinite)(vl_n) ? vl_n - end_tol_lo : vl_n;
    const RealScalar keep_hi = (numext::isfinite)(vu_n) ? vu_n - end_tol_hi : vu_n;
    Index k = 0;
    for (Index j = 0; j < m; ++j)
      if (out[j] >= keep_lo && out[j] < keep_hi) out[k++] = out[j];
    m_out = k;
    if (m_out != m) eivalues.derived().resize(m_out);
  }

  // Undo the normalization.
  eivalues = (mid_all.head(m_out) * scale).matrix();
  return m_out;
}

}  // namespace internal
}  // namespace Eigen

#endif  // EIGEN_TRIDIAGONAL_BISECTION_H