bpo-33089: Add math.dist() for computing the Euclidean distance between two points (GH-8561)

5 files changed, 236 insertions, 1 deletions

diff --git a/Doc/library/math.rst b/Doc/library/math.rst
index e60d029..76226c2 100644
--- a/Doc/library/math.rst
+++ b/Doc/library/math.rst
@@ -330,6 +330,18 @@ Trigonometric functions
    Return the cosine of *x* radians.
 
 
+.. function:: dist(p, q)
+
+   Return the Euclidean distance between two points *p* and *q*, each
+   given as a tuple of coordinates.  The two tuples must be the same size.
+
+   Roughly equivalent to::
+
+       sqrt(sum((px - qx) ** 2.0 for px, qx in zip(p, q)))
+
+   .. versionadded:: 3.8
+
+
 .. function:: hypot(*coordinates)
 
    Return the Euclidean norm, ``sqrt(sum(x**2 for x in coordinates))``.
diff --git a/Lib/test/test_math.py b/Lib/test/test_math.py
index 4843899..448110b 100644
--- a/Lib/test/test_math.py
+++ b/Lib/test/test_math.py
@@ -4,9 +4,11 @@
 from test.support import run_unittest, verbose, requires_IEEE_754
 from test import support
 import unittest
+import itertools
 import math
 import os
 import platform
+import random
 import struct
 import sys
 import sysconfig
@@ -787,6 +789,107 @@ class MathTests(unittest.TestCase):
             scale = FLOAT_MIN / 2.0 ** exp
             self.assertEqual(math.hypot(4*scale, 3*scale), 5*scale)
 
+    def testDist(self):
+        from decimal import Decimal as D
+        from fractions import Fraction as F
+
+        dist = math.dist
+        sqrt = math.sqrt
+
+        # Simple exact case
+        self.assertEqual(dist((1, 2, 3), (4, 2, -1)), 5.0)
+
+        # Test different numbers of arguments (from zero to nine)
+        # against a straightforward pure python implementation
+        for i in range(9):
+            for j in range(5):
+                p = tuple(random.uniform(-5, 5) for k in range(i))
+                q = tuple(random.uniform(-5, 5) for k in range(i))
+                self.assertAlmostEqual(
+                    dist(p, q),
+                    sqrt(sum((px - qx) ** 2.0 for px, qx in zip(p, q)))
+                )
+
+        # Test allowable types (those with __float__)
+        self.assertEqual(dist((14.0, 1.0), (2.0, -4.0)), 13.0)
+        self.assertEqual(dist((14, 1), (2, -4)), 13)
+        self.assertEqual(dist((D(14), D(1)), (D(2), D(-4))), D(13))
+        self.assertEqual(dist((F(14, 32), F(1, 32)), (F(2, 32), F(-4, 32))),
+                         F(13, 32))
+        self.assertEqual(dist((True, True, False, True, False),
+                              (True, False, True, True, False)),
+                         sqrt(2.0))
+
+        # Test corner cases
+        self.assertEqual(dist((13.25, 12.5, -3.25),
+                              (13.25, 12.5, -3.25)),
+                         0.0)                      # Distance with self is zero
+        self.assertEqual(dist((), ()), 0.0)        # Zero-dimensional case
+        self.assertEqual(1.0,                      # Convert negative zero to positive zero
+            math.copysign(1.0, dist((-0.0,), (0.0,)))
+        )
+        self.assertEqual(1.0,                      # Convert negative zero to positive zero
+            math.copysign(1.0, dist((0.0,), (-0.0,)))
+        )
+
+        # Verify tuple subclasses are allowed
+        class T(tuple):     # tuple subclas
+            pass
+        self.assertEqual(dist(T((1, 2, 3)), ((4, 2, -1))), 5.0)
+
+        # Test handling of bad arguments
+        with self.assertRaises(TypeError):         # Reject keyword args
+            dist(p=(1, 2, 3), q=(4, 5, 6))
+        with self.assertRaises(TypeError):         # Too few args
+            dist((1, 2, 3))
+        with self.assertRaises(TypeError):         # Too many args
+            dist((1, 2, 3), (4, 5, 6), (7, 8, 9))
+        with self.assertRaises(TypeError):         # Scalars not allowed
+            dist(1, 2)
+        with self.assertRaises(TypeError):         # Lists not allowed
+            dist([1, 2, 3], [4, 5, 6])
+        with self.assertRaises(TypeError):         # Reject values without __float__
+            dist((1.1, 'string', 2.2), (1, 2, 3))
+        with self.assertRaises(ValueError):        # Check dimension agree
+            dist((1, 2, 3, 4), (5, 6, 7))
+        with self.assertRaises(ValueError):        # Check dimension agree
+            dist((1, 2, 3), (4, 5, 6, 7))
+
+
+        # Verify that the one dimensional case equivalent to abs()
+        for i in range(20):
+            p, q = random.random(), random.random()
+            self.assertEqual(dist((p,), (q,)), abs(p - q))
+
+        # Test special values
+        values = [NINF, -10.5, -0.0, 0.0, 10.5, INF, NAN]
+        for p in itertools.product(values, repeat=3):
+            for q in itertools.product(values, repeat=3):
+                diffs = [px - qx for px, qx in zip(p, q)]
+                if any(map(math.isinf, diffs)):
+                    # Any infinite difference gives positive infinity.
+                    self.assertEqual(dist(p, q), INF)
+                elif any(map(math.isnan, diffs)):
+                    # If no infinity, any NaN gives a Nan.
+                    self.assertTrue(math.isnan(dist(p, q)))
+
+        # Verify scaling for extremely large values
+        fourthmax = FLOAT_MAX / 4.0
+        for n in range(32):
+            p = (fourthmax,) * n
+            q = (0.0,) * n
+            self.assertEqual(dist(p, q), fourthmax * math.sqrt(n))
+            self.assertEqual(dist(q, p), fourthmax * math.sqrt(n))
+
+        # Verify scaling for extremely small values
+        for exp in range(32):
+            scale = FLOAT_MIN / 2.0 ** exp
+            p = (4*scale, 3*scale)
+            q = (0.0, 0.0)
+            self.assertEqual(math.dist(p, q), 5*scale)
+            self.assertEqual(math.dist(q, p), 5*scale)
+
+
     def testLdexp(self):
         self.assertRaises(TypeError, math.ldexp)
         self.ftest('ldexp(0,1)', math.ldexp(0,1), 0)
