bug-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

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


From: Zhoulai Fu
Subject: [Bug-gsl] Unexpected results in GSL's cosine function and airy_Ai function
Date: Tue, 7 Nov 2017 10:46:11 +0000

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

Attachment: airy_divbyzero.c
Description: airy_divbyzero.c

Attachment: gsl_cos_unexpected.c
Description: gsl_cos_unexpected.c


reply via email to

[Prev in Thread] Current Thread [Next in Thread]