[Help-gsl] Jacobi polynomials

[Jonny Taylor]

Hi all,

I hope this is not considered too off-topic: my question refers to a 
third-party extension to GSL (for Jacobi polynomials), but I am hoping the 
mathematical experts here may be able to point me in the right direction.

I have been using the contributed extension found here:
http://www.network-theory.co.uk/download/gslextras/Jacobi/

This works ok up to a point, but the recurrence relation suffers from 
catastrophic cancelation under some circumstances, such as x->1 for large n and 
negative a. I attempted to enhance the code by implementing an alternative 
method based on series summation. The original reference for this method is as 
follows (though I have not yet been able to acquire this primary reference):
G. Szego, Orthogonal Polynomials, Amer. Math. Soc. Colloq. Publ., Vol. 23, 
Amer. Math. Soc., p62

Unfortunately, there are some circumstances when BOTH methods suffer from 
cancelation e.g. P_34^(-16,16) for x~0.8055. I am therefore rather stuck as to 
how I can obtain an accurate answer for cases such as this. I notice that 
Mathematica is able to give an accurate answer in such cases, but I have no 
idea how to go about calculating that answer myself (I am mathematically out of 
my depth with this sort of thing). Should I be looking for yet another way of 
generating the Jacobi function, or do I simply need to resort to higher 
precision arithmetic in these pathological cases?

Thanks in advance for any suggestions!
Jonny