From e085c10cc334db9ccafbc8c35590f30e31e47584 Mon Sep 17 00:00:00 2001
From: Serhiy Storchaka <storchaka@gmail.com>
Date: Mon, 6 Mar 2017 15:08:05 +0200
Subject: [PATCH 1/2] bpo-26121: Use C library implementation for math
 functions: tgamma(), lgamma(), erf() and erfc().

---
 Doc/whatsnew/3.7.rst |  5 +++++
 Misc/NEWS            |  3 +++
 Modules/mathmodule.c | 37 ++++++++++++++++++++++++++++++++++++-
 3 files changed, 44 insertions(+), 1 deletion(-)

diff --git a/Doc/whatsnew/3.7.rst b/Doc/whatsnew/3.7.rst
index 5c5ca147e3adf5..6f284ae2287bb1 100644
--- a/Doc/whatsnew/3.7.rst
+++ b/Doc/whatsnew/3.7.rst
@@ -127,6 +127,11 @@ Optimizations
   in method calls being faster up to 20%.
   (Contributed by Yury Selivanov and INADA Naoki in :issue:`26110`.)
 
+* Fast implementation from standard C library is now used for functions
+  :func:`~math.tgamma`, :func:`~math.lgamma`, :func:`~math.erf` and
+  :func:`~math.erfc` in the :mod:`math` module.
+  (Contributed by Serhiy Storchaka in :issue:`26121`.)
+
 
 Build and C API Changes
 =======================
diff --git a/Misc/NEWS b/Misc/NEWS
index 81b02fdd4a48b3..fa08732dec0a31 100644
--- a/Misc/NEWS
+++ b/Misc/NEWS
@@ -266,6 +266,9 @@ Extension Modules
 Library
 -------
 
+- bpo-26121: Use C library implementation for math functions:
+  tgamma(), lgamma(), erf() and erfc().
+
 - bpo-29623: Allow use of path-like object as a single argument in
   ConfigParser.read().  Patch by David Ellis.
 
diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c
index 8bd38d0eb865c5..7380b099bbbdd5 100644
--- a/Modules/mathmodule.c
+++ b/Modules/mathmodule.c
@@ -74,6 +74,8 @@ static const double pi = 3.141592653589793238462643383279502884197;
 static const double sqrtpi = 1.772453850905516027298167483341145182798;
 static const double logpi = 1.144729885849400174143427351353058711647;
 
+#if !defined(HAVE_TGAMMA) || !defined(HAVE_LGAMMA)
+
 static double
 sinpi(double x)
 {
@@ -230,6 +232,7 @@ lanczos_sum(double x)
     }
     return num/den;
 }
+#endif /* !defined(HAVE_TGAMMA) || !defined(HAVE_LGAMMA) */
 
 /* Constant for +infinity, generated in the same way as float('inf'). */
 
@@ -263,6 +266,14 @@ m_nan(void)
 static double
 m_tgamma(double x)
 {
+#ifdef HAVE_TGAMMA
+    if (x == 0.0) {
+        errno = EDOM;
+        /* tgamma(+-0.0) = +-inf, divide-by-zero */
+        return copysign(Py_HUGE_VAL, x);
+    }
+    return tgamma(x);
+#else
     double absx, r, y, z, sqrtpow;
 
     /* special cases */
@@ -354,6 +365,7 @@ m_tgamma(double x)
     if (Py_IS_INFINITY(r))
         errno = ERANGE;
     return r;
+#endif
 }
 
 /*
@@ -364,7 +376,17 @@ m_tgamma(double x)
 static double
 m_lgamma(double x)
 {
-    double r, absx;
+    double r;
+
+#ifdef HAVE_LGAMMA
+    r = lgamma(x);
+    if (errno == ERANGE && x == floor(x) && x <= 0.0) {
+        errno = EDOM;           /* lgamma(n) = inf, divide-by-zero for */
+        return Py_HUGE_VAL;     /* integers n <= 0 */
+    }
+    return r;
+#else
+    double absx;
 
     /* special cases */
     if (!Py_IS_FINITE(x)) {
@@ -402,8 +424,11 @@ m_lgamma(double x)
     if (Py_IS_INFINITY(r))
         errno = ERANGE;
     return r;
+#endif
 }
 
