bpo-33089: Multidimensional math.hypot() (GH-8474)

5 files changed, 143 insertions, 81 deletions

diff --git a/Doc/library/math.rst b/Doc/library/math.rst
index 33aec57..e60d029 100644
--- a/Doc/library/math.rst
+++ b/Doc/library/math.rst
@@ -330,10 +330,19 @@ Trigonometric functions
    Return the cosine of *x* radians.
 
 
-.. function:: hypot(x, y)
+.. function:: hypot(*coordinates)
 
-   Return the Euclidean norm, ``sqrt(x*x + y*y)``. This is the length of the vector
-   from the origin to point ``(x, y)``.
+   Return the Euclidean norm, ``sqrt(sum(x**2 for x in coordinates))``.
+   This is the length of the vector from the origin to the point
+   given by the coordinates.
+
+   For a two dimensional point ``(x, y)``, this is equivalent to computing
+   the hypotenuse of a right triangle using the Pythagorean theorem,
+   ``sqrt(x*x + y*y)``.
+
+   .. versionchanged:: 3.8
+      Added support for n-dimensional points. Formerly, only the two
+      dimensional case was supported.
 
 
 .. function:: sin(x)
diff --git a/Lib/test/test_math.py b/Lib/test/test_math.py
index 44785d3..4843899 100644
--- a/Lib/test/test_math.py
+++ b/Lib/test/test_math.py
@@ -16,6 +16,7 @@ NAN = float('nan')
 INF = float('inf')
 NINF = float('-inf')
 FLOAT_MAX = sys.float_info.max
+FLOAT_MIN = sys.float_info.min
 
 # detect evidence of double-rounding: fsum is not always correctly
 # rounded on machines that suffer from double rounding.
@@ -720,16 +721,71 @@ class MathTests(unittest.TestCase):
         self.assertEqual(gcd(MyIndexable(120), MyIndexable(84)), 12)
 
     def testHypot(self):
-        self.assertRaises(TypeError, math.hypot)
-        self.ftest('hypot(0,0)', math.hypot(0,0), 0)
-        self.ftest('hypot(3,4)', math.hypot(3,4), 5)
-        self.assertEqual(math.hypot(NAN, INF), INF)
-        self.assertEqual(math.hypot(INF, NAN), INF)
-        self.assertEqual(math.hypot(NAN, NINF), INF)
-        self.assertEqual(math.hypot(NINF, NAN), INF)
-        self.assertRaises(OverflowError, math.hypot, FLOAT_MAX, FLOAT_MAX)
-        self.assertTrue(math.isnan(math.hypot(1.0, NAN)))
-        self.assertTrue(math.isnan(math.hypot(NAN, -2.0)))
+        from decimal import Decimal
+        from fractions import Fraction
+
+        hypot = math.hypot
+
+        # Test different numbers of arguments (from zero to five)
+        # against a straightforward pure python implementation
+        args = math.e, math.pi, math.sqrt(2.0), math.gamma(3.5), math.sin(2.1)
+        for i in range(len(args)+1):
+            self.assertAlmostEqual(
+                hypot(*args[:i]),
+                math.sqrt(sum(s**2 for s in args[:i]))
+            )
+
+        # Test allowable types (those with __float__)
+        self.assertEqual(hypot(12.0, 5.0), 13.0)
+        self.assertEqual(hypot(12, 5), 13)
+        self.assertEqual(hypot(Decimal(12), Decimal(5)), 13)
+        self.assertEqual(hypot(Fraction(12, 32), Fraction(5, 32)), Fraction(13, 32))
+        self.assertEqual(hypot(bool(1), bool(0), bool(1), bool(1)), math.sqrt(3))
+
+        # Test corner cases
+        self.assertEqual(hypot(0.0, 0.0), 0.0)     # Max input is zero
+        self.assertEqual(hypot(-10.5), 10.5)       # Negative input
+        self.assertEqual(hypot(), 0.0)             # Negative input
+        self.assertEqual(1.0,
+            math.copysign(1.0, hypot(-0.0))        # Convert negative zero to positive zero
+        )
+
+        # Test handling of bad arguments
+        with self.assertRaises(TypeError):         # Reject keyword args
+            hypot(x=1)
+        with self.assertRaises(TypeError):         # Reject values without __float__
+            hypot(1.1, 'string', 2.2)
+
+        # Any infinity gives positive infinity.
+        self.assertEqual(hypot(INF), INF)
+        self.assertEqual(hypot(0, INF), INF)
+        self.assertEqual(hypot(10, INF), INF)
+        self.assertEqual(hypot(-10, INF), INF)
+        self.assertEqual(hypot(NAN, INF), INF)
+        self.assertEqual(hypot(INF, NAN), INF)
+        self.assertEqual(hypot(NINF, NAN), INF)
+        self.assertEqual(hypot(NAN, NINF), INF)
+        self.assertEqual(hypot(-INF, INF), INF)
+        self.assertEqual(hypot(-INF, -INF), INF)
+        self.assertEqual(hypot(10, -INF), INF)
+
+        # If no infinity, any NaN gives a Nan.
+        self.assertTrue(math.isnan(hypot(NAN)))
+        self.assertTrue(math.isnan(hypot(0, NAN)))
+        self.assertTrue(math.isnan(hypot(NAN, 10)))
+        self.assertTrue(math.isnan(hypot(10, NAN)))
+        self.assertTrue(math.isnan(hypot(NAN, NAN)))
+        self.assertTrue(math.isnan(hypot(NAN)))
+
+        # Verify scaling for extremely large values
+        fourthmax = FLOAT_MAX / 4.0
+        for n in range(32):
+            self.assertEqual(hypot(*([fourthmax]*n)), fourthmax * math.sqrt(n))
+
+        # Verify scaling for extremely small values
+        for exp in range(32):
+            scale = FLOAT_MIN / 2.0 ** exp
+            self.assertEqual(math.hypot(4*scale, 3*scale), 5*scale)
 
     def testLdexp(self):
         self.assertRaises(TypeError, math.ldexp)
