scipy · rgommers · Apr 20, 2024 · Feb 7, 2024 · Feb 15, 2024 · Jan 28, 2024
@@ -4,7 +4,9 @@ import cython
 
 from libc.math cimport log, exp, fabs, sqrt, isnan, isinf, NAN, M_PI
 
-from ._cephes cimport ellpk
+
+cdef extern from "special_wrappers.h" nogil:
+    double cephes_ellpk_wrap(double x)
 
 
 cdef inline double _agm_iter(double a, double b) noexcept nogil:
@@ -64,7 +66,7 @@ cdef inline double agm(double a, double b) noexcept nogil:
 
     if (invsqrthalfmax < a < sqrthalfmax) and (invsqrthalfmax < b < sqrthalfmax):
         e = 4*a*b/(a + b)**2
-        return sgn*(M_PI/4)*(a + b)/ellpk(e)
+        return sgn*(M_PI/4)*(a + b)/cephes_ellpk_wrap(e)
     else:
         # At least one value is "extreme" (very big or very small).
         # Use the iteration to avoid overflow or underflow.

@@ -1,9 +1,10 @@
 from . cimport sf_error
 
 from libc.math cimport NAN, isnan, isinf, isfinite
-cdef extern from "cephes.h" nogil:
-    double ndtr(double a)
-    double ndtri(double y0)
+
+cdef extern from "special_wrappers.h" nogil:
+    double cephes_ndtr_wrap(double a)
+    double cephes_ndtri_wrap(double y0)
 
 cdef extern from "cdflib.h" nogil:
     cdef struct TupleDDI:
@@ -671,7 +672,7 @@ cdef inline double stdtr(double df, double t) noexcept nogil:
     argnames[1] = "df"
 
     if isinf(df) and df > 0:
-        return NAN if isnan(t) else ndtr(t)
+        return NAN if isnan(t) else cephes_ndtr_wrap(t)
 
     if isnan(df) or isnan(t):
       return NAN
@@ -710,7 +711,7 @@ cdef inline double stdtrit(double df, double p) noexcept nogil:
         TupleDID ret
 
     if isinf(df) and df > 0:
-        return NAN if isnan(p) else ndtri(p)
+        return NAN if isnan(p) else cephes_ndtri_wrap(p)
 
     if isnan(p) or isnan(df):
       return NAN

@@ -14,8 +14,8 @@
  *  The ufuncs are used by the class scipy.stats.cosine_gen.
  */
 
+#include "special_wrappers.h"
 #include <math.h>
-#include "cephes/polevl.h"
 
 // M_PI64 is the 64 bit floating point representation of π, e.g.
 //   >>> math.pi.hex()
@@ -101,9 +101,9 @@ double cosine_cdf_pade_approx_at_neg_pi(double x)
     h = (x + M_PI64) + 1.2246467991473532e-16;
     h2 = h*h;
     h3 = h2*h;
-    numer = h3*polevl(h2, numer_coeffs,
+    numer = h3*cephes_polevl_wrap(h2, numer_coeffs,
                       sizeof(numer_coeffs)/sizeof(numer_coeffs[0]) - 1);
-    denom = polevl(h2, denom_coeffs,
+    denom = cephes_polevl_wrap(h2, denom_coeffs,
                    sizeof(denom_coeffs)/sizeof(denom_coeffs[0]) - 1);
     return numer / denom;
 }
@@ -178,7 +178,7 @@ double _p2(double t)
                        0.5};
     double v;
 
-    v = polevl(t, coeffs, sizeof(coeffs) / sizeof(coeffs[0]) - 1);
+    v = cephes_polevl_wrap(t, coeffs, sizeof(coeffs) / sizeof(coeffs[0]) - 1);
     return v;
 }
 
@@ -193,7 +193,7 @@ double _q2(double t)
                        1.0};
     double v;
 
-    v = polevl(t, coeffs, sizeof(coeffs) / sizeof(coeffs[0]) - 1);
+    v = cephes_polevl_wrap(t, coeffs, sizeof(coeffs) / sizeof(coeffs[0]) - 1);
     return v;
 }
 
@@ -225,7 +225,7 @@ double _poly_approx(double s)
     // Here we include terms up to s**13.
     //
     s2 = s*s;
