From 79c4a7fec90e697d1a5896c7883ff24d562bad6d Mon Sep 17 00:00:00 2001
From: Paul Khuong <pvk@pvk.ca>
Date: Sun, 13 Nov 2011 15:00:21 -0500
Subject: [PATCH] More numerically stable %hypot (ABS of complex floats) on
 win32

Based on a patch by Robert Smith.
---
 NEWS                  |    2 ++
 src/code/irrat.lisp   |   21 ++++++++++++++++++---
 tests/float.pure.lisp |    4 ++++
 3 files changed, 24 insertions(+), 3 deletions(-)

diff --git a/NEWS b/NEWS
index 70dcf9b..b041a7a 100644
--- a/NEWS
+++ b/NEWS
@@ -3,6 +3,8 @@ changes relative to sbcl-1.0.53:
   * enchancement: on CHENEYGC targets, SB-KERNEL:MAKE-LISP-OBJ now does
     the same validation of pointer objects as GENCGC does, instead of a
     comparatively weak bounds-check against the heap spaces.
+  * enhancement: on win32, ABS of complex floats guards better against
+    overflows. (lp#888410)
   * bug fix: on 64-bit targets, atomic-incf/aref does index computation
     correctly, even on wide-fixnum builds. (lp#887220)
   * bug fix: (directory "foo/*/*.*") did not follow symlinks in foo/ that
diff --git a/src/code/irrat.lisp b/src/code/irrat.lisp
index 05410a6..93c957b 100644
--- a/src/code/irrat.lisp
+++ b/src/code/irrat.lisp
@@ -127,11 +127,26 @@
 
 #!+win32
 (progn
-  ;; FIXME: libc hypot "computes the sqrt(x*x+y*y) without undue overflow or underflow"
-  ;; ...we just do the stupid simple thing.
+  ;; This is written in a peculiar way to avoid overflow. Note that in
+  ;; sqrt(x^2 + y^2), either square or the sum can overflow.
+  ;;
+  ;; Factoring x^2 out of sqrt(x^2 + y^2) gives us the expression
+  ;; |x|sqrt(1 + (y/x)^2), which, assuming |x| >= |y|, can only overflow
+  ;; if |x| is sufficiently large.
+  ;;
+  ;; The ZEROP test suffices (y is non-negative) to guard against
+  ;; divisions by zero: x >= y > 0.
   (declaim (inline %hypot))
   (defun %hypot (x y)
-    (sqrt (+ (* x x) (* y y)))))
+    (declare (type double-float x y))
+    (let ((x (abs x))
+          (y (abs y)))
+      (when (> y x)
+        (rotatef x y))
+      (if (zerop y)
+          x
+          (let ((y/x (/ y x)))
+            (* x (sqrt (1+ (* y/x y/x)))))))))
 
 ;;;; power functions
 
diff --git a/tests/float.pure.lisp b/tests/float.pure.lisp
index 64f601e..4f28680 100644
--- a/tests/float.pure.lisp
+++ b/tests/float.pure.lisp
@@ -415,3 +415,7 @@
 (with-test (:name :round-single-to-bignum)
   (assert (= (round 1e14) 100000000376832))
   (assert (= (round 1e19) 9999999980506447872)))
+
+(with-test (:name :scaled-%hypot)
+  (assert (<= (abs (complex most-positive-double-float 1d0))
+              (1+ most-positive-double-float))))
-- 
1.7.10.4