Merge pull request scipy#8409 from person142/cephes-pxd

MAINT: special: add a `.pxd` file for Cephes functions

scipy/special/_agm.pxd

@@ -7,9 +7,7 @@ cdef extern from "_c99compat.h":
     int sc_isnan(double x) nogil
     int sc_isinf(double x) nogil
 
-cdef extern from "cephes.h":
-    double ellpk(double m) nogil
-
+from ._cephes cimport ellpk
 from ._complexstuff cimport pi, nan
 
 
@@ -55,7 +53,7 @@ cdef inline double agm(double a, double b) nogil:
 
     if a == 0 or b == 0:
         # At least one of the arguments is 0.
-        return 0.0 
+        return 0.0
 
     if a == b:
         return a

scipy/special/_cephes.pxd

@@ -0,0 +1,98 @@
+cdef extern from "cephes.h" nogil:
+    int airy(double x, double *ai, double *aip, double *bi, double *bip);
+    double bdtrc(int k, int n, double p);
+    double bdtr(int k, int n, double p);
+    double bdtri(int k, int n, double y);
+    double beta(double a, double b);
+    double lbeta(double a, double b);
+    double btdtr(double a, double b, double x);
+    double cbrt(double x);
+    double chdtrc(double df, double x);
+    double chdtr(double df, double x);
+    double chdtri(double df, double y);
+    double dawsn(double xx);
+    double ellie(double phi, double m);
+    double ellik(double phi, double m);
+    double ellpe(double x);
+    int ellpj(double u, double m, double *sn, double *cn, double *dn, double *ph);
+    double ellpk(double x);
+    double exp10(double x);
+    double exp1m(double x);
+    double exp2(double x);
+    double expn(int n, double x);
+    double fdtrc(double a, double b, double x);
+    double fdtr(double a, double b, double x);
+    double fdtri(double a, double b, double y);
+    int fresnl(double xxa, double *ssa, double *cca);
+    double Gamma(double x);
+    double lgam(double x);
+    double gdtr(double a, double b, double x);
+    double gdtrc(double a, double b, double x);
+    double gdtri(double a, double b, double y);
+    double hyp2f1(double a, double b, double c, double x);
+    double hyperg(double a, double b, double x);
+    double hyp2f0(double a, double b, double x, int type, double *err);
+    double onef2(double a, double b, double c, double x, double *err);
+    double threef0(double a, double b, double c, double x, double *err);
+    double i0(double x);
+    double i0e(double x);
+    double i1(double x);
+    double i1e(double x);
+    double igamc(double a, double x);
+    double igam(double a, double x);
+    double igam_fac( double a, double x);
+    double igamci( double a, double q);
+    double igami(double a, double p);
+    double incbet(double aa, double bb, double xx);
+    double incbi(double aa, double bb, double yy0);
+    double iv(double v, double x);
+    double j0(double x);
+    double y0(double x);
+    double j1(double x);
+    double y1(double x);
+    double jn(int n, double x);
+    double jv(double n, double x);
+    double k0(double x);
+    double k0e(double x);
+    double k1(double x);
+    double k1e(double x);
+    double kn(int nn, double x);
+    int levnsn(int n, double r[], double a[], double e[], double refl[]);
+    double nbdtrc(int k, int n, double p);
+    double nbdtr(int k, int n, double p);
+    double nbdtri(int k, int n, double p);
+    double ndtr(double a);
+    double log_ndtr(double a);
+    double erfc(double a);
+    double erf(double x);
+    double ndtri(double y0);
+    double pdtrc(int k, double m);
+    double pdtr(int k, double m);
+    double pdtri(int k, double y);
+    double psi(double x);
+    double rgamma(double x);
+    double round(double x);
+    int shichi(double x, double *si, double *ci);
+    int sici(double x, double *si, double *ci);
+    double radian(double d, double m, double s);
+    double sindg(double x);
+    double cosdg(double x);
+    double spence(double x);
+    double stdtr(int k, double t);
+    double stdtri(int k, double p);
+    double yv(double v, double x);
+    double tandg(double x);
+    double cotdg(double x);
+    double log1p(double x);
+    double log1pmx( double x);
+    double expm1(double x);
+    double cosm1(double x);
+    double lgam1p(double x);
+    double yn(int n, double x);
+    double zeta(double x, double q);
+    double zetac(double x);
+    double smirnov (int n, double e);
+    double smirnovi (int n, double p);
+    double kolmogorov(double x);
+    double kolmogi(double p);
+    double lanczos_sum_expg_scaled(double x);