-    p = s*polevl(s2, coeffs, sizeof(coeffs)/sizeof(coeffs[0]) - 1);
+    p = s*cephes_polevl_wrap(s2, coeffs, sizeof(coeffs)/sizeof(coeffs[0]) - 1);
     return p;
 }
 

@@ -10,17 +10,20 @@ cdef extern from "_complexstuff.h":
     np.npy_cdouble npy_cexp(np.npy_cdouble x) nogil
 
 
-cdef extern from "cephes/dd_real.h":
+cdef extern from "dd_real_wrappers.h":
     ctypedef struct double2:
-        double x[2]
+        double hi
+        double lo
 
     double2 dd_create_d(double x) nogil
-    double2 dd_add(const double2 a, const double2 b) nogil
-    double2 dd_mul(const double2 a, const double2 b) nogil
-    double dd_to_double(const double2 a) nogil
-
-from ._cephes cimport log1p, expm1, cosm1
+    double2 dd_add(const double2* a, const double2* b) nogil
+    double2 dd_mul(const double2* a, const double2* b) nogil
+    double dd_to_double(const double2* a) nogil
 
+cdef extern from "special_wrappers.h" nogil:
+    double cephes_cosm1_wrap(double x)
+    double cephes_expm1_wrap(double x)
+    double cephes_log1p_wrap(double x)
 
 # log(z + 1) = log(x + 1 + 1j*y)
 #             = log(sqrt((x+1)**2 + y**2)) + 1j*atan2(y, x+1)
@@ -49,15 +52,15 @@ cdef inline double complex clog1p(double complex z) noexcept nogil:
     zi = z.imag
 
     if zi == 0.0 and zr >= -1.0:
-        return zpack(log1p(zr), 0.0)
+        return zpack(cephes_log1p_wrap(zr), 0.0)
 
     az = zabs(z)
     if az < 0.707:
         azi = fabs(zi)
         if zr < 0 and fabs(-zr - azi*azi/2)/(-zr) < 0.5:
             return clog1p_ddouble(zr, zi)
         else:
-            x = 0.5 * log1p(az*(az + 2*zr/az))
+            x = 0.5 * cephes_log1p_wrap(az*(az + 2*zr/az))
             y = atan2(zi, zr + 1.0)
             return zpack(x, y)
 
@@ -73,13 +76,13 @@ cdef inline double complex clog1p_ddouble(double zr, double zi) noexcept nogil:
     i = dd_create_d(zi)
     two = dd_create_d(2.0)
 
-    rsqr = dd_mul(r, r)
-    isqr = dd_mul(i, i)
-    rtwo = dd_mul(two, r)
-    absm1 = dd_add(rsqr, isqr)
-    absm1 = dd_add(absm1, rtwo)
+    rsqr = dd_mul(&r,& r)
+    isqr = dd_mul(&i, &i)
+    rtwo = dd_mul(&two, &r)
    absm1 = dd_add(&rsqr, &isqr)
+    absm1 = dd_add(&absm1, &rtwo)
 
-    x = 0.5 * log1p(dd_to_double(absm1))
+    x = 0.5 * cephes_log1p_wrap(dd_to_double(&absm1))
     y = atan2(zi, zr+1.0)
     return zpack(x, y)
 
@@ -104,8 +107,8 @@ cdef inline double complex cexpm1(double complex z) noexcept nogil:
     if zr <= -40:
         x = -1.0
     else:
-        ezr = expm1(zr)
-        x = ezr*cos(zi) + cosm1(zi)
+        ezr = cephes_expm1_wrap(zr)
+        x = ezr*cos(zi) + cephes_cosm1_wrap(zi)
     # don't compute exp(zr) too, unless necessary
     if zr > -1.0:
         y = (ezr + 1.0)*sin(zi)

@@ -1,5 +1,6 @@
-from ._cephes cimport ellpk
+cdef extern from "special_wrappers.h" nogil:
+    double cephes_ellpk_wrap(double x)
 
 
 cdef inline double ellipk(double m) noexcept nogil:
-    return ellpk(1.0 - m)
+    return cephes_ellpk_wrap(1.0 - m)
@@ -1,8 +1,9 @@
-from ._cephes cimport Gamma
+cdef extern from "special_wrappers.h" nogil:
+    double cephes_gamma_wrap(double x)
 
 
 cdef inline double _factorial(double n) noexcept nogil:
     if n < 0:
         return 0
     else:
