[Bug-gsl] Inaccuracy of the Airy function due to invocation of GSL's cosine function with large input parameters

[address@hidden]

Dear GSL developers,



This is a follow-up email of my previous report of regarding GSL’s Airy and
cosine functions (https://lists.gnu.org/archive/html/bug-gsl/2017-11/
msg00011.html).


Here I would like to report another kind of unexpected results when using
the Airy function gsl_sf_Ai_e.  Code to reproduce the issue:


#include <gsl/gsl_math.h>

#include <gsl/gsl_errno.h>

#include <gsl/gsl_sf_trig.h>

#include <gsl/gsl_sf_airy.h>



int main(){

  gsl_sf_result res;

  double   x = -1.14e34;

  gsl_mode_t mode = GSL_PREC_DOUBLE;

  int status = gsl_sf_airy_Ai_e(x, mode, &res);

  printf("status = %d, val = %g, err = %g\n", status,res.val,res.err);

}



Running the code gives:  status = 0, val = -inf, err = inf



The calculation result here, val = -inf, is incorrect: Airy functions are
damped oscillatory for negative x; using Mathematica you get AiryAi
[-1.1e34]=-1.36e-9.  Besides, the error estimate, inf, is obviously
imprecise, and the returned status, GSL_SUCCESS 0, may be somewhat
misleading.



Cause of the issue in my humble opinion: Airy_ai (-1.1e34) invokes
gsl_sf_cos_err_e (theta.val, theta.err, &cos_result) where theta.val =
-8.1145794715437919e+50, theta.err = 7.4985953321244595e+35 (according to
gdb). After this invocation, cos_result.val and cos_result.err become -inf
and inf respectively, which is clearly wrong, which stays away from its
[-1,1] bound.  Since “it is known that the GSL gsl_sf_cos and gsl_sf_sin
functions fail for large inputs” (from https://lists.gnu.org/archive/
html/bug-gsl/2017-11/msg00012.html), maybe GSL's Airy functions, both the
first and the second kinds, should consider use libm’s cosine instead.


Attached is the same code as above.


Thanks.


Zhoulai Fu

airy_inaccuratecosine.c
Description: Text Data