scipy/special/_cunity.pxd

@@ -8,7 +8,7 @@ cdef extern from "_complexstuff.h":
     double npy_atan2(double y, double x) nogil
     np.npy_cdouble npy_clog(np.npy_cdouble x) nogil
     np.npy_cdouble npy_cexp(np.npy_cdouble x) nogil
-    
+
 
 cdef extern from "c_misc/double2.h":
     ctypedef struct double2_t:
@@ -19,22 +19,20 @@ cdef extern from "c_misc/double2.h":
     void double2_mul(double2_t* a, double2_t* b, double2_t* c) nogil
     double double2_double(double2_t* a) nogil
 
-cdef extern from "cephes.h":
-    double log1p(double x) nogil
-    double expm1(double x) nogil
-    double cosm1(double x) nogil
+from ._cephes cimport log1p, expm1, cosm1
+
 
 # log(z + 1) = log(x + 1 + 1j*y)
 #             = log(sqrt((x+1)**2 + y**2)) + 1j*atan2(y, x+1)
-# 
+#
 # Using atan2(y, x+1) for the imaginary part is always okay.  The real part
 # needs to be calculated more carefully.  For |z| large, the naive formula
-# log(z + 1) can be used.  When |z| is small, rewrite as 
+# log(z + 1) can be used.  When |z| is small, rewrite as
 #
 # log(sqrt((x+1)**2 + y**2)) = 0.5*log(x**2 + 2*x +1 + y**2)
 #       = 0.5 * log1p(x**2 + y**2 + 2*x)
 #       = 0.5 * log1p(hypot(x,y) * (hypot(x, y) + 2*x/hypot(x,y)))
-# 
+#
 # This expression suffers from cancellation when x < 0 and
 # y = +/-sqrt(2*fabs(x)). To get around this cancellation problem, we use
 # double-double precision when necessary.
@@ -43,15 +41,15 @@ cdef inline double complex clog1p(double complex z) nogil:
     cdef np.npy_cdouble ret
 
     if not zisfinite(z):
-        z = z + 1 
+        z = z + 1
         ret = npy_clog(npy_cdouble_from_double_complex(z))
         return double_complex_from_npy_cdouble(ret)
 
     zr = z.real
     zi = z.imag
 
     if zi == 0.0 and zr >= -1.0:
-        return zpack(log1p(zr), 0.0) 
+        return zpack(log1p(zr), 0.0)
 
     az = zabs(z)
     if az < 0.707:
@@ -74,7 +72,7 @@ cdef inline double complex clog1p_ddouble(double zr, double zi) nogil:
     double2_init(&r, zr)
     double2_init(&i, zi)
     double2_init(&two, 2.0)
-    
+
     double2_mul(&r, &r, &rsqr)
     double2_mul(&i, &i, &isqr)
     double2_mul(&two, &r, &rtwo)
@@ -83,10 +81,10 @@ cdef inline double complex clog1p_ddouble(double zr, double zi) nogil:
 
     x = 0.5 * log1p(double2_double(&absm1))
     y = npy_atan2(zi, zr+1.0)
-    return zpack(x, y) 
- 
+    return zpack(x, y)
+
 # cexpm1(z) = cexp(z) - 1
-# 
+#
 # The imaginary part of this is easily computed via exp(z.real)*sin(z.imag)
 # The real part is difficult to compute when there is cancellation e.g. when
 # z.real = -log(cos(z.imag)).  There isn't a way around this problem  that
@@ -115,4 +113,3 @@ cdef inline double complex cexpm1(double complex z) nogil:
         y = exp(zr)*sin(zi)
 
     return zpack(x, y)
-

scipy/special/_digamma.pxd

@@ -14,12 +14,9 @@ import cython
 from libc.math cimport ceil, fabs, M_PI
 from ._complexstuff cimport number_t, nan, zlog, zabs
 from ._trig cimport sinpi, cospi
+from ._cephes cimport zeta, psi
 from . cimport sf_error
 
