<table cellspacing="0" cellpadding="0" border="0" ><tr><td valign="top" style="font: inherit;">I recall that stable routines to calculate Bessel functions use descending recurrence relations (descending in indices), rather than ascending recurrence which is unstable. Only trick is that the values have to be renormalized at the end (or as you go to keep the numbers in floating point range, I think) since descending recurrence does not set the scale initially. I wonder if this is the case for other recurrence relations. That is, if one recurrence is unstable (e.g. ascending), then the opposite (descending) will be stable. Perhaps the scipy routines can be easily reworked, if so. Just a thought.<br><br>-- Lou Pecora, my views are my own.<br><br>--- On <b>Mon, 4/13/09, Anne Archibald <i><firstname.lastname@example.org></i></b> wrote:<br><blockquote style="border-left: 2px solid rgb(16, 16, 255); margin-left: 5px; padding-left: 5px;"><br>From: Anne
Archibald <email@example.com><br><div class="plainMail"><br>Yes, unfortunately all scipy's orthogonal polynomials are computed<br>using generic recurrence relations. This allows a quick flexible<br>implementation, but as you and I both discovered it suffers from<br>numerical instability.<br><br>> Another Scipy user seems to have run into a similar problem:<br>><br>> <a href="http://thread.gmane.org/gmane.comp.python.scientific.user/11206" target="_blank">http://thread.gmane.org/gmane.comp.python.scientific.user/11206</a><br>><br>> So, it seems the correct course of action is clear: I should identify<br>> a stable algorithm to numerically evaluate Jacobi polynomials, I<br>> should implement this algorithm in something fast like C, and then<br>> call it from within python. I found some promising information on a<br>> possible algorithm:<br>><br>> <a href="http://portal.acm.org/citation.cfm?id=363393"
target="_blank">http://portal.acm.org/citation.cfm?id=363393</a><br>> <a href="http://portal.acm.org/citation.cfm?id=362686.362700" target="_blank">http://portal.acm.org/citation.cfm?id=362686.362700</a><br>> <a href="https://people.sc.fsu.edu/~burkardt/f77_src/toms332/toms332.html" target="_blank">https://people.sc.fsu.edu/~burkardt/f77_src/toms332/toms332.html</a><br>><br>> But I do not speak FORTRAN. I am also new to python and not yet<br>> familiar with how to integrate a fast low-level language with python,<br>> nor am I especially proficient in any low-level language. I suspect if<br>> numerical evaluation of these functions was faster or more robust,<br>> other scipy users would also benefit. Should I attempt this myself (in<br>> C with weave?)? If I do, would others be interested in the results?<br>><br>> Is this an easy task for anyone else? I would greatly appreciate any<br>> help or guidance.<br><br>If you are comfortable with python but not
FORTRAN or<br>python-to-other-language linking, then I recommend using cython to<br>implement compiled extensions. It can be very fast, and the path to a<br>working implementation is all modest steps:<br>* write a python point-at-a-time implementation based on the articles<br>* write a python unit testing module (makes the rest much easier!)<br>* compile the point-at-a-time implementation with cython<br>* add enough type annotations that you can see that the generated C<br>code contains no python calls<br>* write a wrapper function in cython that calls the point-at-a-time<br>function on each element of an array<br><br>Ideally, I'd like to see the special functions module modified so that<br>the orthogonal polynomials could use better implementations where<br>available, falling back to the recurrence relation. But it's kind of a<br>maze.<br><br>Good luck,<br>Anne<br>_______________________________________________<br>SciPy-user mailing list<br><a
ymailto="mailto:SciPyfirstname.lastname@example.org" href="/mc/compose?to=SciPyemail@example.com">SciPyfirstname.lastname@example.org</a><br><a href="http://mail.scipy.org/mailman/listinfo/scipy-user" target="_blank">http://mail.scipy.org/mailman/listinfo/scipy-user</a><br></div></blockquote></td></tr></table><br>