Thu 15 Apr 2010 02:26:24 PM UTC, original submission:
From: Wolfgang Ehrhardt <Wolfgang.Ehrhardt@munich.netsurf.de>
To: Brian Gough <bjg@network-theory.co.uk>
Subject: Bug in GSL 1.14 gsl_sf_atanint_e
Date: Wed, 14 Apr 2010 18:2:39 +0100
Dear Brian Gough,
During the development of my Pascal special function library I had a look at the
GSL/specfunc/atanint.c code. Since I am no GSL user and you seem (according to the
changelog) to be the main developer for GSL/SpecFunc I sent the following
observation
to you via email:
I think gsl_sf_atanint_e in atanint.c is slightly buggy for arguments with absolute
values >= 1.0/GSL_SQRT_DBL_EPSILON because an 1/x term is missing (in line 102).
The corresponding source code branch should be something like this (not sure about
result->err)
...
else if(ax < 1.0/GSL_SQRT_DBL_EPSILON) {
...
}
else {
result->val = sgn * (0.5M_PIlog(ax) + 1.0/ax);
result->err = 2.0 * fabs(result->val * GSL_DBL_EPSILON);
return GSL_SUCCESS;
}
And may I suggest that you add a line to specfunc/test_sf.c
TEST_SF(s, gsl_sf_atanint_e, (1.0e+9, &r), 32.552029856869591656, TEST_TOL0,
GSL_SUCCESS);
The check value is calculated with Pari/GP 2.3.4 using the definition
atanint(x)=1/2I(dilog(-Ix)-dilog(Ix))
Although I am quite sure about the issue, please note, that this is a purely
theoretical observation.
Best regards
Wolfgang Ehrhardt
(Munich, Germany)
|