From d86b785b6d67e2bac78a5d14a01cd806b3151b31 Mon Sep 17 00:00:00 2001
From: Luc Maisonobe <luc@apache.org>
Date: Sun, 23 Jan 2011 22:37:50 +0000
Subject: [PATCH] added FastMath.hypot JIRA: MATH-478

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1062559 13f79535-47bb-0310-9956-ffa450edef68
---
 .../apache/commons/math/util/FastMath.java    | 35 +++++++++++++++++--
 src/site/xdoc/changes.xml                     |  6 ++--
 2 files changed, 36 insertions(+), 5 deletions(-)

diff --git a/src/main/java/org/apache/commons/math/util/FastMath.java b/src/main/java/org/apache/commons/math/util/FastMath.java
index 0e55eef66..d22232752 100644
--- a/src/main/java/org/apache/commons/math/util/FastMath.java
+++ b/src/main/java/org/apache/commons/math/util/FastMath.java
@@ -3847,8 +3847,39 @@ public class FastMath {
      * @param y a value
      * @return sqrt(<i>x</i><sup>2</sup>&nbsp;+<i>y</i><sup>2</sup>)
      */
-    public static double hypot(double x, double y) {
-        return StrictMath.hypot(x, y); // TODO provide our own implementation
+    public static double hypot(final double x, final double y) {
+        if (Double.isInfinite(x) || Double.isInfinite(y)) {
+            return Double.POSITIVE_INFINITY;
+        } else if (Double.isNaN(x) || Double.isNaN(y)) {
+            return Double.NaN;
+        } else {
+
+            final int expX = getExponent(x);
+            final int expY = getExponent(y);
+            if (expX > expY + 27) {
+                // y is neglectible with respect to x
+                return abs(x);
+            } else if (expY > expX + 27) {
+                // x is neglectible with respect to y
+                return abs(y);
+            } else {
+
+                // find an intermediate scale to avoid both overflow and underflow
+                final int middleExp = (expX + expY) / 2;
+
+                // scale parameters without losing precision
+                final double scaledX = scalb(x, -middleExp);
+                final double scaledY = scalb(y, -middleExp);
+
+                // compute scaled hypotenuse
+                final double scaledH = sqrt(scaledX * scaledX + scaledY * scaledY);
+
+                // remove scaling
+                return scalb(scaledH, middleExp);
+
+            }
+
+        }
     }
 
     /**
diff --git a/src/site/xdoc/changes.xml b/src/site/xdoc/changes.xml
index f5675b96f..951d58ec4 100644
--- a/src/site/xdoc/changes.xml
+++ b/src/site/xdoc/changes.xml
@@ -192,9 +192,9 @@ The <action> type attribute can be add,update,fix,remove.
       </action>
       <action dev="luc" type="fix" issue="MATH-478">
           FastMath is not an exact replacement for StrictMath
-          (partially fixed) added nextAfter(double, double) and nextAfter(float,double)
-          (beware of the strange double second argument) so that they handle
-          special values in the way as StrictMath
+          (partially fixed) added hypot(double, double), nextAfter(double, double)
+          and nextAfter(float,double) (beware of the strange double second argument)
+          so that they handle special values in the way as StrictMath
       </action>
       <action dev="luc" type="fix" issue="MATH-497">
           FastMath is not an exact replacement for StrictMath