Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Browse files

MAINT: special: move also orthogonal_eval ufuncs to the new ufunc system

  • Loading branch information...
commit e76714514a06009e99690bd6962e5c821a3abc45 1 parent 401e69d
@pv authored
View
1  .gitattributes
@@ -19,6 +19,5 @@ scipy/spatial/ckdtree.c binary
scipy/spatial/qhull.c binary
scipy/special/_ufuncs.c binary
scipy/special/_ufuncs_cxx.c binary
-scipy/special/orthogonal_eval.c binary
scipy/stats/_rank.c binary
scipy/stats/vonmises_cython.c binary
View
5 scipy/special/add_newdocs.py
@@ -21,6 +21,11 @@ def add_newdoc(place, name, doc):
Internal function, use `lambertw` instead.
""")
+add_newdoc("scipy.special", "_eval_chebyt",
+ """
+ Internal function, use `eval_chebyt` instead.
+ """)
+
add_newdoc("scipy.special", "airy",
"""
(Ai,Aip,Bi,Bip)=airy(z) calculates the Airy functions and their derivatives
View
1  scipy/special/generate_ufuncs.py
@@ -14,6 +14,7 @@
# Ufuncs without C++
UFUNCS = """
_lambertw -- lambertw_scalar: Dld->D -- lambertw.pxd
+_eval_chebyt -- eval_poly_chebyt: ld->d -- orthogonal_eval.pxd
logit -- logitf: f->f, logit: d->d, logitl: g->g -- _logit.h
expit -- expitf: f->f, expit: d->d, expitl: g->g -- _logit.h
bdtrc -- bdtrc: iid->d -- cephes.h
View
6,327 scipy/special/orthogonal_eval.c
0 additions, 6,327 deletions not shown
View
36 scipy/special/orthogonal_eval.pxd
@@ -0,0 +1,36 @@
+"""
+Evaluate orthogonal polynomial values using recurrence relations.
+
+References
+----------
+
+.. [AMS55] Abramowitz & Stegun, Section 22.5.
+
+.. [MH] Mason & Handscombe, Chebyshev Polynomials, CRC Press (2003).
+
+"""
+#
+# Copyright (C) 2009 Pauli Virtanen
+# Distributed under the same license as Scipy.
+#
+
+#------------------------------------------------------------------------------
+# Direct evaluation of polynomials
+#------------------------------------------------------------------------------
+
+from libc.math import sqrt
+
+cdef inline double eval_poly_chebyt(long k, double x) nogil:
@charris
charris added a note

I'm curious as to why the Chebyshev functions of the first kind are computed using Chebyshev functions of the second kind and the identity T_n = (U_n - U_{n-2})/2. I would think it would be less accurate as you can end up differencing large numbers. Am I missing something?

@charris
charris added a note

Nevermind, I got it figured out.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
+ # Use Chebyshev T recurrence directly, see [MH]
+ cdef long m
+ cdef double b2, b1, b0
+
+ b2 = 0
+ b1 = -1
+ b0 = 0
+ x = 2*x
+ for m in range(k+1):
+ b2 = b1
+ b1 = b0
+ b0 = x*b1 - b2
+ return (b0 - b2)/2.0
View
83 scipy/special/orthogonal_eval.pyx → scipy/special/orthogonal_eval.py
@@ -1,6 +1,6 @@
"""
-Evaluate orthogonal polynomial values using recurrence relations
-or by calling special functions.
+Evaluate orthogonal polynomial values using special functions or
+recurrence relations.
References
----------
@@ -16,83 +16,12 @@
#
#------------------------------------------------------------------------------
-# Direct evaluation of polynomials
-#------------------------------------------------------------------------------
-
-cdef extern from "math.h":
- double sqrt(double x) nogil
-
-cdef double eval_poly_chebyt(long k, double x) nogil:
- # Use Chebyshev T recurrence directly, see [MH]
- cdef long m
- cdef double b2, b1, b0
-
- b2 = 0
- b1 = -1
- b0 = 0
- x = 2*x
- for m in range(k+1):
- b2 = b1
- b1 = b0
- b0 = x*b1 - b2
- return (b0 - b2)/2.0
-
-#------------------------------------------------------------------------------
-# Ufunc boilerplate
-#------------------------------------------------------------------------------
-
-cdef extern from "numpy/arrayobject.h":
- void import_array()
- ctypedef int npy_intp
- cdef enum NPY_TYPES:
- NPY_LONG
- NPY_DOUBLE
-
-cdef extern from "numpy/ufuncobject.h":
- void import_ufunc()
- ctypedef void (*PyUFuncGenericFunction)(char**, npy_intp*, npy_intp*, void*)
- object PyUFunc_FromFuncAndData(PyUFuncGenericFunction* func, void** data,
- char* types, int ntypes, int nin, int nout,
- int identity, char* name, char* doc, int c)
-
-cdef void _loop_id_d(char **args, npy_intp *dimensions, npy_intp *steps,
- void *func) nogil:
- cdef int i
- cdef double x
- cdef char *ip1=args[0], *ip2=args[1], *op=args[2]
- for i in range(0, dimensions[0]):
- (<double*>op)[0] = (<double(*)(long,double) nogil>func)(
- (<long*>ip1)[0], (<double*>ip2)[0])
- ip1 += steps[0]; ip2 += steps[1]; op += steps[2]
-
-cdef char _id_d_types[3]
-
-cdef PyUFuncGenericFunction _id_d_funcs[1]
-
-_id_d_types[0] = NPY_LONG
-_id_d_types[1] = NPY_DOUBLE
-_id_d_types[2] = NPY_DOUBLE
-
-_id_d_funcs[0] = _loop_id_d
-
-import_array()
-import_ufunc()
-
-#--
-
-cdef void *chebyt_data[1]
-chebyt_data[0] = <void*>eval_poly_chebyt
-_eval_chebyt = PyUFunc_FromFuncAndData(_id_d_funcs, chebyt_data,
- _id_d_types, 1, 2, 1, 0, "", "", 0)
-
-
-#------------------------------------------------------------------------------
# Actual evaluation functions
#------------------------------------------------------------------------------
import numpy as np
-from scipy.special._ufuncs import gamma, hyp2f1, hyp1f1, gammaln
-from numpy import exp
+from scipy.special._ufuncs import gamma, hyp2f1, hyp1f1, gammaln, _eval_chebyt
+from numpy import exp, sqrt
def binom(n, k):
"""
@@ -100,7 +29,7 @@ def binom(n, k):
Binomial coefficient
"""
- return np.exp(gammaln(1+n) - gammaln(1+k) - gammaln(1+n-k))
+ return exp(gammaln(1+n) - gammaln(1+k) - gammaln(1+n-k))
def eval_jacobi(n, alpha, beta, x, out=None):
"""
@@ -121,7 +50,7 @@ def eval_sh_jacobi(n, p, q, x, out=None):
Evaluate shifted Jacobi polynomial at a point.
"""
- factor = np.exp(gammaln(1+n) + gammaln(n+p) - gammaln(2*n+p))
+ factor = exp(gammaln(1+n) + gammaln(n+p) - gammaln(2*n+p))
return factor * eval_jacobi(n, p-q, q-1, 2*x-1)
def eval_gegenbauer(n, alpha, x, out=None):
View
6 scipy/special/setup.py
@@ -46,12 +46,6 @@ def configuration(parent_package='',top_path=None):
define_macros=[],
libraries=['sc_specfun'])
- # Extension orthogonal_eval
- config.add_extension('orthogonal_eval',
- sources=['orthogonal_eval.c'],
- define_macros=[],
- extra_info=get_info("npymath"))
-
# Extension _ufuncs
curdir = os.path.abspath(os.path.dirname(__file__))
config.add_extension('_ufuncs',
Please sign in to comment.
Something went wrong with that request. Please try again.