diff --git a/mpmath/matrices/linalg.py b/mpmath/matrices/linalg.py
index 88703df7..2dad1fde 100644
--- a/mpmath/matrices/linalg.py
+++ b/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):
diff --git a/mpmath/tests/test_linalg.py b/mpmath/tests/test_linalg.py
index 6b929110..58f29069 100644
--- a/mpmath/tests/test_linalg.py
+++ b/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))