From 5f729d28a7040d587fe6a339e362dec9360505ab Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Thu, 26 Mar 2020 17:14:17 +0200 Subject: [PATCH] Avoid using Stead method from GSL above certain threshold. Former-commit-id: aa8012deef69ecce95f203b1b5746cfaa0980806 --- qpms/bessel.c | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/qpms/bessel.c b/qpms/bessel.c index 996f6fc..6c8dc80 100644 --- a/qpms/bessel.c +++ b/qpms/bessel.c @@ -48,13 +48,17 @@ static inline complex double spherical_yn(qpms_l_t l, complex double z) { } #endif +// Don't use Stead algorithm from GSL above this threshold (it fails for x's around 1000000 and higher) +static const double stead_threshold = 1000; + // There is a big issue with gsl's precision of spherical bessel function; these have to be implemented differently qpms_errno_t qpms_sph_bessel_realx_fill(qpms_bessel_t typ, qpms_l_t lmax, double x, complex double *result_array) { int retval; double tmparr[lmax+1]; switch(typ) { case QPMS_BESSEL_REGULAR: - retval = gsl_sf_bessel_jl_steed_array(lmax, x, tmparr); + retval = (x < stead_threshold) ? gsl_sf_bessel_jl_steed_array(lmax, x, tmparr) + : gsl_sf_bessel_jl_array(lmax, x, tmparr); for (int l = 0; l <= lmax; ++l) result_array[l] = tmparr[l]; return retval; break; @@ -70,7 +74,8 @@ qpms_errno_t qpms_sph_bessel_realx_fill(qpms_bessel_t typ, qpms_l_t lmax, double break; case QPMS_HANKEL_PLUS: case QPMS_HANKEL_MINUS: - retval = gsl_sf_bessel_jl_steed_array(lmax, x, tmparr); + retval = (x < stead_threshold) ? gsl_sf_bessel_jl_steed_array(lmax, x, tmparr) + : gsl_sf_bessel_jl_array(lmax, x, tmparr); for (int l = 0; l <= lmax; ++l) result_array[l] = tmparr[l]; if(retval) return retval; retval = gsl_sf_bessel_yl_array(lmax, x, tmparr);