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 _______________________________________________ Help-gsl mailing list [email protected] http://lists.gnu.org/mailman/listinfo/help-gsl
