Hello,
I found a problem in gsl (assuming the compiler/linker is not broken),
which I would consider a bug.
A minimal working example is attached.
I compiled using g++ 6.3.0 and gsl 2.3.0. The compile command is given
in the file.
The code uses an extension to cause the code to stop on floating point
exceptions (found this via stack exchange
https://linux.die.net/man/3/feenableexcept).
I would expect the code to just run through with output something like
result: 4.39016e-164
abserr: 0
nevals: 33
Instead i get the message
Gleitkomma-Ausnahme (german for floating point exception)
Running the program with gdb I get
Program received signal SIGFPE, Arithmetic exception.
0x00007ffff7a1b42f in gsl_integration_cquad () from
/usr/lib/x86_64-linux-gnu/libgsl.so.19
(gdb) backtrace
#0 0x00007ffff7a1b42f in gsl_integration_cquad () from
/usr/lib/x86_64-linux-gnu/libgsl.so.19
#1 0x0000555555554cee in main () at main.cpp:42
Changing the lower bound to 19.2 will cause the code to work, i.e. no
floating point exception.
If I do not trap the exception the code works as expected, but if there
is no real problem, then no exception should be raised.
With kind regards
Rico Buchholz
/**
*
* g++ main.cpp -g -Wall -o main.x -L`gsl-config --libs` -lm
*
*/
#include <iostream>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <fenv.h>
#include <gsl/gsl_integration.h>
double f(double x, void * params) {
return std::exp(-x*x);
}
int main (void)
{
size_t const n(100);
double const a(19.3);
double const b(20.0);
double const epsabs(1.0e-13);
double const epsrel(1.0e-13);
double result(0.0);
double abserr(0.0);
size_t nevals(0);
/* FE_DIVBYZERO, FE_INEXACT, FE_INVALID, FE_OVERFLOW, FE_UNDERFLOW */
feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
gsl_integration_cquad_workspace * workspace(NULL);
gsl_function F;
F.function = &f;
F.params = NULL;
workspace = gsl_integration_cquad_workspace_alloc(n);
gsl_integration_cquad(&F, a, b, epsabs, epsrel, workspace, &result, &abserr, &nevals);
std::cout<<"result: " << result << "\n";
std::cout<<"abserr: " << abserr << "\n";
std::cout<<"nevals: " << nevals << "\n";
gsl_integration_cquad_workspace_free(workspace);
return 0;
}