diff --git a/Misc/NEWS.d/next/Library/2018-07-25-22-38-54.bpo-33089.C3CB7e.rst b/Misc/NEWS.d/next/Library/2018-07-25-22-38-54.bpo-33089.C3CB7e.rst
new file mode 100644
index 0000000..83152a7
--- /dev/null
+++ b/Misc/NEWS.d/next/Library/2018-07-25-22-38-54.bpo-33089.C3CB7e.rst
@@ -0,0 +1 @@
+Enhanced math.hypot() to support more than two dimensions.
diff --git a/Modules/clinic/mathmodule.c.h b/Modules/clinic/mathmodule.c.h
index ffebb07..3a487bd 100644
--- a/Modules/clinic/mathmodule.c.h
+++ b/Modules/clinic/mathmodule.c.h
@@ -269,35 +269,6 @@ exit:
     return return_value;
 }
 
-PyDoc_STRVAR(math_hypot__doc__,
-"hypot($module, x, y, /)\n"
-"--\n"
-"\n"
-"Return the Euclidean distance, sqrt(x*x + y*y).");
-
-#define MATH_HYPOT_METHODDEF    \
-    {"hypot", (PyCFunction)math_hypot, METH_FASTCALL, math_hypot__doc__},
-
-static PyObject *
-math_hypot_impl(PyObject *module, double x, double y);
-
-static PyObject *
-math_hypot(PyObject *module, PyObject *const *args, Py_ssize_t nargs)
-{
-    PyObject *return_value = NULL;
-    double x;
-    double y;
-
-    if (!_PyArg_ParseStack(args, nargs, "dd:hypot",
-        &x, &y)) {
-        goto exit;
-    }
-    return_value = math_hypot_impl(module, x, y);
-
-exit:
-    return return_value;
-}
-
 PyDoc_STRVAR(math_pow__doc__,
 "pow($module, x, y, /)\n"
 "--\n"
@@ -516,4 +487,4 @@ math_isclose(PyObject *module, PyObject *const *args, Py_ssize_t nargs, PyObject
 exit:
     return return_value;
 }
-/*[clinic end generated code: output=e554bad553045546 input=a9049054013a1b77]*/
+/*[clinic end generated code: output=1c35516a10443902 input=a9049054013a1b77]*/
diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c
index 5d9fe5a..726fc56 100644
--- a/Modules/mathmodule.c
+++ b/Modules/mathmodule.c
@@ -2031,49 +2031,74 @@ math_fmod_impl(PyObject *module, double x, double y)
         return PyFloat_FromDouble(r);
 }
 
