Fix qr_solve() failure on well-conditioned matrices with zero pivot (#1083)

* Fix qr_solve() failure on well-conditioned matrices with zero pivot

In householder(), the sign convention
    p[j] = -sign(Re(A[j,j])) * sqrt(s)
collapses to zero when A[j,j] is exactly zero, because ctx.sign(0) == 0.
A zero p[j] makes kappa = 1/(s - p[j]*A[j,j]) = 1/s instead of the
correct 2/||v||**2, corrupting the Householder reflection. The damage
propagates to subsequent columns and eventually trips the
"matrix is numerically singular" guard.

Default sign to ctx.one when A[j,j] is zero (matching LAPACK's dlarfg
convention) so the reflection is computed correctly. The existing
singularity check is left in place to catch genuinely zero column
slices.

Fixes #983.

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>

mpmath/matrices/linalg.py

@@ -371,7 +371,10 @@ class LinearAlgebraMethods:
             s = ctx.fsum(abs(A[i,j])**2 for i in range(j, m))
             if not abs(s) > ctx.eps:
                 raise ValueError('matrix is numerically singular')
-            p.append(-ctx.sign(ctx.re(A[j,j])) * ctx.sqrt(s))
+            sign = ctx.sign(ctx.re(A[j,j]))
+            if sign == 0:
+                sign = ctx.one
+            p.append(-sign * ctx.sqrt(s))
             kappa = ctx.one / (s - p[j] * A[j,j])
             A[j,j] -= p[j]
             for k in range(j+1, n):

mpmath/tests/test_linalg.py

@@ -201,6 +201,18 @@ def test_solve_overdet_complex():
     b = matrix([1 + j, 2, -j])
     assert norm(residual(A, lu_solve(A, b), b)) < 1.0208
 
+def test_qr_solve_issue_983():
+    A = matrix([[1, -pi/20, (-pi/20)**2, (-pi/20)**3],
+                [1, 0, 0, 0],
+                [1, pi / 20, (pi/20)**2, (pi/20)**3],
+                [1, pi/10, (pi/10)**2, (pi/10)**3]])
+    b = matrix([[mp.sin(-pi/20)],
+                [0],
+                [mp.sin(pi/20)],
+                [mp.sin(pi/20)]])
+    x, _ = qr_solve(A, b)
+    assert norm(residual(A, x, b), inf) < 1e-14
+
 def test_singular():
     A = [[5.6, 1.2], [7./15, .1]]
     B = repr(zeros(2))