[SciPy-user] Orthogonal polynomial evaluation

[Andrew York]

 I am running a simulation with Scipy that needs to numerically
evaluate Jacobi polynomials. I have been using the function
scipy.special.orthogonal.jacobi() heavily. Evaluations of this
function are >90% of my simulation running time. Much more troubling,
I begin to see unphysical results when I need to evaluate high-order
Jacobi polynomials, and I suspect I am running into the problem warned
against here:

http://docs.scipy.org/doc/scipy/reference/special.html

"Warning Evaluating large-order polynomials using these functions can
be numerically unstable."

Another Scipy user seems to have run into a similar problem:

http://thread.gmane.org/gmane.comp.python.scientific.user/11206

 So, it seems the correct course of action is clear: I should identify
a stable algorithm to numerically evaluate Jacobi polynomials, I
should implement this algorithm in something fast like C, and then
call it from within python. I found some promising information on a
possible algorithm:

http://portal.acm.org/citation.cfm?id=363393
http://portal.acm.org/citation.cfm?id=362686.362700
https://people.sc.fsu.edu/~burkardt/f77_src/toms332/toms332.html

But I do not speak FORTRAN. I am also new to python and not yet
familiar with how to integrate a fast low-level language with python,
nor am I especially proficient in any low-level language. I suspect if
numerical evaluation of these functions was faster or more robust,
other scipy users would also benefit. Should I attempt this myself (in
C with weave?)? If I do, would others be interested in the results?

Is this an easy task for anyone else? I would greatly appreciate any
help or guidance.

Thank you for your attention.