-cdef extern from "cephes.h":
-    double zeta(double x, double q) nogil
-    double psi(double x) nogil
-
 # Use the asymptotic series for z away from the negative real axis
 # with abs(z) > smallabsz.
 DEF smallabsz = 16

scipy/special/_hyp0f1.pxd

@@ -10,11 +10,7 @@ from ._complexstuff cimport (
 cdef extern from "float.h":
     double DBL_MAX, DBL_MIN
 
-cdef extern from "cephes.h":
-    double iv(double v, double x) nogil
-    double jv(double n, double x) nogil
-    double Gamma(double x) nogil
-    double lgam(double x) nogil
+from ._cephes cimport iv, jv, Gamma, lgam
 
 cdef extern from "c_misc/misc.h":
     double gammasgn(double x) nogil
@@ -30,7 +26,7 @@ cdef extern from "amos_wrappers.h":
 cdef inline double _hyp0f1_real(double v, double z) nogil:
     cdef double arg, v1, arg_exp, bess_val
 
-    # handle poles, zeros 
+    # handle poles, zeros
     if v <= 0.0 and v == floor(v):
         return NAN
     if z == 0.0 and v != 0.0:
@@ -44,7 +40,7 @@ cdef inline double _hyp0f1_real(double v, double z) nogil:
         arg = sqrt(z)
         arg_exp = xlogy(1.0-v, arg) + lgam(v)
         bess_val = iv(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)
@@ -56,7 +52,7 @@ cdef inline double _hyp0f1_real(double v, double z) nogil:
 
 
 cdef inline double _hyp0f1_asy(double v, double z) nogil:
-    r"""Asymptotic expansion for I_{v-1}(2*sqrt(z)) * Gamma(v) 
+    r"""Asymptotic expansion for I_{v-1}(2*sqrt(z)) * Gamma(v)
     for real $z > 0$ and $v\to +\infty$.
 
     Based off DLMF 10.41
@@ -109,7 +105,7 @@ cdef inline double complex _hyp0f1_cmplx(double v, double complex z) nogil:
         double complex arg, s
         double complex t1, t2
 
-    # handle poles, zeros 
+    # handle poles, zeros
     if v <= 0.0 and v == floor(v):
         return NAN
     if z.real == 0.0 and z.imag == 0.0 and v != 0.0:
@@ -133,4 +129,3 @@ cdef inline double complex _hyp0f1_cmplx(double v, double complex z) nogil:
         r = cbesj_wrap(v-1.0, npy_cdouble_from_double_complex(s))
 
     return double_complex_from_npy_cdouble(r) * Gamma(v) * zpow(arg, 1.0 - v)
-

scipy/special/_legacy.pxd

@@ -17,22 +17,9 @@ cdef extern from "numpy/npy_math.h" nogil:
     double npy_isnan(double)
     double nan "NPY_NAN"
 
-cdef extern from "cephes.h":
-    double bdtrc(int k, int n, double p) nogil
-    double bdtr(int k, int n, double p) nogil
-    double bdtri(int k, int n, double y) nogil
-    double expn(int n, double x) nogil
-    double hyp2f0(double a, double b, double x, int type, double *err) nogil
-    double nbdtrc(int k, int n, double p) nogil
-    double nbdtr(int k, int n, double p) nogil
-    double nbdtri(int k, int n, double p) nogil
-    double pdtrc(int k, double m) nogil
-    double pdtr(int k, double m) nogil
-    double pdtri(int k, double y) nogil
-    double kn(int n, double x) nogil
-    double yn(int n, double x) nogil
-    double smirnov(int n, double e) nogil
-    double smirnovi(int n, double p) nogil
+from ._cephes cimport (bdtrc, bdtr, bdtri, expn, hyp2f0, nbdtrc,
+                       nbdtr, nbdtri, pdtrc, pdtr, pdtri, kn, yn,
+                       smirnov, smirnovi)
 
 cdef extern from "amos_wrappers.h":
     double cbesk_wrap_real_int(int n, double z) nogil

scipy/special/_spherical_bessel.pxd

@@ -39,9 +39,7 @@ cdef extern from "amos_wrappers.h":
     npy_cdouble cbesy_wrap(double v, npy_cdouble z) nogil
     double cbesy_wrap_real(double v, double x) nogil
 
-cdef extern from "cephes.h":
-    double iv(double v, double x) nogil
-
+from ._cephes cimport iv
 
 # Fused type wrappers
 

scipy/special/orthogonal_eval.pxd

@@ -11,7 +11,7 @@ References
 
 .. [LP] P. Levrie & R. Piessens, A note on the evaluation of orthogonal
         polynomials using recurrence relations, Internal Report TW74 (1985)
-        Dept. of Computer Science, K.U. Leuven, Belgium 
+        Dept. of Computer Science, K.U. Leuven, Belgium
         https://lirias.kuleuven.be/handle/123456789/131600
 
 """
@@ -31,17 +31,12 @@ from ._complexstuff cimport (
     double_complex_from_npy_cdouble)
 
 from . cimport sf_error
-
-cdef extern from "cephes.h":
-    double Gamma(double x) nogil
-    double lgam(double x) nogil
-    double beta (double a, double b) nogil
-    double lbeta (double a, double b) nogil
-    double hyp2f1_wrap "hyp2f1" (double a, double b, double c, double x) nogil 
+from ._cephes cimport Gamma, lgam, beta, lbeta
+from ._cephes cimport hyp2f1 as hyp2f1_wrap
 
 cdef extern from "specfun_wrappers.h":
     double hyp1f1_wrap(double a, double b, double x) nogil
-    npy_cdouble chyp2f1_wrap( double a, double b, double c, npy_cdouble z) nogil 
+    npy_cdouble chyp2f1_wrap( double a, double b, double c, npy_cdouble z) nogil
     npy_cdouble chyp1f1_wrap( double a, double b, npy_cdouble z) nogil
 
 cdef extern from "c_misc/misc.h":
@@ -135,9 +130,9 @@ cdef inline double binom(double n, double k) nogil:
 #-----------------------------------------------------------------------------
 
 cdef inline number_t eval_jacobi(double n, double alpha, double beta, number_t x) nogil:
-    cdef double a, b, c, d 
+    cdef double a, b, c, d
     cdef number_t g
-    
+
     d = binom(n+alpha, n)
     a = -n
     b = n + alpha + beta + 1
@@ -156,10 +151,10 @@ cdef inline double eval_jacobi_l(long n, double alpha, double beta, double x) no
     elif n == 0:
         return 1.0
     elif n == 1:
-        return 0.5*(2*(alpha+1)+(alpha+beta+2)*(x-1)) 
+        return 0.5*(2*(alpha+1)+(alpha+beta+2)*(x-1))
     else:
         d = (alpha+beta+2)*(x - 1) / (2*(alpha+1))
-        p = d + 1 
+        p = d + 1
         for kk in range(n-1):
             k = kk+1.0
             t = 2*k+alpha+beta
@@ -233,7 +228,7 @@ cdef inline double eval_gegenbauer_l(long n, double alpha, double x) nogil:
         return p
     else:
         d = x - 1
-        p = x 
+        p = x
         for kk in range(n-1):
             k = kk+1.0
             d = (2*(k+alpha)/(k+2*alpha))*(x-1)*p + (k/(k+2*alpha)) * d
@@ -409,7 +404,7 @@ cdef inline double eval_legendre_l(long n, double x) nogil:
         return p
     else:
         d = x - 1
-        p = x 
+        p = x
         for kk in range(n-1):
             k = kk+1.0
             d = ((2*k+1)/(k+1))*(x-1)*p + (k/(k+1)) * d
@@ -463,8 +458,8 @@ cdef inline double eval_genlaguerre_l(long n, double alpha, double x) nogil:
     elif n == 1:
         return -x+alpha+1
     else:
-        d = -x/(alpha+1) 
-        p = d + 1 
+        d = -x/(alpha+1)
+        p = d + 1
         for kk in range(n-1):
             k = kk+1.0
             d = -x/(k+alpha+1)*p + (k/(k+alpha+1)) * d
@@ -511,4 +506,3 @@ cdef inline double eval_hermitenorm(long n, double x) nogil:
 @cython.cdivision(True)
 cdef inline double eval_hermite(long n, double x) nogil:
     return eval_hermitenorm(n, sqrt(2)*x) * 2**(n/2.0)
-