Re: [Bug-gsl] Unexpected results in GSL's cosine function and airy_Ai function

[Patrick Alken]

Hello,

  It is known that the GSL gsl_sf_cos and gsl_sf_sin functions fail for
large inputs. It is recommended to use the libc/libm cos/sin functions.

I will log your airy bug report on the bug tracker, but it may take a
while for me or someone else to look into it. My time is quite limited
these days and there aren't so many people working on GSL anymore

Thanks,
Patrick

On 11/07/2017 11:46 AM, Zhoulai Fu wrote:
> Dear GSL developers,
>
> I have recently got some unexpected results when using GSL. Below are the 
> symptoms:
>
> 1. The GSL’s cosine function gsl_sf_cos_e  gives unexpected results on large 
> inputs. For example, it returns -inf for the input 8.1e50, or -nan for the 
> input 4.2e70. The error estimation and returned status of GSL’s cosine 
> function are also strange. Please check gsl_cos_unexpected.c (pasted below 
> and attached) to reproduce the symptom.
>
> 2., The GSL’s airy_Ai function gsl_sf_airy_Ai_e  triggers a floating-point 
> division-by-zero exception with the input -1.842761151977744. The exception 
> can be trapped by using “feenableexcept”, but it will disappear if you 
> slightly disturb the input, say, to -1.84276115198. By digging into the 
> issue, I have observed that in the airy_mod_phase function that airy_Ai 
> invokes, the variable result_m is divided whereas it has vanished following a 
> non-trivial computation (with loops, etc., in function cheb_eval_mode_e). 
> Please check the code airy_divbyzero.c below.
>
> Attached are code that reproduces the symptoms described above. They are also 
> pasted below for your convenience.  Could you please let me know whether I 
> misuse or misunderstand something? For your info, my GSL version is 2.4. It 
> is compiled with gcc 4.8.2. All experiments are conducted on Ubuntu 14.04.
>
> Thanks.
> Zhoulai Fu
>
> Code gsl_cos_unexpected.c:
>
> #include <gsl/gsl_math.h>
> #include <gsl/gsl_errno.h>
> #include <gsl/gsl_sf_trig.h>
> #include <math.h>
>
> int main(){
>   static double inputs [] = {8.1e50, 4.2e70};
>
>   for (int i = 0; i < 2; i++){
>     gsl_sf_result res;
>     double x = inputs[i];
>     int status = gsl_sf_cos_e(x, &res);
>     printf("With input %g:\n", x);
>     printf("  gsl_cos status = %d\n",    status);
>     printf("  gsl_cos val    = %.16g\n", res.val);
>     printf("  gsl_cos err    = %g\n",    res.err);
>     printf("  But, glibc_cos = %g\n",    cos(x));
>     printf("\n");
>   }
> }
> //To reproduce the symptom, use:
> //    g++ -g -std=c++11 gsl_cos_unexpected.c -lm -lgsl -lgslcblas -o a.out -I 
> /home/parallels/Work/gsl-2.4/ -I /home/parallels/Work/gsl-2.4/specfunc/; 
> ./a.out
>
>
> Code airy_divbyzero.c:
>
> #ifndef _GNU_SOURCE
>   #define _GNU_SOURCE
> #endif
> #include <fenv.h>
> #include <stdlib.h>
> #include <gsl/gsl_math.h>
> #include <gsl/gsl_errno.h>
> #include <gsl/gsl_sf_trig.h>
> #include <gsl/gsl_sf_airy.h>
>
> int main(){
>   if (getenv("DIV"))
>     feenableexcept(FE_DIVBYZERO);
>   gsl_sf_result res;
>   double x = -1.842761151977744;
>   gsl_mode_t mode = GSL_PREC_DOUBLE;
>   int status = gsl_sf_airy_Ai_e(x, mode, &res);
>
>   printf(" status = %d\n", status);
>   printf(" result = %.16g\n", res.val);
>   printf(" error = %g\n", res.err);
> }
> //To reproduce the symptom, use:
> // g++ -g -std=c++11 airy_divbyzero.c -lm -lgsl -lgslcblas -o a.out -I 
> /home/parallels/Work/gsl-2.4/ -I /home/parallels/Work/gsl-2.4/specfunc/; 
> DIV=1 ./a.out