[bug #58064] Inconsistent evaluation of limits in special functions

[Jackson Vanover]

URL:
  <https://savannah.gnu.org/bugs/?58064>

                 Summary: Inconsistent evaluation of limits in special
functions
                 Project: GNU Scientific Library
            Submitted by: jacksonvanover
            Submitted on: Sat 28 Mar 2020 05:20:50 PM UTC
                Category: Accuracy problem
                Severity: 3 - Normal
        Operating System: 18.04
                  Status: None
             Assigned to: None
             Open/Closed: Open
                 Release: 2.6
         Discussion Lock: Any

    _______________________________________________________

Details:

Using version 2.6 of GSL from ftp://ftp.gnu.org/gnu/gsl/gsl-2.6.tar.gz

OS is Ubuntu 18.04.2 LTS

Hardware is a Dell workstation with a 3.60 GHz Intel i7-4790 and 32 GB
of RAM

Compiler is gcc 7.4.0 with the -w option

All of the below pertain to evaluation of limits; I bring these
up as issues because in other cases, the special functions in the library
handle such limits approaching infinity properly, e.g. gsl_sf_eta(inf)
returns 1. The documentation is not clear about which functions support this
and which do not.


- When the first argument is infinity, gsl_sf_gamma_inc_Q consistently
returns 1. Meanwhile, the implementations in SciPy, mpmath, and the JMAT
library all recognize this as an erroneous computation and throw exceptions
or return NaNs; e.g. gsl_sf_gamma_inc_Q(GSL_POSINF, 22)

- When evaluating gsl_sf_gamma_inc_Q(x, inf) where x is some floating point
number, a NaN is returned even though the limit approaches 0;
e.g. gsl_sf_gamma_inc_Q(2.533265554575144, GSL_POSINF)

- When evaluating gsl_sf_erf(inf), a NaN is returned even though the limit
approaches 1, e.g. gsl_sf_erf(GSL_POSINF)

- When evaluating gsl_sf_gamma_inc(x, inf) where x is some floating point
number, a NaN is returned even though the limit approaches 0;
e.g. gsl_sf_gamma_inc(2.533265554575144, GSL_POSINF)

- When the first argument is infinity, gsl_sf_gamma_inc_P consistently
returns 0. Meanwhile, the implementations in SciPy, mpmath, and the JMAT
library all recognize this as an erroneous computation and throw exceptions
or return NaNs; e.g. gsl_sf_gamma_inc_P(GSL_POSINF, 22)

- When evaluating gsl_sf_gamma_inc_P(x, inf) where x is some floating point
number, a NaN is returned even though the limit approaches 1;
e.g. gsl_sf_gamma_inc_P(2.533265554575144, GSL_POSINF)

- When the first argument to gsl_sf_poch is a negative
infinity, that negative infinity is dropped and the function returns a
zero. The docs state that if the first argument and the sum of the two
arguments are both negative, the limiting case is returned; that
doesn't happen here; e.g. gsl_sf_poch(GSL_NEGINF, 2.2738632088209076)

- When evaluating the Si(inf), we get a return of NaN though the limit
approaches 1.570796326794896619231322 (
http://mpmath.org/doc/1.1.0/functions/expintegrals.html?highlight=si#mpmath.si
); e.g. gsl_sf_Si(GSL_POSINF)







    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?58064>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/