-
-/*[clinic input]
-math.hypot
-
-    x: double
-    y: double
-    /
-
-Return the Euclidean distance, sqrt(x*x + y*y).
-[clinic start generated code]*/
-
+/* AC: cannot convert yet, waiting for *args support */
 static PyObject *
-math_hypot_impl(PyObject *module, double x, double y)
-/*[clinic end generated code: output=b7686e5be468ef87 input=7f8eea70406474aa]*/
+math_hypot(PyObject *self, PyObject *args)
 {
-    double r;
-    /* hypot(x, +/-Inf) returns Inf, even if x is a NaN. */
-    if (Py_IS_INFINITY(x))
-        return PyFloat_FromDouble(fabs(x));
-    if (Py_IS_INFINITY(y))
-        return PyFloat_FromDouble(fabs(y));
-    errno = 0;
-    PyFPE_START_PROTECT("in math_hypot", return 0);
-    r = hypot(x, y);
-    PyFPE_END_PROTECT(r);
-    if (Py_IS_NAN(r)) {
-        if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
-            errno = EDOM;
-        else
-            errno = 0;
+    Py_ssize_t i, n;
+    PyObject *item;
+    double *coordinates;
+    double max = 0.0;
+    double csum = 0.0;
+    double x, result;
+    int found_nan = 0;
+
+    n = PyTuple_GET_SIZE(args);
+    coordinates = (double *) PyObject_Malloc(n * sizeof(double));
+    if (coordinates == NULL)
+        return NULL;
+    for (i=0 ; i<n ; i++) {
+        item = PyTuple_GET_ITEM(args, i);
+        x = PyFloat_AsDouble(item);
+        if (x == -1.0 && PyErr_Occurred()) {
+            PyObject_Free(coordinates);
+            return NULL;
+        }
+        x = fabs(x);
+        coordinates[i] = x;
+        found_nan |= Py_IS_NAN(x);
+        if (x > max) {
+            max = x;
+        }
     }
-    else if (Py_IS_INFINITY(r)) {
-        if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
-            errno = ERANGE;
-        else
-            errno = 0;
+    if (Py_IS_INFINITY(max)) {
+        result = max;
+        goto done;
     }
-    if (errno && is_error(r))
-        return NULL;
-    else
-        return PyFloat_FromDouble(r);
+    if (found_nan) {
+        result = Py_NAN;
+        goto done;
+    }
+    if (max == 0.0) {
+        result = 0.0;
+        goto done;
+    }
+    for (i=0 ; i<n ; i++) {
+        x = coordinates[i] / max;
+        csum += x * x;
+    }
+    result = max * sqrt(csum);
+
+  done:
+    PyObject_Free(coordinates);
+    return PyFloat_FromDouble(result);
 }
 
+PyDoc_STRVAR(math_hypot_doc,
+             "hypot(*coordinates) -> value\n\n\
+Multidimensional Euclidean distance from the origin to a point.\n\
+\n\
+Roughly equivalent to:\n\
+    sqrt(sum(x**2 for x in coordinates))\n\
+\n\
+For a two dimensional point (x, y), gives the hypotenuse\n\
+using the Pythagorean theorem:  sqrt(x*x + y*y).\n\
+\n\
+For example, the hypotenuse of a 3/4/5 right triangle is:\n\
+\n\
+    >>> hypot(3.0, 4.0)\n\
+    5.0\n\
+");
 
 /* pow can't use math_2, but needs its own wrapper: the problem is
    that an infinite result can arise either as a result of overflow
@@ -2345,7 +2370,7 @@ static PyMethodDef math_methods[] = {
     MATH_FSUM_METHODDEF
     {"gamma",           math_gamma,     METH_O,         math_gamma_doc},
     MATH_GCD_METHODDEF
-    MATH_HYPOT_METHODDEF
+    {"hypot",           math_hypot,     METH_VARARGS,   math_hypot_doc},
     MATH_ISCLOSE_METHODDEF
     MATH_ISFINITE_METHODDEF
     MATH_ISINF_METHODDEF