From 992ac445b936c8da7f90789cfd917ad9727440e7 Mon Sep 17 00:00:00 2001
From: Raymond Hettinger <raymond.hettinger@gmail.com>
Date: Sat, 11 Aug 2018 11:48:58 -0700
Subject: [PATCH 1/4] Small speed-up (approx 9ns) by preallocating space on the
 stack.

Also, put the result computation in an if-else chain to show
parallel structure and to eliminate the need for a "done" target.
---
 Modules/mathmodule.c | 77 ++++++++++++++++++++++++++------------------
 1 file changed, 46 insertions(+), 31 deletions(-)

diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c
index cd390f240da22cc..0c5d7b127f51280 100644
--- a/Modules/mathmodule.c
+++ b/Modules/mathmodule.c
@@ -2057,7 +2057,7 @@ scaled_vector_squared(Py_ssize_t n, double *vec, double max)
         return 0.0;
     }
     assert(n > 0);
-    for (i=0 ; i<n-1 ; i++) {
+    for (i=0 ; i < n-1 ; i++) {
         x = vec[i];
         if (x == max) {
             x = vec[n-1];
@@ -2074,6 +2074,8 @@ scaled_vector_squared(Py_ssize_t n, double *vec, double max)
     return csum;
 }
 
+#define NUM_STACK_ELEMS 16
+
 /*[clinic input]
 math.dist
 
@@ -2095,11 +2097,12 @@ 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 x, px, qx, result;
     Py_ssize_t i, m, n;
     int found_nan = 0;
+    double diffs_on_stack[NUM_STACK_ELEMS];
+    double *diffs = diffs_on_stack;
 
     m = PyTuple_GET_SIZE(p);
     n = PyTuple_GET_SIZE(q);
@@ -2109,22 +2112,22 @@ math_dist_impl(PyObject *module, PyObject *p, PyObject *q)
         return NULL;
 
     }
-    diffs = (double *) PyObject_Malloc(n * sizeof(double));
-    if (diffs == NULL) {
-        return NULL;
+    if (n > NUM_STACK_ELEMS) {
+        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;
+            goto error_exit;
         }
         item = PyTuple_GET_ITEM(q, i);
         qx = PyFloat_AsDouble(item);
         if (qx == -1.0 && PyErr_Occurred()) {
-            PyObject_Free(diffs);
-            return NULL;
+            goto error_exit;
         }
         x = fabs(px - qx);
         diffs[i] = x;
@@ -2135,17 +2138,21 @@ math_dist_impl(PyObject *module, PyObject *p, PyObject *q)
     }
     if (Py_IS_INFINITY(max)) {
         result = max;
-        goto done;
-    }
-    if (found_nan) {
+    } else if (found_nan) {
         result = Py_NAN;
-        goto done;
+    } else {
+        result = max * sqrt(scaled_vector_squared(n, diffs, max));
+    }
+    if (diffs != diffs_on_stack) {
+        PyObject_Free(diffs);
     }
-    result = max * sqrt(scaled_vector_squared(n, diffs, max));
-
-  done:
-    PyObject_Free(diffs);
     return PyFloat_FromDouble(result);
+
+  error_exit:
+    if (diffs != diffs_on_stack) {
+        PyObject_Free(diffs);
+    }
+    return NULL;
 }
 
 /* AC: cannot convert yet, waiting for *args support */
@@ -2154,21 +2161,23 @@ math_hypot(PyObject *self, PyObject *args)
 {
     Py_ssize_t i, n;
     PyObject *item;
-    double *coordinates;
     double max = 0.0;
     double x, result;
     int found_nan = 0;
+    double coord_on_stack[NUM_STACK_ELEMS];
+    double *coordinates = coord_on_stack;
 
     n = PyTuple_GET_SIZE(args);
-    coordinates = (double *) PyObject_Malloc(n * sizeof(double));
-    if (coordinates == NULL)
-        return NULL;
+    if (n > NUM_STACK_ELEMS) {
+        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;
+            goto error_exit;
         }
         x = fabs(x);
         coordinates[i] = x;
@@ -2179,19 +2188,25 @@ math_hypot(PyObject *self, PyObject *args)
     }
     if (Py_IS_INFINITY(max)) {
         result = max;
-        goto done;
-    }
-    if (found_nan) {
+    } else if (found_nan) {
         result = Py_NAN;
-        goto done;
+    } else {
+        result = max * sqrt(scaled_vector_squared(n, coordinates, max));
+    }
+    if (coordinates != coord_on_stack) {
+        PyObject_Free(coordinates);
     }
-    result = max * sqrt(scaled_vector_squared(n, coordinates, max));
-
-  done:
-    PyObject_Free(coordinates);
     return PyFloat_FromDouble(result);
+
+  error_exit:
+    if (coordinates != coord_on_stack) {
+        PyObject_Free(coordinates);
+    }
+    return NULL;
 }
 
+#undef NUM_STACK_ELEMS
+
 PyDoc_STRVAR(math_hypot_doc,
              "hypot(*coordinates) -> value\n\n\
 Multidimensional Euclidean distance from the origin to a point.\n\

From 1fda3c3d07e03b68364e099c4dd9a4ffe8b7bc02 Mon Sep 17 00:00:00 2001
From: Raymond Hettinger <raymond.hettinger@gmail.com>
Date: Sat, 11 Aug 2018 12:52:21 -0700
Subject: [PATCH 2/4] Additional factoring to improve clarity and reduce
 redundancy

---
 Modules/mathmodule.c | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c
index 0c5d7b127f51280..007a80f527581df 100644
--- a/Modules/mathmodule.c
+++ b/Modules/mathmodule.c
@@ -2035,7 +2035,7 @@ math_fmod_impl(PyObject *module, double x, double y)
 Given an *n* length *vec* of non-negative, non-nan, non-inf values
 where *max* is the largest value in the vector, compute:
 
-    sum((x / max) ** 2 for x in vec)
+    max * sqrt(sum((x / max) ** 2 for x in vec))
 
 When a maximum value is found, it is swapped to the end.  This
 lets us skip one loop iteration and just add 1.0 at the end.
@@ -2048,7 +2048,7 @@ fractional round-off error for the most recent addition.
 */
 
 static inline double
-scaled_vector_squared(Py_ssize_t n, double *vec, double max)
+vector_norm(Py_ssize_t n, double *vec, double max)
 {
     double x, csum = 0.0, oldcsum, frac = 0.0;
     Py_ssize_t i;
@@ -2071,7 +2071,7 @@ scaled_vector_squared(Py_ssize_t n, double *vec, double max)
     }
     assert(vec[n-1] == max);
     csum += 1.0 - frac;
-    return csum;
+    return max * sqrt(csum);
 }
 
 #define NUM_STACK_ELEMS 16
@@ -2141,7 +2141,7 @@ math_dist_impl(PyObject *module, PyObject *p, PyObject *q)
     } else if (found_nan) {
         result = Py_NAN;
     } else {
-        result = max * sqrt(scaled_vector_squared(n, diffs, max));
+        result = vector_norm(n, diffs, max);
     }
     if (diffs != diffs_on_stack) {
         PyObject_Free(diffs);
@@ -2191,7 +2191,7 @@ math_hypot(PyObject *self, PyObject *args)
     } else if (found_nan) {
         result = Py_NAN;
     } else {
-        result = max * sqrt(scaled_vector_squared(n, coordinates, max));
+        result = vector_norm(n, coordinates, max);
     }
     if (coordinates != coord_on_stack) {
         PyObject_Free(coordinates);

From ac3dea31cbd89ec1ea526412dd4fbcfeeac4eda1 Mon Sep 17 00:00:00 2001
From: Raymond Hettinger <raymond.hettinger@gmail.com>
Date: Sat, 11 Aug 2018 16:25:11 -0700
Subject: [PATCH 3/4] Refactor by moving special value handling inside the
 helper function.

---
 Modules/mathmodule.c | 29 +++++++++++++----------------
 1 file changed, 13 insertions(+), 16 deletions(-)

diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c
index 007a80f527581df..b63655ef3ebfef4 100644
--- a/Modules/mathmodule.c
+++ b/Modules/mathmodule.c
@@ -2032,7 +2032,7 @@ math_fmod_impl(PyObject *module, double x, double y)
 }
 
 /*
-Given an *n* length *vec* of non-negative, non-nan, non-inf values
+Given an *n* length *vec* of non-negative values
 where *max* is the largest value in the vector, compute:
 
     max * sqrt(sum((x / max) ** 2 for x in vec))
@@ -2045,14 +2045,23 @@ Kahan summation is used to improve accuracy.  The *csum*
 variable tracks the cumulative sum and *frac* tracks
 fractional round-off error for the most recent addition.
 
+The *max* variable must be present in *vec* or equal to 0.0
+when n==0.  Likewise, *max* will be INF if an infinity is
+present in the vec.  The *found_nan* variable must be true
+if some member of the *vec* is a NaN.
 */
 
 static inline double
