[thirdparty/glibc.git] / sysdeps / ieee754 / ldbl-96 / e_gammal_r.c

/* Implementation of gamma function according to ISO C.
   Copyright (C) 1997-2021 Free Software Foundation, Inc.
   This file is part of the GNU C Library.
   Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.

   The GNU C 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.

   The GNU C 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 the GNU C Library; if not, see
   <https://www.gnu.org/licenses/>.  */

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

/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's
   approximation to gamma function.  */

static const long double gamma_coeff[] =
  {
    0x1.5555555555555556p-4L,
    -0xb.60b60b60b60b60bp-12L,
    0x3.4034034034034034p-12L,
    -0x2.7027027027027028p-12L,
    0x3.72a3c5631fe46aep-12L,
    -0x7.daac36664f1f208p-12L,
    0x1.a41a41a41a41a41ap-8L,
    -0x7.90a1b2c3d4e5f708p-8L,
  };

#define NCOEFF (sizeof (gamma_coeff) / sizeof (gamma_coeff[0]))

/* Return gamma (X), for positive X less than 1766, in the form R *
   2^(*EXP2_ADJ), where R is the return value and *EXP2_ADJ is set to
   avoid overflow or underflow in intermediate calculations.  */

static long double
gammal_positive (long double x, int *exp2_adj)
{
  int local_signgam;
  if (x < 0.5L)
    {
      *exp2_adj = 0;
      return __ieee754_expl (__ieee754_lgammal_r (x + 1, &local_signgam)) / x;
    }
  else if (x <= 1.5L)
    {
      *exp2_adj = 0;
      return __ieee754_expl (__ieee754_lgammal_r (x, &local_signgam));
    }
  else if (x < 7.5L)
    {
      /* Adjust into the range for using exp (lgamma).  */
      *exp2_adj = 0;
      long double n = ceill (x - 1.5L);
      long double x_adj = x - n;
      long double eps;
      long double prod = __gamma_productl (x_adj, 0, n, &eps);
      return (__ieee754_expl (__ieee754_lgammal_r (x_adj, &local_signgam))
	      * prod * (1.0L + eps));
    }
  else
    {
      long double eps = 0;
      long double x_eps = 0;
      long double x_adj = x;
      long double prod = 1;
      if (x < 13.0L)
	{
	  /* Adjust into the range for applying Stirling's
	     approximation.  */
	  long double n = ceill (13.0L - x);
	  x_adj = x + n;
	  x_eps = (x - (x_adj - n));
	  prod = __gamma_productl (x_adj - n, x_eps, n, &eps);
	}
      /* The result is now gamma (X_ADJ + X_EPS) / (PROD * (1 + EPS)).
	 Compute gamma (X_ADJ + X_EPS) using Stirling's approximation,
	 starting by computing pow (X_ADJ, X_ADJ) with a power of 2
	 factored out.  */
      long double exp_adj = -eps;
      long double x_adj_int = roundl (x_adj);
      long double x_adj_frac = x_adj - x_adj_int;
      int x_adj_log2;
      long double x_adj_mant = __frexpl (x_adj, &x_adj_log2);
      if (x_adj_mant < M_SQRT1_2l)
	{
	  x_adj_log2--;
	  x_adj_mant *= 2.0L;
	}
      *exp2_adj = x_adj_log2 * (int) x_adj_int;
      long double ret = (__ieee754_powl (x_adj_mant, x_adj)
			 * __ieee754_exp2l (x_adj_log2 * x_adj_frac)
			 * __ieee754_expl (-x_adj)
			 * sqrtl (2 * M_PIl / x_adj)
			 / prod);
      exp_adj += x_eps * __ieee754_logl (x_adj);
      long double bsum = gamma_coeff[NCOEFF - 1];
      long double x_adj2 = x_adj * x_adj;
      for (size_t i = 1; i <= NCOEFF - 1; i++)
	bsum = bsum / x_adj2 + gamma_coeff[NCOEFF - 1 - i];
      exp_adj += bsum / x_adj;
      return ret + ret * __expm1l (exp_adj);
    }
}

