diff --git a/Doc/whatsnew/3.7.rst b/Doc/whatsnew/3.7.rst
index 1c737e8..71ae10b 100644
--- a/Doc/whatsnew/3.7.rst
+++ b/Doc/whatsnew/3.7.rst
@@ -133,6 +133,11 @@
   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 c4bdbef..59ec700 100644
--- a/Misc/NEWS
+++ b/Misc/NEWS
@@ -270,6 +270,9 @@
 Library
 -------
 
+- bpo-26121: Use C library implementation for math functions:
+  tgamma(), lgamma(), erf() and erfc().
+
 - bpo-29619: os.stat() and os.DirEntry.inode() now convert inode (st_ino) using
   unsigned integers.
 
diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c
index 8bd38d0..bfc028e 100644
--- a/Modules/mathmodule.c
+++ b/Modules/mathmodule.c
@@ -74,6 +74,17 @@
 static const double sqrtpi = 1.772453850905516027298167483341145182798;
 static const double logpi = 1.144729885849400174143427351353058711647;
 
+#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)
 {
@@ -230,6 +241,7 @@
     }
     return num/den;
 }
+#endif /* !defined(USE_TGAMMA) || !defined(USE_LGAMMA) */
 
 /* Constant for +infinity, generated in the same way as float('inf'). */
 
@@ -263,6 +275,14 @@
 static double
 m_tgamma(double x)
 {
+#ifdef USE_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 +374,7 @@
     if (Py_IS_INFINITY(r))
         errno = ERANGE;
     return r;
+#endif
 }
 
 /*
@@ -364,7 +385,17 @@
 static double
 m_lgamma(double x)
 {
-    double r, absx;
+    double r;
+
+#ifdef USE_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 +433,11 @@
     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 +547,16 @@
     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 +568,7 @@
         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 +576,9 @@
 static double
 m_erfc(double x)
 {
+#ifdef HAVE_ERFC
+    return erfc(x);
+#else
     double absx, cf;
 
     if (Py_IS_NAN(x))
@@ -547,6 +590,7 @@
         cf = m_erfc_contfrac(absx);
         return x > 0.0 ? cf : 2.0 - cf;
     }
+#endif
 }
 
 /*