From 773866fe3502d3748bd224f2af0dd3e53f932a92 Mon Sep 17 00:00:00 2001 From: Josh Wilson Date: Sun, 11 Feb 2018 22:13:02 -0600 Subject: [PATCH] MAINT: special: add a `.pxd` file for Cephes functions --- scipy/special/_agm.pxd | 6 +- scipy/special/_cephes.pxd | 98 +++++++++++++++++++++++++++++ scipy/special/_cunity.pxd | 27 ++++---- scipy/special/_digamma.pxd | 5 +- scipy/special/_hyp0f1.pxd | 15 ++--- scipy/special/_legacy.pxd | 19 +----- scipy/special/_spherical_bessel.pxd | 4 +- scipy/special/orthogonal_eval.pxd | 30 ++++----- 8 files changed, 134 insertions(+), 70 deletions(-) create mode 100644 scipy/special/_cephes.pxd diff --git a/scipy/special/_agm.pxd b/scipy/special/_agm.pxd index b6910df366a3..0c3eb9890778 100644 --- a/scipy/special/_agm.pxd +++ b/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 diff --git a/scipy/special/_cephes.pxd b/scipy/special/_cephes.pxd new file mode 100644 index 000000000000..857073b99b81 --- /dev/null +++ b/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); diff --git a/scipy/special/_cunity.pxd b/scipy/special/_cunity.pxd index 15763ee6e578..59566ede50cf 100644 --- a/scipy/special/_cunity.pxd +++ b/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,7 +41,7 @@ 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) @@ -51,7 +49,7 @@ cdef inline double complex clog1p(double complex z) nogil: 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) - diff --git a/scipy/special/_digamma.pxd b/scipy/special/_digamma.pxd index 1fdf6d1b2faf..2880eccc824c 100644 --- a/scipy/special/_digamma.pxd +++ b/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 diff --git a/scipy/special/_hyp0f1.pxd b/scipy/special/_hyp0f1.pxd index 2a008f6f2992..3a4967f36df9 100644 --- a/scipy/special/_hyp0f1.pxd +++ b/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) - diff --git a/scipy/special/_legacy.pxd b/scipy/special/_legacy.pxd index 70287ef01952..a55e85793d46 100644 --- a/scipy/special/_legacy.pxd +++ b/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 diff --git a/scipy/special/_spherical_bessel.pxd b/scipy/special/_spherical_bessel.pxd index dbf0a97b8c4d..c05abe0b431c 100644 --- a/scipy/special/_spherical_bessel.pxd +++ b/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 diff --git a/scipy/special/orthogonal_eval.pxd b/scipy/special/orthogonal_eval.pxd index 57821f8ae785..6f8e9258401a 100644 --- a/scipy/special/orthogonal_eval.pxd +++ b/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) -