long double
__ieee754_gammal_r (long double x, int *signgamp)
{
  uint32_t es, hx, lx;
  long double ret;

  GET_LDOUBLE_WORDS (es, hx, lx, x);

  if (__glibc_unlikely (((es & 0x7fff) | hx | lx) == 0))
    {
      /* Return value for x == 0 is Inf with divide by zero exception.  */
      *signgamp = 0;
      return 1.0 / x;
    }
  if (__glibc_unlikely (es == 0xffffffff && ((hx & 0x7fffffff) | lx) == 0))
    {
      /* x == -Inf.  According to ISO this is NaN.  */
      *signgamp = 0;
      return x - x;
    }
  if (__glibc_unlikely ((es & 0x7fff) == 0x7fff))
    {
      /* Positive infinity (return positive infinity) or NaN (return
	 NaN).  */
      *signgamp = 0;
      return x + x;
    }
  if (__builtin_expect ((es & 0x8000) != 0, 0) && rintl (x) == x)
    {
      /* Return value for integer x < 0 is NaN with invalid exception.  */
      *signgamp = 0;
      return (x - x) / (x - x);
    }

  if (x >= 1756.0L)
    {
      /* Overflow.  */
      *signgamp = 0;
      return LDBL_MAX * LDBL_MAX;
    }
  else
    {
      SET_RESTORE_ROUNDL (FE_TONEAREST);
      if (x > 0.0L)
	{
	  *signgamp = 0;
	  int exp2_adj;
	  ret = gammal_positive (x, &exp2_adj);
	  ret = __scalbnl (ret, exp2_adj);
	}
      else if (x >= -LDBL_EPSILON / 4.0L)
	{
	  *signgamp = 0;
	  ret = 1.0L / x;
	}
      else
	{
	  long double tx = truncl (x);
	  *signgamp = (tx == 2.0L * truncl (tx / 2.0L)) ? -1 : 1;
	  if (x <= -1766.0L)
	    /* Underflow.  */
	    ret = LDBL_MIN * LDBL_MIN;
	  else
	    {
	      long double frac = tx - x;
	      if (frac > 0.5L)
		frac = 1.0L - frac;
	      long double sinpix = (frac <= 0.25L
				    ? __sinl (M_PIl * frac)
				    : __cosl (M_PIl * (0.5L - frac)));
	      int exp2_adj;
	      ret = M_PIl / (-x * sinpix
			     * gammal_positive (-x, &exp2_adj));
	      ret = __scalbnl (ret, -exp2_adj);
	      math_check_force_underflow_nonneg (ret);
	    }
	}
    }
  if (isinf (ret) && x != 0)
    {
      if (*signgamp < 0)
	return -(-copysignl (LDBL_MAX, ret) * LDBL_MAX);
      else
	return copysignl (LDBL_MAX, ret) * LDBL_MAX;
    }
  else if (ret == 0)
    {
      if (*signgamp < 0)
	return -(-copysignl (LDBL_MIN, ret) * LDBL_MIN);
      else
	return copysignl (LDBL_MIN, ret) * LDBL_MIN;
    }
  else
    return ret;
}
libm_alias_finite (__ieee754_gammal_r, __gammal_r)
Commit	Line	Data
d705269e	1	/* Implementation of gamma function according to ISO C.
2b778ceb	2	Copyright (C) 1997-2021 Free Software Foundation, Inc.
c131718c UD	3	This file is part of the GNU C Library.
	4	Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
	5
	6	The GNU C Library is free software; you can redistribute it and/or
41bdb6e2 AJ	7	modify it under the terms of the GNU Lesser General Public
	8	License as published by the Free Software Foundation; either
	9	version 2.1 of the License, or (at your option) any later version.
c131718c UD	10
	11	The GNU C Library is distributed in the hope that it will be useful,
	12	but WITHOUT ANY WARRANTY; without even the implied warranty of
	13	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
41bdb6e2	14	Lesser General Public License for more details.
c131718c	15
41bdb6e2	16	You should have received a copy of the GNU Lesser General Public
59ba27a6	17	License along with the GNU C Library; if not, see
5a82c748	18	<https://www.gnu.org/licenses/>. */
c131718c	19
d705269e UD	20	#include <math.h>
d705269e UD	21	#include <math_private.h>
70e2ba33	22	#include <fenv_private.h>
8f5b00d3	23	#include <math-underflow.h>
d8cd06db	24	#include <float.h>
220622dd	25	#include <libm-alias-finite.h>
d705269e	26
d8cd06db JM	27	/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's
	28	approximation to gamma function. */
	29
	30	static const long double gamma_coeff[] =
	31	{
	32	0x1.5555555555555556p-4L,
	33	-0xb.60b60b60b60b60bp-12L,
	34	0x3.4034034034034034p-12L,
	35	-0x2.7027027027027028p-12L,
	36	0x3.72a3c5631fe46aep-12L,
	37	-0x7.daac36664f1f208p-12L,
	38	0x1.a41a41a41a41a41ap-8L,
	39	-0x7.90a1b2c3d4e5f708p-8L,
	40	};
	41
	42	#define NCOEFF (sizeof (gamma_coeff) / sizeof (gamma_coeff[0]))
	43
	44	/* Return gamma (X), for positive X less than 1766, in the form R *
	45	2^(EXP2_ADJ), where R is the return value and EXP2_ADJ is set to
	46	avoid overflow or underflow in intermediate calculations. */
	47
	48	static long double
	49	gammal_positive (long double x, int *exp2_adj)
	50	{
	51	int local_signgam;
	52	if (x < 0.5L)
	53	{
	54	*exp2_adj = 0;
	55	return __ieee754_expl (__ieee754_lgammal_r (x + 1, &local_signgam)) / x;
	56	}
	57	else if (x <= 1.5L)
	58	{
	59	*exp2_adj = 0;
	60	return __ieee754_expl (__ieee754_lgammal_r (x, &local_signgam));
	61	}
	62	else if (x < 7.5L)
	63	{
	64	/* Adjust into the range for using exp (lgamma). */
	65	*exp2_adj = 0;
71223ef9	66	long double n = ceill (x - 1.5L);
d8cd06db JM	67	long double x_adj = x - n;
	68	long double eps;
	69	long double prod = __gamma_productl (x_adj, 0, n, &eps);
	70	return (__ieee754_expl (__ieee754_lgammal_r (x_adj, &local_signgam))
	71	* prod * (1.0L + eps));
	72	}
	73	else
	74	{
	75	long double eps = 0;
	76	long double x_eps = 0;
	77	long double x_adj = x;
	78	long double prod = 1;
	79	if (x < 13.0L)
	80	{
	81	/* Adjust into the range for applying Stirling's
	82	approximation. */
71223ef9	83	long double n = ceill (13.0L - x);
d8cd06db JM	84	x_adj = x + n;
	85	x_eps = (x - (x_adj - n));
	86	prod = __gamma_productl (x_adj - n, x_eps, n, &eps);
	87	}
	88	/* The result is now gamma (X_ADJ + X_EPS) / (PROD * (1 + EPS)).
	89	Compute gamma (X_ADJ + X_EPS) using Stirling's approximation,
	90	starting by computing pow (X_ADJ, X_ADJ) with a power of 2
	91	factored out. */
	92	long double exp_adj = -eps;
9755bc46	93	long double x_adj_int = roundl (x_adj);
d8cd06db JM	94	long double x_adj_frac = x_adj - x_adj_int;
	95	int x_adj_log2;
	96	long double x_adj_mant = __frexpl (x_adj, &x_adj_log2);
	97	if (x_adj_mant < M_SQRT1_2l)
	98	{
	99	x_adj_log2--;
	100	x_adj_mant *= 2.0L;
	101	}
	102	exp2_adj = x_adj_log2 (int) x_adj_int;
	103	long double ret = (__ieee754_powl (x_adj_mant, x_adj)
	104	* __ieee754_exp2l (x_adj_log2 * x_adj_frac)
	105	* __ieee754_expl (-x_adj)
f67a8147	106	* sqrtl (2 * M_PIl / x_adj)
d8cd06db	107	/ prod);
e02920bc	108	exp_adj += x_eps * __ieee754_logl (x_adj);
d8cd06db JM	109	long double bsum = gamma_coeff[NCOEFF - 1];
	110	long double x_adj2 = x_adj * x_adj;
	111	for (size_t i = 1; i <= NCOEFF - 1; i++)
	112	bsum = bsum / x_adj2 + gamma_coeff[NCOEFF - 1 - i];
	113	exp_adj += bsum / x_adj;
	114	return ret + ret * __expm1l (exp_adj);
	115	}
	116	}
