commit 7d7576f3262fa15c34d5575637bd8d7ff4a83f16
author Rasmus Munk Larsen <rmlarsen@google.com> Tue Jun 06 14:06:19 2023 -0700
committer Rasmus Munk Larsen <rmlarsen@google.com> Tue Jun 06 14:06:19 2023 -0700
tree a934eca40dba75f87d58e12836ee9955a6cef9fd
parent b7151ffaab545c812e047686d979f28ce5fef9c7

Avoid underflow in prsqrt.

Eigen/src/Core/MathFunctionsImpl.h

diff --git a/Eigen/src/Core/MathFunctionsImpl.h b/Eigen/src/Core/MathFunctionsImpl.h
index 5782db4..86c06b2 100644
--- a/Eigen/src/Core/MathFunctionsImpl.h
+++ b/Eigen/src/Core/MathFunctionsImpl.h
@@ -87,11 +87,11 @@
     Packet inv_sqrt = approx_rsqrt;
     for (int step = 0; step < Steps; ++step) {
       // Refine the approximation using one Newton-Raphson step:
-      // h_n = x * (inv_sqrt * inv_sqrt) - 1 (so that h_n is nearly 0).
+      // h_n = (x * inv_sqrt) * inv_sqrt - 1 (so that h_n is nearly 0).
       // inv_sqrt = inv_sqrt - 0.5 * inv_sqrt * h_n
-      Packet r2 = pmul(inv_sqrt, inv_sqrt);
+      Packet r2 = pmul(a, inv_sqrt);
       Packet half_r = pmul(inv_sqrt, cst_minus_half);
-      Packet h_n = pmadd(a, r2, cst_minus_one);
+      Packet h_n = pmadd(r2, inv_sqrt, cst_minus_one);
       inv_sqrt = pmadd(half_r, h_n, inv_sqrt);
     }