-        return Gamma(n + 1)
+        return cephes_gamma_wrap(n + 1)
@@ -1,12 +1,18 @@
 from libc.math cimport pow, sqrt, floor, log, log1p, exp, M_PI, NAN, fabs, isinf
 cimport numpy as np
 
-from ._cephes cimport iv, jv, Gamma, lgam, gammasgn
 from ._xlogy cimport xlogy
 from ._complexstuff cimport (
     zsqrt, zpow, zabs, npy_cdouble_from_double_complex,
     double_complex_from_npy_cdouble)
 
+cdef extern from "special_wrappers.h" nogil:
+    double cephes_iv_wrap(double v, double x)
+    double cephes_jv_wrap(double v, double x)
+    double cephes_gamma_wrap(double x)
+    double cephes_lgam_wrap(double x)
+    double cephes_gammasgn_wrap(double x)
+
 cdef extern from "float.h":
     double DBL_MAX, DBL_MIN
 
@@ -37,17 +43,17 @@ cdef inline double _hyp0f1_real(double v, double z) noexcept nogil:
 
     if z > 0:
         arg = sqrt(z)
-        arg_exp = xlogy(1.0-v, arg) + lgam(v)
-        bess_val = iv(v-1, 2.0*arg)
+        arg_exp = xlogy(1.0-v, arg) + cephes_lgam_wrap(v)
+        bess_val = cephes_iv_wrap(v-1, 2.0*arg)
 
         if (arg_exp > log(DBL_MAX) or bess_val == 0 or   # overflow
             arg_exp < log(DBL_MIN) or isinf(bess_val)):  # underflow
             return _hyp0f1_asy(v, z)
         else:
-            return exp(arg_exp) * gammasgn(v) * bess_val
+            return exp(arg_exp) * cephes_gammasgn_wrap(v) * bess_val
     else:
         arg = sqrt(-z)
-        return pow(arg, 1.0 - v) * Gamma(v) * jv(v - 1, 2*arg)
+        return pow(arg, 1.0 - v) * cephes_gamma_wrap(v) * cephes_jv_wrap(v - 1, 2*arg)
 
 
 cdef inline double _hyp0f1_asy(double v, double z) noexcept nogil:
@@ -68,8 +74,8 @@ cdef inline double _hyp0f1_asy(double v, double z) noexcept nogil:
 
     arg_exp_i = -0.5*log(p1)
     arg_exp_i -= 0.5*log(2.0*M_PI*v1)
-    arg_exp_i += lgam(v)
-    gs = gammasgn(v)
+    arg_exp_i += cephes_lgam_wrap(v)
+    gs = cephes_gammasgn_wrap(v)
 
     arg_exp_k = arg_exp_i
     arg_exp_i += v1 * eta
@@ -127,4 +133,4 @@ cdef inline double complex _hyp0f1_cmplx(double v, double complex z) noexcept no
         s = 2.0 * arg
         r = special_ccyl_bessel_j(v-1.0, npy_cdouble_from_double_complex(s))
 
-    return double_complex_from_npy_cdouble(r) * Gamma(v) * zpow(arg, 1.0 - v)
+    return double_complex_from_npy_cdouble(r) * cephes_gamma_wrap(v) * zpow(arg, 1.0 - v)
@@ -3,10 +3,10 @@ from libc.math cimport fabs, exp, floor, isnan, M_PI, NAN, INFINITY
 import cython
 
 from . cimport sf_error
-from ._cephes cimport expm1, poch
 
 cdef extern from 'special_wrappers.h':
     double hypU_wrap(double, double, double) nogil
+    double cephes_poch_wrap(double x, double m) nogil
 
 
 @cython.cdivision(True)
@@ -25,6 +25,6 @@ cdef inline double hyperu(double a, double b, double x) noexcept nogil:
             return INFINITY
         else:
             # DLMF 13.2.14-15 and 13.2.19-21
-            return poch(1.0 - b + a, -a)
+            return cephes_poch_wrap(1.0 - b + a, -a)
 
     return hypU_wrap(a, b, x)