d705269e UD	117
	118	long double
	119	__ieee754_gammal_r (long double x, int *signgamp)
	120	{
24ab7723	121	uint32_t es, hx, lx;
e02920bc	122	long double ret;
d705269e UD	123
	124	GET_LDOUBLE_WORDS (es, hx, lx, x);
	125
a1ffb40e	126	if (__glibc_unlikely (((es & 0x7fff) \| hx \| lx) == 0))
b3fc5f84	127	{
52495f29	128	/* Return value for x == 0 is Inf with divide by zero exception. */
b3fc5f84	129	*signgamp = 0;
52495f29	130	return 1.0 / x;
b3fc5f84	131	}
a1ffb40e	132	if (__glibc_unlikely (es == 0xffffffff && ((hx & 0x7fffffff) \| lx) == 0))
3bde1a69 UD	133	{
	134	/* x == -Inf. According to ISO this is NaN. */
	135	*signgamp = 0;
	136	return x - x;
	137	}
a1ffb40e	138	if (__glibc_unlikely ((es & 0x7fff) == 0x7fff))
abefbc51	139	{
d8cd06db JM	140	/* Positive infinity (return positive infinity) or NaN (return
d8cd06db JM	141	NaN). */
abefbc51	142	*signgamp = 0;
d8cd06db	143	return x + x;
abefbc51	144	}
f29b6f17	145	if (__builtin_expect ((es & 0x8000) != 0, 0) && rintl (x) == x)
276ae1f2 UD	146	{
	147	/* Return value for integer x < 0 is NaN with invalid exception. */
	148	*signgamp = 0;
	149	return (x - x) / (x - x);
	150	}
d705269e	151
d8cd06db JM	152	if (x >= 1756.0L)
	153	{
	154	/* Overflow. */
	155	*signgamp = 0;
	156	return LDBL_MAX * LDBL_MAX;
	157	}
e02920bc	158	else
d8cd06db	159	{
e02920bc JM	160	SET_RESTORE_ROUNDL (FE_TONEAREST);
	161	if (x > 0.0L)
	162	{
	163	*signgamp = 0;
	164	int exp2_adj;
	165	ret = gammal_positive (x, &exp2_adj);
	166	ret = __scalbnl (ret, exp2_adj);
	167	}
	168	else if (x >= -LDBL_EPSILON / 4.0L)
	169	{
	170	*signgamp = 0;
	171	ret = 1.0L / x;
	172	}
	173	else
	174	{
7abf97be JM	175	long double tx = truncl (x);
7abf97be JM	176	signgamp = (tx == 2.0L truncl (tx / 2.0L)) ? -1 : 1;
e02920bc JM	177	if (x <= -1766.0L)
	178	/* Underflow. */
	179	ret = LDBL_MIN * LDBL_MIN;
	180	else
	181	{
	182	long double frac = tx - x;
	183	if (frac > 0.5L)
	184	frac = 1.0L - frac;
	185	long double sinpix = (frac <= 0.25L
	186	? __sinl (M_PIl * frac)
	187	: __cosl (M_PIl * (0.5L - frac)));
	188	int exp2_adj;
	189	ret = M_PIl / (-x * sinpix
	190	* gammal_positive (-x, &exp2_adj));
	191	ret = __scalbnl (ret, -exp2_adj);
d96164c3	192	math_check_force_underflow_nonneg (ret);
e02920bc JM	193	}
e02920bc JM	194	}
d8cd06db	195	}
e02920bc	196	if (isinf (ret) && x != 0)
d8cd06db	197	{
e02920bc	198	if (*signgamp < 0)
81dca813	199	return -(-copysignl (LDBL_MAX, ret) * LDBL_MAX);
e02920bc	200	else
81dca813	201	return copysignl (LDBL_MAX, ret) * LDBL_MAX;
d8cd06db	202	}
e02920bc	203	else if (ret == 0)
d8cd06db	204	{
e02920bc	205	if (*signgamp < 0)
81dca813	206	return -(-copysignl (LDBL_MIN, ret) * LDBL_MIN);
e02920bc	207	else
81dca813	208	return copysignl (LDBL_MIN, ret) * LDBL_MIN;
d8cd06db	209	}
e02920bc JM	210	else
e02920bc JM	211	return ret;
d705269e	212	}
220622dd	213	libm_alias_finite (__ieee754_gammal_r, __gammal_r)