+#if !defined(HAVE_ERF) || !defined(HAVE_ERFC)
+
 /*
    Implementations of the error function erf(x) and the complementary error
    function erfc(x).
@@ -513,11 +538,16 @@ m_erfc_contfrac(double x)
     return result;
 }
 
+#endif /* !defined(HAVE_ERF) || !defined(HAVE_ERFC) */
+
 /* Error function erf(x), for general x */
 
 static double
 m_erf(double x)
 {
+#ifdef HAVE_ERF
+    return erf(x);
+#else
     double absx, cf;
 
     if (Py_IS_NAN(x))
@@ -529,6 +559,7 @@ m_erf(double x)
         cf = m_erfc_contfrac(absx);
         return x > 0.0 ? 1.0 - cf : cf - 1.0;
     }
+#endif
 }
 
 /* Complementary error function erfc(x), for general x. */
@@ -536,6 +567,9 @@ m_erf(double x)
 static double
 m_erfc(double x)
 {
+#ifdef HAVE_ERFC
+    return erfc(x);
+#else
     double absx, cf;
 
     if (Py_IS_NAN(x))
@@ -547,6 +581,7 @@ m_erfc(double x)
         cf = m_erfc_contfrac(absx);
         return x > 0.0 ? cf : 2.0 - cf;
     }
+#endif
 }
 
 /*

From bfe71106306ab1fe655c6267c75c319dd915f9f8 Mon Sep 17 00:00:00 2001
From: Serhiy Storchaka <storchaka@gmail.com>
Date: Sat, 11 Mar 2017 23:10:18 +0200
Subject: [PATCH 2/2] Don't use tgamma() and lgamma() from libc on OS X.

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

diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c
index 7380b099bbbdd5..bfc028edddfa83 100644
--- a/Modules/mathmodule.c
+++ b/Modules/mathmodule.c
@@ -74,7 +74,16 @@ static const double pi = 3.141592653589793238462643383279502884197;
 static const double sqrtpi = 1.772453850905516027298167483341145182798;
 static const double logpi = 1.144729885849400174143427351353058711647;
 
-#if !defined(HAVE_TGAMMA) || !defined(HAVE_LGAMMA)
+#ifndef __APPLE__
+#  ifdef HAVE_TGAMMA
+#    define USE_TGAMMA
+#  endif
+#  ifdef HAVE_LGAMMA
+#    define USE_LGAMMA
+#  endif
+#endif
+
+#if !defined(USE_TGAMMA) || !defined(USE_LGAMMA)
 
 static double
 sinpi(double x)
@@ -232,7 +241,7 @@ lanczos_sum(double x)
     }
     return num/den;
 }
-#endif /* !defined(HAVE_TGAMMA) || !defined(HAVE_LGAMMA) */
+#endif /* !defined(USE_TGAMMA) || !defined(USE_LGAMMA) */
 
 /* Constant for +infinity, generated in the same way as float('inf'). */
 
@@ -266,7 +275,7 @@ m_nan(void)
 static double
 m_tgamma(double x)
 {
-#ifdef HAVE_TGAMMA
+#ifdef USE_TGAMMA
     if (x == 0.0) {
         errno = EDOM;
         /* tgamma(+-0.0) = +-inf, divide-by-zero */
@@ -378,7 +387,7 @@ m_lgamma(double x)
 {
     double r;
 
-#ifdef HAVE_LGAMMA
+#ifdef USE_LGAMMA
     r = lgamma(x);
     if (errno == ERANGE && x == floor(x) && x <= 0.0) {
         errno = EDOM;           /* lgamma(n) = inf, divide-by-zero for */