help-gsl
[Top][All Lists]

## Re: [Help-gsl] Jacobi polynomials

 From: Jonny Taylor Subject: Re: [Help-gsl] Jacobi polynomials Date: Wed, 2 Jun 2010 15:05:18 +0100

```Hi Francesco,

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

```