diff --git a/benchmarks/LU/CMakeLists.txt b/benchmarks/LU/CMakeLists.txt
index e28693382..86483e664 100644
--- a/benchmarks/LU/CMakeLists.txt
+++ b/benchmarks/LU/CMakeLists.txt
@@ -1 +1,2 @@
 eigen_add_benchmark(bench_lu bench_lu.cpp)
+eigen_add_benchmark(bench_rcond bench_rcond.cpp)
diff --git a/benchmarks/LU/bench_rcond.cpp b/benchmarks/LU/bench_rcond.cpp
new file mode 100644
index 000000000..e1f5ea376
--- /dev/null
+++ b/benchmarks/LU/bench_rcond.cpp
@@ -0,0 +1,79 @@
+// 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/.
+
+// Benchmarks for reciprocal condition number estimation.
+//
+// Times the rcond() call on pre-computed factorizations.
+// The rcond estimator is O(n^2), so this isolates its cost from the O(n^3) factorization.
+
+#include <benchmark/benchmark.h>
+#include <Eigen/Dense>
+
+using namespace Eigen;
+
+// --- PartialPivLU ---
+
+template <typename Scalar>
+static void BM_PartialPivLU_Rcond(benchmark::State& state) {
+  const Index n = state.range(0);
+  using Mat = Matrix<Scalar, Dynamic, Dynamic>;
+  Mat A = Mat::Random(n, n);
+  A.diagonal().array() += Scalar(2 * n);
+  PartialPivLU<Mat> lu(A);
+  for (auto _ : state) {
+    auto rc = lu.rcond();
+    benchmark::DoNotOptimize(rc);
+  }
+  state.SetItemsProcessed(state.iterations());
+}
+
+// --- FullPivLU ---
+
+template <typename Scalar>
+static void BM_FullPivLU_Rcond(benchmark::State& state) {
+  const Index n = state.range(0);
+  using Mat = Matrix<Scalar, Dynamic, Dynamic>;
+  Mat A = Mat::Random(n, n);
+  A.diagonal().array() += Scalar(2 * n);
+  FullPivLU<Mat> lu(A);
+  for (auto _ : state) {
+    auto rc = lu.rcond();
+    benchmark::DoNotOptimize(rc);
+  }
+  state.SetItemsProcessed(state.iterations());
+}
+
+// --- LLT ---
+
+template <typename Scalar>
+static void BM_LLT_Rcond(benchmark::State& state) {
+  const Index n = state.range(0);
+  using Mat = Matrix<Scalar, Dynamic, Dynamic>;
+  Mat A = Mat::Random(n, n);
+  Mat SPD = A.adjoint() * A + Mat::Identity(n, n);
+  LLT<Mat> llt(SPD);
+  for (auto _ : state) {
+    auto rc = llt.rcond();
+    benchmark::DoNotOptimize(rc);
+  }
+  state.SetItemsProcessed(state.iterations());
+}
+
+// --- Size configurations ---
+
+static void SquareSizes(::benchmark::Benchmark* b) {
+  for (int n : {8, 32, 64, 128, 256, 512, 1024}) b->Arg(n);
+}
+
+BENCHMARK(BM_PartialPivLU_Rcond<float>)->Apply(SquareSizes)->Name("PartialPivLU_Rcond_float");
+BENCHMARK(BM_PartialPivLU_Rcond<double>)->Apply(SquareSizes)->Name("PartialPivLU_Rcond_double");
+BENCHMARK(BM_FullPivLU_Rcond<float>)->Apply(SquareSizes)->Name("FullPivLU_Rcond_float");
+BENCHMARK(BM_FullPivLU_Rcond<double>)->Apply(SquareSizes)->Name("FullPivLU_Rcond_double");
+BENCHMARK(BM_LLT_Rcond<float>)->Apply(SquareSizes)->Name("LLT_Rcond_float");
+BENCHMARK(BM_LLT_Rcond<double>)->Apply(SquareSizes)->Name("LLT_Rcond_double");
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index a70c1fc4c..285977443 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -256,6 +256,7 @@ ei_add_test(stable_norm)
 ei_add_test(permutationmatrices)
 ei_add_test(bandmatrix)
 ei_add_test(cholesky)