-vector_norm(Py_ssize_t n, double *vec, double max)
+vector_norm(Py_ssize_t n, double *vec, double max, int found_nan)
 {
     double x, csum = 0.0, oldcsum, frac = 0.0;
     Py_ssize_t i;
 
+    if (Py_IS_INFINITY(max)) {
+        return max;
+    } else if (found_nan) {
+        return Py_NAN;
+    }
     if (max == 0.0) {
         return 0.0;
     }
@@ -2136,13 +2145,7 @@ math_dist_impl(PyObject *module, PyObject *p, PyObject *q)
             max = x;
         }
     }
-    if (Py_IS_INFINITY(max)) {
-        result = max;
-    } else if (found_nan) {
-        result = Py_NAN;
-    } else {
-        result = vector_norm(n, diffs, max);
-    }
+    result = vector_norm(n, diffs, max, found_nan);
     if (diffs != diffs_on_stack) {
         PyObject_Free(diffs);
     }
@@ -2186,13 +2189,7 @@ math_hypot(PyObject *self, PyObject *args)
             max = x;
         }
     }
-    if (Py_IS_INFINITY(max)) {
-        result = max;
-    } else if (found_nan) {
-        result = Py_NAN;
-    } else {
-        result = vector_norm(n, coordinates, max);
-    }
+    result = vector_norm(n, coordinates, max, found_nan);
     if (coordinates != coord_on_stack) {
         PyObject_Free(coordinates);
     }

