path: root/sysdeps/ieee754/ldbl-128ibm/e_lgammal_r.c

/* Natural logarithm of gamma function.  IBM Extended Precision version.
   Copyright 2001 by Stephen L. Moshier <moshier@na-net.ornl.gov>

   This library is free software; you can redistribute it and/or
   modify it under the terms of the GNU Lesser General Public
   License as published by the Free Software Foundation; either
   version 2.1 of the License, or (at your option) any later version.

   This library is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   Lesser General Public License for more details.

   You should have received a copy of the GNU Lesser General Public
   License along with this library; if not, see
   <https://www.gnu.org/licenses/>.  */

/* This file was copied from sysdeps/ieee754/ldbl-128/e_lgammal_r.c.  */


#include <math.h>
#include <math_private.h>
#include <float.h>
#include <libm-alias-finite.h>

static const long double PIL = 3.1415926535897932384626433832795028841972E0L;
static const long double MAXLGM = 0x5.d53649e2d469dbc1f01e99fd66p+1012L;
static const long double one = 1;
static const long double huge = LDBL_MAX;

/* log gamma(x) = ( x - 0.5 ) * log(x) - x + LS2PI + 1/x P(1/x^2)
   1/x <= 0.0741 (x >= 13.495...)
   Peak relative error 1.5e-36  */
static const long double ls2pi = 9.1893853320467274178032973640561763986140E-1L;
#define NRAS