[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
## 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,
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

[Prev in Thread] |
**Current Thread** |
[Next in Thread] |

**Re: [Help-gsl] Jacobi polynomials**,
*Jonny Taylor* **<=**