From 60986636b7c2e68dea2edc89782e2db2f2c6e895 Mon Sep 17 00:00:00 2001
From: Raymond Hettinger <raymond.hettinger@gmail.com>
Date: Sat, 11 Aug 2018 16:45:51 -0700
Subject: [PATCH 4/4] Improve comment

---
 Modules/mathmodule.c | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c
index b63655ef3ebfef4..896cf8dcd5bdeb7 100644
--- a/Modules/mathmodule.c
+++ b/Modules/mathmodule.c
@@ -2045,10 +2045,12 @@ Kahan summation is used to improve accuracy.  The *csum*
 variable tracks the cumulative sum and *frac* tracks
 fractional round-off error for the most recent addition.
 
-The *max* variable must be present in *vec* or equal to 0.0
-when n==0.  Likewise, *max* will be INF if an infinity is
-present in the vec.  The *found_nan* variable must be true
-if some member of the *vec* is a NaN.
+The value of the *max* variable must be present in *vec*
+or should equal to 0.0 when n==0.  Likewise, *max* will
+be INF if an infinity is present in the vec.
+
+The *found_nan* variable indicates whether some member of
+the *vec* is a NaN.
 */
 
 static inline double
@@ -2059,7 +2061,8 @@ vector_norm(Py_ssize_t n, double *vec, double max, int found_nan)
 
     if (Py_IS_INFINITY(max)) {
         return max;
-    } else if (found_nan) {
+    }
+    if (found_nan) {
         return Py_NAN;
     }
     if (max == 0.0) {