[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
airy_divbyzero.c
Description: airy_divbyzero.c
gsl_cos_unexpected.c
Description: gsl_cos_unexpected.c
- [Bug-gsl] Unexpected results in GSL's cosine function and airy_Ai function,
Zhoulai Fu <=