diff --git a/Misc/NEWS.d/next/Library/2018-07-29-21-53-15.bpo-33089.hxbp3g.rst b/Misc/NEWS.d/next/Library/2018-07-29-21-53-15.bpo-33089.hxbp3g.rst
new file mode 100644
index 0000000..9b253d0
--- /dev/null
+++ b/Misc/NEWS.d/next/Library/2018-07-29-21-53-15.bpo-33089.hxbp3g.rst
@@ -0,0 +1 @@
+Add math.dist() to compute the Euclidean distance between two points.
diff --git a/Modules/clinic/mathmodule.c.h b/Modules/clinic/mathmodule.c.h
index 3a487bd..c4d2786 100644
--- a/Modules/clinic/mathmodule.c.h
+++ b/Modules/clinic/mathmodule.c.h
@@ -269,6 +269,41 @@ exit:
     return return_value;
 }
 
+PyDoc_STRVAR(math_dist__doc__,
+"dist($module, p, q, /)\n"
+"--\n"
+"\n"
+"Return the Euclidean distance between two points p and q.\n"
+"\n"
+"The points should be specified as tuples of coordinates.\n"
+"Both tuples must be the same size.\n"
+"\n"
+"Roughly equivalent to:\n"
+"    sqrt(sum((px - qx) ** 2.0 for px, qx in zip(p, q)))");
+
+#define MATH_DIST_METHODDEF    \
+    {"dist", (PyCFunction)math_dist, METH_FASTCALL, math_dist__doc__},
+
+static PyObject *
+math_dist_impl(PyObject *module, PyObject *p, PyObject *q);
+
+static PyObject *
+math_dist(PyObject *module, PyObject *const *args, Py_ssize_t nargs)
+{
+    PyObject *return_value = NULL;
+    PyObject *p;
+    PyObject *q;
+
+    if (!_PyArg_ParseStack(args, nargs, "O!O!:dist",
+        &PyTuple_Type, &p, &PyTuple_Type, &q)) {
+        goto exit;
+    }
+    return_value = math_dist_impl(module, p, q);
+
+exit:
+    return return_value;
+}
+
 PyDoc_STRVAR(math_pow__doc__,
 "pow($module, x, y, /)\n"
 "--\n"
@@ -487,4 +522,4 @@ math_isclose(PyObject *module, PyObject *const *args, Py_ssize_t nargs, PyObject
 exit:
     return return_value;
 }
-/*[clinic end generated code: output=1c35516a10443902 input=a9049054013a1b77]*/
+/*[clinic end generated code: output=d936137c1189b89b input=a9049054013a1b77]*/
diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c
index 726fc56..8f11f2c 100644
--- a/Modules/mathmodule.c
+++ b/Modules/mathmodule.c
@@ -2031,6 +2031,89 @@ math_fmod_impl(PyObject *module, double x, double y)
         return PyFloat_FromDouble(r);
 }
 