+ei_add_test(condition_estimator)
 ei_add_test(lu)
 ei_add_test(determinant)
 ei_add_test(inverse)
diff --git a/test/condition_estimator.cpp b/test/condition_estimator.cpp
new file mode 100644
index 000000000..94e9e1917
--- /dev/null
+++ b/test/condition_estimator.cpp
@@ -0,0 +1,265 @@
+// 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/.
+
+#include "main.h"
+#include <Eigen/Dense>
+
+template <typename MatrixType>
+typename MatrixType::RealScalar matrix_l1_norm(const MatrixType& m) {
+  return m.cwiseAbs().colwise().sum().maxCoeff();
+}
+
+template <typename MatrixType>
+void rcond_partial_piv_lu() {
+  typedef typename MatrixType::RealScalar RealScalar;
+  Index size = MatrixType::RowsAtCompileTime;
+  if (size == Dynamic) size = internal::random<Index>(2, EIGEN_TEST_MAX_SIZE);
+
+  // Create a random diagonally dominant (thus invertible) matrix.
+  MatrixType m = MatrixType::Random(size, size);
+  m.diagonal().array() += RealScalar(2 * size);
+
+  PartialPivLU<MatrixType> lu(m);
+  MatrixType m_inverse = lu.inverse();
+  RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m)) / matrix_l1_norm(m_inverse);
+  RealScalar rcond_est = lu.rcond();
+  // Verify the estimate is within a factor of 10 of the truth.
+  VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10);
+}
+
+template <typename MatrixType>
+void rcond_full_piv_lu() {
+  typedef typename MatrixType::RealScalar RealScalar;
+  Index size = MatrixType::RowsAtCompileTime;
+  if (size == Dynamic) size = internal::random<Index>(2, EIGEN_TEST_MAX_SIZE);
+
+  // Create a random diagonally dominant (thus invertible) matrix.
+  MatrixType m = MatrixType::Random(size, size);
+  m.diagonal().array() += RealScalar(2 * size);
+
+  FullPivLU<MatrixType> lu(m);
+  MatrixType m_inverse = lu.inverse();
+  RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m)) / matrix_l1_norm(m_inverse);
+  RealScalar rcond_est = lu.rcond();
+  VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10);
+}
+
+template <typename MatrixType>
+void rcond_llt() {
+  typedef typename MatrixType::RealScalar RealScalar;
+  Index size = MatrixType::RowsAtCompileTime;
+  if (size == Dynamic) size = internal::random<Index>(2, EIGEN_TEST_MAX_SIZE);
+
+  // Create a random SPD matrix: A^T * A + I.
+  MatrixType a = MatrixType::Random(size, size);
+  MatrixType m = a.adjoint() * a + MatrixType::Identity(size, size);
+
+  LLT<MatrixType> llt(m);
+  VERIFY(llt.info() == Success);
+  MatrixType m_inverse = llt.solve(MatrixType::Identity(size, size));
+  RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m)) / matrix_l1_norm(m_inverse);
+  RealScalar rcond_est = llt.rcond();
+  VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10);
+}
+
+template <typename MatrixType>
+void rcond_ldlt() {
+  typedef typename MatrixType::RealScalar RealScalar;
+  Index size = MatrixType::RowsAtCompileTime;
+  if (size == Dynamic) size = internal::random<Index>(2, EIGEN_TEST_MAX_SIZE);
+
+  // Create a random SPD matrix: A^T * A + I.
+  MatrixType a = MatrixType::Random(size, size);
+  MatrixType m = a.adjoint() * a + MatrixType::Identity(size, size);
+
+  LDLT<MatrixType> ldlt(m);
+  VERIFY(ldlt.info() == Success);
+  MatrixType m_inverse = ldlt.solve(MatrixType::Identity(size, size));
+  RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m)) / matrix_l1_norm(m_inverse);
+  RealScalar rcond_est = ldlt.rcond();
+  VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10);
+}
+
+template <typename MatrixType>
+void rcond_singular() {
+  typedef typename MatrixType::Scalar Scalar;
+  Index size = MatrixType::RowsAtCompileTime;
+  if (size == Dynamic) size = internal::random<Index>(2, EIGEN_TEST_MAX_SIZE);
+
+  // Create a rank-deficient matrix: first row is zero.
+  MatrixType m = MatrixType::Random(size, size);
+  m.row(0).setZero();
+
+  FullPivLU<MatrixType> lu(m);
+  VERIFY_IS_EQUAL(lu.rcond(), Scalar(0));
+}
+
+template <typename MatrixType>
+void rcond_identity() {
+  typedef typename MatrixType::RealScalar RealScalar;
+  Index size = MatrixType::RowsAtCompileTime;
+  if (size == Dynamic) size = internal::random<Index>(2, EIGEN_TEST_MAX_SIZE);
+
+  MatrixType m = MatrixType::Identity(size, size);
+
+  // All decompositions should give rcond ~= 1 for the identity.
+  {
+    PartialPivLU<MatrixType> lu(m);
+    VERIFY(lu.rcond() > RealScalar(0.5));
+  }
+  {
+    FullPivLU<MatrixType> lu(m);
+    VERIFY(lu.rcond() > RealScalar(0.5));
+  }
+  {
+    LLT<MatrixType> llt(m);
+    VERIFY(llt.rcond() > RealScalar(0.5));
+  }
+  {
+    LDLT<MatrixType> ldlt(m);
+    VERIFY(ldlt.rcond() > RealScalar(0.5));
+  }
+}
+
+template <typename MatrixType>
+void rcond_ill_conditioned() {
+  typedef typename MatrixType::RealScalar RealScalar;
+  Index size = MatrixType::RowsAtCompileTime;
+  if (size == Dynamic) size = internal::random<Index>(4, EIGEN_TEST_MAX_SIZE);
+
+  // Create a diagonal matrix with known large condition number.
+  // Use 1e-3 to stay well within single-precision range.
+  MatrixType m = MatrixType::Zero(size, size);
+  m(0, 0) = RealScalar(1);
+  for (Index i = 1; i < size; ++i) {
+    m(i, i) = RealScalar(1e-3);
+  }
+  // True condition number = 1e3, so rcond = 1e-3.
+
+  {
+    PartialPivLU<MatrixType> lu(m);
+    RealScalar rcond_est = lu.rcond();
+    VERIFY(rcond_est < RealScalar(1e-1));
+    VERIFY(rcond_est > RealScalar(1e-5));
+  }
+  {
+    FullPivLU<MatrixType> lu(m);
+    RealScalar rcond_est = lu.rcond();
+    VERIFY(rcond_est < RealScalar(1e-1));
+    VERIFY(rcond_est > RealScalar(1e-5));
+  }
+}
+
+template <typename MatrixType>
+void rcond_1x1() {
+  typedef typename MatrixType::RealScalar RealScalar;
+  typedef Matrix<typename MatrixType::Scalar, 1, 1> Mat1;
+  Mat1 m;
+  m(0, 0) = internal::random<RealScalar>(RealScalar(1), RealScalar(100));
+
+  {
+    PartialPivLU<Mat1> lu(m);
+    VERIFY_IS_APPROX(lu.rcond(), RealScalar(1));
+  }
+  {
+    FullPivLU<Mat1> lu(m);
+    VERIFY_IS_APPROX(lu.rcond(), RealScalar(1));
+  }
+  {
+    LLT<Mat1> llt(m);
+    VERIFY_IS_APPROX(llt.rcond(), RealScalar(1));
+  }
+  {
+    LDLT<Mat1> ldlt(m);
+    VERIFY_IS_APPROX(ldlt.rcond(), RealScalar(1));
+  }
+}
+
+template <typename MatrixType>
+void rcond_2x2() {
+  typedef typename MatrixType::RealScalar RealScalar;
+  typedef Matrix<typename MatrixType::Scalar, 2, 2> Mat2;
+
+  // Well-conditioned 2x2 matrix.
+  Mat2 m;
+  m << RealScalar(2), RealScalar(1), RealScalar(1), RealScalar(3);
+
+  {
+    PartialPivLU<Mat2> lu(m);
+    Mat2 m_inverse = lu.inverse();
+    RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m)) / matrix_l1_norm(m_inverse);
+    RealScalar rcond_est = lu.rcond();
+    VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10);
+  }
+  {
+    FullPivLU<Mat2> lu(m);
+    Mat2 m_inverse = lu.inverse();
+    RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m)) / matrix_l1_norm(m_inverse);
+    RealScalar rcond_est = lu.rcond();
+    VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10);
+  }
+  {
+    LLT<Mat2> llt(m);
+    Mat2 m_inverse = llt.solve(Mat2::Identity());
+    RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m)) / matrix_l1_norm(m_inverse);
+    RealScalar rcond_est = llt.rcond();
+    VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10);
+  }
+}
+
+EIGEN_DECLARE_TEST(condition_estimator) {
+  for (int i = 0; i < g_repeat; i++) {
+    CALL_SUBTEST_1(rcond_partial_piv_lu<Matrix3f>());
+    CALL_SUBTEST_1(rcond_full_piv_lu<Matrix3f>());
+    CALL_SUBTEST_1(rcond_llt<Matrix3f>());
+    CALL_SUBTEST_1(rcond_ldlt<Matrix3f>());
+    CALL_SUBTEST_1(rcond_singular<Matrix3f>());
+    CALL_SUBTEST_1(rcond_identity<Matrix3f>());
+    CALL_SUBTEST_1(rcond_1x1<Matrix3f>());
+    CALL_SUBTEST_1(rcond_2x2<Matrix3f>());
+
+    CALL_SUBTEST_2(rcond_partial_piv_lu<Matrix4d>());
+    CALL_SUBTEST_2(rcond_full_piv_lu<Matrix4d>());
+    CALL_SUBTEST_2(rcond_llt<Matrix4d>());
+    CALL_SUBTEST_2(rcond_ldlt<Matrix4d>());
+    CALL_SUBTEST_2(rcond_singular<Matrix4d>());
+    CALL_SUBTEST_2(rcond_identity<Matrix4d>());
+    CALL_SUBTEST_2(rcond_2x2<Matrix4d>());
+
+    CALL_SUBTEST_3(rcond_partial_piv_lu<MatrixXf>());
+    CALL_SUBTEST_3(rcond_full_piv_lu<MatrixXf>());
+    CALL_SUBTEST_3(rcond_llt<MatrixXf>());
+    CALL_SUBTEST_3(rcond_ldlt<MatrixXf>());
+    CALL_SUBTEST_3(rcond_singular<MatrixXf>());
+    CALL_SUBTEST_3(rcond_identity<MatrixXf>());
+    CALL_SUBTEST_3(rcond_ill_conditioned<MatrixXf>());
+
+    CALL_SUBTEST_4(rcond_partial_piv_lu<MatrixXd>());
+    CALL_SUBTEST_4(rcond_full_piv_lu<MatrixXd>());
+    CALL_SUBTEST_4(rcond_llt<MatrixXd>());
+    CALL_SUBTEST_4(rcond_ldlt<MatrixXd>());
+    CALL_SUBTEST_4(rcond_singular<MatrixXd>());
+    CALL_SUBTEST_4(rcond_identity<MatrixXd>());
+    CALL_SUBTEST_4(rcond_ill_conditioned<MatrixXd>());
+
+    CALL_SUBTEST_5(rcond_partial_piv_lu<MatrixXcf>());
+    CALL_SUBTEST_5(rcond_full_piv_lu<MatrixXcf>());
+    CALL_SUBTEST_5(rcond_llt<MatrixXcf>());
+    CALL_SUBTEST_5(rcond_ldlt<MatrixXcf>());
+    CALL_SUBTEST_5(rcond_singular<MatrixXcf>());
+    CALL_SUBTEST_5(rcond_identity<MatrixXcf>());
+
+    CALL_SUBTEST_6(rcond_partial_piv_lu<MatrixXcd>());
+    CALL_SUBTEST_6(rcond_full_piv_lu<MatrixXcd>());
+    CALL_SUBTEST_6(rcond_llt<MatrixXcd>());
+    CALL_SUBTEST_6(rcond_ldlt<MatrixXcd>());
+    CALL_SUBTEST_6(rcond_singular<MatrixXcd>());
+    CALL_SUBTEST_6(rcond_identity<MatrixXcd>());
+  }
+}