Re: [Help-gsl] Jacobi polynomials

[Jonny Taylor]

Hi Francesco,

Thanks for your reply. As you noted in your followup email, though, 
unfortunately your suggested method also suffers from cancellation, as I think 
any sum-based method is liable to.

I gave this some thought over the weekend, and have come up with a solution I 
am happy with. I realised that since it is a polynomial it can of course be 
factorized. Having factorized it, rather than a sum over terms we now have a 
product over terms - and therefore there is no problem with cancellation 
between terms. Considered in terms of the coefficients a_m and their associated 
roundoff errors err_m, summation methods have the form:

sum_m( (a_m + err_m) * x^m ) = sum_m( a_m * x^m ) + sum_m( err_m * x^m )

In cases where the answer is much smaller than the constituent a_m we can end 
up with the second (error) sum being of similar or even greater magnitude than 
the first sum.

This contrasts with evaluating a factorized expression:

prod_m( x^m - (b_m + err_m) )

Now (when multiplied out) the error terms will always be much smaller than the 
"signal" term. This method will therefore not suffer from cancellation issues 
when the b_m coefficients are represented as double-precision FP.


The challenge of course is factorizing the polynomial. I believe this *does* 
require extremely high precision intermediate variables (and a lot of 
computation time!). Fortunately, since I am evaluating a given polynomial for 
many, many different x, this slowness is not a problem. "All" I needed to do 
was to do a very tedious find-and-replace in order to make a clone of 
gsl_eigen_nonsymm (and its many dependencies) that uses the GMP 
extended-precision library for all its underlying variables. This is then able 
to find the roots of the polynomial, which can be stored in a lookup and used 
when subsequently evaluating the Jacobi polynomial at any given x.

Best regards,
Jonny