+/*[clinic input]
+math.dist
+
+    p: object(subclass_of='&PyTuple_Type')
+    q: object(subclass_of='&PyTuple_Type')
+    /
+
+Return the Euclidean distance between two points p and q.
+
+The points should be specified as tuples of coordinates.
+Both tuples must be the same size.
+
+Roughly equivalent to:
+    sqrt(sum((px - qx) ** 2.0 for px, qx in zip(p, q)))
+[clinic start generated code]*/
+
+static PyObject *
+math_dist_impl(PyObject *module, PyObject *p, PyObject *q)
+/*[clinic end generated code: output=56bd9538d06bbcfe input=937122eaa5f19272]*/
+{
+    PyObject *item;
+    double *diffs;
+    double max = 0.0;
+    double csum = 0.0;
+    double x, px, qx, result;
+    Py_ssize_t i, m, n;
+    int found_nan = 0;
+
+    m = PyTuple_GET_SIZE(p);
+    n = PyTuple_GET_SIZE(q);
+    if (m != n) {
+        PyErr_SetString(PyExc_ValueError,
+                        "both points must have the same number of dimensions");
+        return NULL;
+
+    }
+    diffs = (double *) PyObject_Malloc(n * sizeof(double));
+    if (diffs == NULL) {
+        return NULL;
+    }
+    for (i=0 ; i<n ; i++) {
+        item = PyTuple_GET_ITEM(p, i);
+        px = PyFloat_AsDouble(item);
+        if (px == -1.0 && PyErr_Occurred()) {
+            PyObject_Free(diffs);
+            return NULL;
+        }
+        item = PyTuple_GET_ITEM(q, i);
+        qx = PyFloat_AsDouble(item);
+        if (qx == -1.0 && PyErr_Occurred()) {
+            PyObject_Free(diffs);
+            return NULL;
+        }
+        x = fabs(px - qx);
+        diffs[i] = x;
+        found_nan |= Py_IS_NAN(x);
+        if (x > max) {
+            max = x;
+        }
+    }
+    if (Py_IS_INFINITY(max)) {
+        result = max;
+        goto done;
+    }
+    if (found_nan) {
+        result = Py_NAN;
+        goto done;
+    }
+    if (max == 0.0) {
+        result = 0.0;
+        goto done;
+    }
+    for (i=0 ; i<n ; i++) {
+        x = diffs[i] / max;
+        csum += x * x;
+    }
+    result = max * sqrt(csum);
+
+  done:
+    PyObject_Free(diffs);
+    return PyFloat_FromDouble(result);
+}
+
 /* AC: cannot convert yet, waiting for *args support */
 static PyObject *
 math_hypot(PyObject *self, PyObject *args)
@@ -2358,6 +2441,7 @@ static PyMethodDef math_methods[] = {
     {"cos",             math_cos,       METH_O,         math_cos_doc},
     {"cosh",            math_cosh,      METH_O,         math_cosh_doc},
     MATH_DEGREES_METHODDEF
+    MATH_DIST_METHODDEF
     {"erf",             math_erf,       METH_O,         math_erf_doc},
     {"erfc",            math_erfc,      METH_O,         math_erfc_doc},
     {"exp",             math_exp,       METH_O,         math_exp_doc},