[Bug-gsl] [bug #29562] missing 1/x term for gsl_sf_atanint_e with large argument

[Brian Gough]

URL:
  <http://savannah.gnu.org/bugs/?29562>

                 Summary: missing 1/x term for gsl_sf_atanint_e with large
argument
                 Project: GNU Scientific Library
            Submitted by: bjg
            Submitted on: Thu 15 Apr 2010 03:26:24 PM BST
                Category: Accuracy problem
                Severity: 3 - Normal
        Operating System: 
                  Status: Confirmed
             Assigned to: bjg
             Open/Closed: Open
                 Release: 1.14
         Discussion Lock: Any

    _______________________________________________________

Details:

From: Wolfgang Ehrhardt <address@hidden>
To: Brian Gough <address@hidden>
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.5*M_PI*log(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/2*I*(dilog(-I*x)-dilog(I*x))


Although I am quite sure about the issue, please note, that this is a purely

theoretical observation.

Best regards

Wolfgang Ehrhardt
(Munich, Germany)





    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?29562>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/