[thirdparty/glibc.git] / sysdeps / ieee754 / ldbl-128 / s_log1pl.c

/*							log1pl.c
 *
 *      Relative error logarithm
 *	Natural logarithm of 1+x, 128-bit long double precision
 *
 *
 *
 * SYNOPSIS:
 *
 * long double x, y, log1pl();
 *
 * y = log1pl( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns the base e (2.718...) logarithm of 1+x.
 *
 * The argument 1+x is separated into its exponent and fractional
 * parts.  If the exponent is between -1 and +1, the logarithm
 * of the fraction is approximated by
 *
 *     log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
 *
 * Otherwise, setting  z = 2(w-1)/(w+1),
 *
 *     log(w) = z + z^3 P(z)/Q(z).
 *
 *
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      -1, 8       100000      1.9e-34     4.3e-35
 */

/* Copyright 2001 by Stephen L. Moshier

    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/>.  */


#include <float.h>
#include <math.h>
#include <math_private.h>
#include <math-underflow.h>

/* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x)
 * 1/sqrt(2) <= 1+x < sqrt(2)
 * Theoretical peak relative error = 5.3e-37,
 * relative peak error spread = 2.3e-14
 */
static const _Float128
  P12 = L(1.538612243596254322971797716843006400388E-6),
  P11 = L(4.998469661968096229986658302195402690910E-1),
  P10 = L(2.321125933898420063925789532045674660756E1),
  P9 = L(4.114517881637811823002128927449878962058E2),
  P8 = L(3.824952356185897735160588078446136783779E3),
  P7 = L(2.128857716871515081352991964243375186031E4),
  P6 = L(7.594356839258970405033155585486712125861E4),
  P5 = L(1.797628303815655343403735250238293741397E5),
  P4 = L(2.854829159639697837788887080758954924001E5),
  P3 = L(3.007007295140399532324943111654767187848E5),
  P2 = L(2.014652742082537582487669938141683759923E5),
  P1 = L(7.771154681358524243729929227226708890930E4),
  P0 = L(1.313572404063446165910279910527789794488E4),
  /* Q12 = 1.000000000000000000000000000000000000000E0L, */
  Q11 = L(4.839208193348159620282142911143429644326E1),
  Q10 = L(9.104928120962988414618126155557301584078E2),
  Q9 = L(9.147150349299596453976674231612674085381E3),
  Q8 = L(5.605842085972455027590989944010492125825E4),
  Q7 = L(2.248234257620569139969141618556349415120E5),
  Q6 = L(6.132189329546557743179177159925690841200E5),
  Q5 = L(1.158019977462989115839826904108208787040E6),
  Q4 = L(1.514882452993549494932585972882995548426E6),
  Q3 = L(1.347518538384329112529391120390701166528E6),
  Q2 = L(7.777690340007566932935753241556479363645E5),
  Q1 = L(2.626900195321832660448791748036714883242E5),
  Q0 = L(3.940717212190338497730839731583397586124E4);

/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
 * where z = 2(x-1)/(x+1)
 * 1/sqrt(2) <= x < sqrt(2)
 * Theoretical peak relative error = 1.1e-35,
 * relative peak error spread 1.1e-9
 */
static const _Float128
  R5 = L(-8.828896441624934385266096344596648080902E-1),
  R4 = L(8.057002716646055371965756206836056074715E1),
  R3 = L(-2.024301798136027039250415126250455056397E3),
  R2 = L(2.048819892795278657810231591630928516206E4),
  R1 = L(-8.977257995689735303686582344659576526998E4),
  R0 = L(1.418134209872192732479751274970992665513E5),
  /* S6 = 1.000000000000000000000000000000000000000E0L, */
  S5 = L(-1.186359407982897997337150403816839480438E2),
  S4 = L(3.998526750980007367835804959888064681098E3),
  S3 = L(-5.748542087379434595104154610899551484314E4),
  S2 = L(4.001557694070773974936904547424676279307E5),
  S1 = L(-1.332535117259762928288745111081235577029E6),
  S0 = L(1.701761051846631278975701529965589676574E6);

/* C1 + C2 = ln 2 */
static const _Float128 C1 = L(6.93145751953125E-1);
static const _Float128 C2 = L(1.428606820309417232121458176568075500134E-6);

static const _Float128 sqrth = L(0.7071067811865475244008443621048490392848);
/* ln (2^16384 * (1 - 2^-113)) */
static const _Float128 zero = 0;

_Float128
__log1pl (_Float128 xm1)
{
  _Float128 x, y, z, r, s;
  ieee854_long_double_shape_type u;
  int32_t hx;
  int e;

  /* Test for NaN or infinity input. */
  u.value = xm1;
  hx = u.parts32.w0;
  if ((hx & 0x7fffffff) >= 0x7fff0000)
    return xm1 + fabsl (xm1);

  /* log1p(+- 0) = +- 0.  */
  if (((hx & 0x7fffffff) == 0)
      && (u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0)
    return xm1;

  if ((hx & 0x7fffffff) < 0x3f8e0000)
    {
      math_check_force_underflow (xm1);
      if ((int) xm1 == 0)
	return xm1;
    }

  if (xm1 >= L(0x1p113))
    x = xm1;
  else
    x = xm1 + 1;

  /* log1p(-1) = -inf */
  if (x <= 0)
    {
      if (x == 0)
	return (-1 / zero);  /* log1p(-1) = -inf */
      else
	return (zero / (x - x));
    }

  /* Separate mantissa from exponent.  */

  /* Use frexp used so that denormal numbers will be handled properly.  */
  x = __frexpl (x, &e);

  /* Logarithm using log(x) = z + z^3 P(z^2)/Q(z^2),
     where z = 2(x-1)/x+1).  */
  if ((e > 2) || (e < -2))
    {
      if (x < sqrth)
	{			/* 2( 2x-1 )/( 2x+1 ) */
	  e -= 1;
	  z = x - L(0.5);
	  y = L(0.5) * z + L(0.5);
	}
      else
	{			/*  2 (x-1)/(x+1)   */
	  z = x - L(0.5);
	  z -= L(0.5);
	  y = L(0.5) * x + L(0.5);
	}
      x = z / y;
      z = x * x;
      r = ((((R5 * z
	      + R4) * z
	     + R3) * z
	    + R2) * z
	   + R1) * z
	+ R0;
      s = (((((z
	       + S5) * z
	      + S4) * z
	     + S3) * z
	    + S2) * z
	   + S1) * z
	+ S0;
      z = x * (z * r / s);
      z = z + e * C2;
      z = z + x;
      z = z + e * C1;
      return (z);
    }


  /* Logarithm using log(1+x) = x - .5x^2 + x^3 P(x)/Q(x). */

  if (x < sqrth)
    {
      e -= 1;
      if (e != 0)
	x = 2 * x - 1;	/*  2x - 1  */
      else
	x = xm1;
    }
  else
    {
      if (e != 0)
	x = x - 1;
      else
	x = xm1;
    }
  z = x * x;
  r = (((((((((((P12 * x
		 + P11) * x
		+ P10) * x
	       + P9) * x
	      + P8) * x
	     + P7) * x
	    + P6) * x
	   + P5) * x
	  + P4) * x
	 + P3) * x
	+ P2) * x
       + P1) * x
    + P0;
  s = (((((((((((x
		 + Q11) * x
		+ Q10) * x
	       + Q9) * x
	      + Q8) * x
	     + Q7) * x
	    + Q6) * x
	   + Q5) * x
	  + Q4) * x
	 + Q3) * x
	+ Q2) * x
       + Q1) * x
    + Q0;
  y = x * (z * r / s);
  y = y + e * C2;
  z = y - L(0.5) * z;
  z = z + x;
  z = z + e * C1;
  return (z);
}
Commit	Line	Data
90b828e6 AJ	1	/* log1pl.c
	2	*
	3	* Relative error logarithm
	4	* Natural logarithm of 1+x, 128-bit long double precision
	5	*
	6	*
	7	*
	8	* SYNOPSIS:
	9	*
	10	* long double x, y, log1pl();
	11	*
	12	* y = log1pl( x );
	13	*
	14	*
	15	*
	16	* DESCRIPTION:
	17	*
	18	* Returns the base e (2.718...) logarithm of 1+x.
	19	*
	20	* The argument 1+x is separated into its exponent and fractional
	21	* parts. If the exponent is between -1 and +1, the logarithm
	22	* of the fraction is approximated by
	23	*
	24	* log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
	25	*
	26	* Otherwise, setting z = 2(w-1)/(w+1),
	27	*
	28	* log(w) = z + z^3 P(z)/Q(z).
	29	*
	30	*
	31	*
	32	* ACCURACY:
	33	*
	34	* Relative error:
	35	* arithmetic domain # trials peak rms
	36	* IEEE -1, 8 100000 1.9e-34 4.3e-35
	37	*/
	38
9c84384c	39	/* Copyright 2001 by Stephen L. Moshier
cc7375ce RM	40
	41	This library is free software; you can redistribute it and/or
	42	modify it under the terms of the GNU Lesser General Public
	43	License as published by the Free Software Foundation; either
	44	version 2.1 of the License, or (at your option) any later version.
	45
	46	This library is distributed in the hope that it will be useful,
	47	but WITHOUT ANY WARRANTY; without even the implied warranty of
	48	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
	49	Lesser General Public License for more details.
	50
	51	You should have received a copy of the GNU Lesser General Public
59ba27a6	52	License along with this library; if not, see
5a82c748	53	<https://www.gnu.org/licenses/>. */
cc7375ce	54
90b828e6	55
0b7a5f92	56	#include <float.h>
1ed0291c RH	57	#include <math.h>
1ed0291c RH	58	#include <math_private.h>
8f5b00d3	59	#include <math-underflow.h>
90b828e6 AJ	60
	61	/* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x)
	62	* 1/sqrt(2) <= 1+x < sqrt(2)
	63	* Theoretical peak relative error = 5.3e-37,
	64	* relative peak error spread = 2.3e-14
	65	*/
15089e04	66	static const _Float128
02bbfb41 PM	67	P12 = L(1.538612243596254322971797716843006400388E-6),
	68	P11 = L(4.998469661968096229986658302195402690910E-1),
	69	P10 = L(2.321125933898420063925789532045674660756E1),
	70	P9 = L(4.114517881637811823002128927449878962058E2),
	71	P8 = L(3.824952356185897735160588078446136783779E3),
	72	P7 = L(2.128857716871515081352991964243375186031E4),
	73	P6 = L(7.594356839258970405033155585486712125861E4),
	74	P5 = L(1.797628303815655343403735250238293741397E5),
	75	P4 = L(2.854829159639697837788887080758954924001E5),
	76	P3 = L(3.007007295140399532324943111654767187848E5),
	77	P2 = L(2.014652742082537582487669938141683759923E5),
	78	P1 = L(7.771154681358524243729929227226708890930E4),
	79	P0 = L(1.313572404063446165910279910527789794488E4),
90b828e6	80	/* Q12 = 1.000000000000000000000000000000000000000E0L, */
02bbfb41 PM	81	Q11 = L(4.839208193348159620282142911143429644326E1),
	82	Q10 = L(9.104928120962988414618126155557301584078E2),
	83	Q9 = L(9.147150349299596453976674231612674085381E3),
	84	Q8 = L(5.605842085972455027590989944010492125825E4),
	85	Q7 = L(2.248234257620569139969141618556349415120E5),
	86	Q6 = L(6.132189329546557743179177159925690841200E5),
	87	Q5 = L(1.158019977462989115839826904108208787040E6),
	88	Q4 = L(1.514882452993549494932585972882995548426E6),
	89	Q3 = L(1.347518538384329112529391120390701166528E6),
	90	Q2 = L(7.777690340007566932935753241556479363645E5),
	91	Q1 = L(2.626900195321832660448791748036714883242E5),
	92	Q0 = L(3.940717212190338497730839731583397586124E4);
90b828e6 AJ	93
	94	/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
	95	* where z = 2(x-1)/(x+1)
	96	* 1/sqrt(2) <= x < sqrt(2)
	97	* Theoretical peak relative error = 1.1e-35,
	98	* relative peak error spread 1.1e-9
	99	*/
15089e04	100	static const _Float128
02bbfb41 PM	101	R5 = L(-8.828896441624934385266096344596648080902E-1),
	102	R4 = L(8.057002716646055371965756206836056074715E1),
	103	R3 = L(-2.024301798136027039250415126250455056397E3),
	104	R2 = L(2.048819892795278657810231591630928516206E4),
	105	R1 = L(-8.977257995689735303686582344659576526998E4),
	106	R0 = L(1.418134209872192732479751274970992665513E5),
90b828e6	107	/* S6 = 1.000000000000000000000000000000000000000E0L, */
02bbfb41 PM	108	S5 = L(-1.186359407982897997337150403816839480438E2),
	109	S4 = L(3.998526750980007367835804959888064681098E3),
	110	S3 = L(-5.748542087379434595104154610899551484314E4),
	111	S2 = L(4.001557694070773974936904547424676279307E5),
	112	S1 = L(-1.332535117259762928288745111081235577029E6),
	113	S0 = L(1.701761051846631278975701529965589676574E6);
90b828e6 AJ	114
90b828e6 AJ	115	/* C1 + C2 = ln 2 */
02bbfb41 PM	116	static const _Float128 C1 = L(6.93145751953125E-1);
02bbfb41 PM	117	static const _Float128 C2 = L(1.428606820309417232121458176568075500134E-6);
90b828e6	118
02bbfb41	119	static const _Float128 sqrth = L(0.7071067811865475244008443621048490392848);
90b828e6	120	/* ln (2^16384 * (1 - 2^-113)) */
02bbfb41	121	static const _Float128 zero = 0;
90b828e6	122
15089e04 PM	123	_Float128
15089e04 PM	124	__log1pl (_Float128 xm1)
90b828e6	125	{
15089e04	126	_Float128 x, y, z, r, s;
90b828e6	127	ieee854_long_double_shape_type u;
52e1b618	128	int32_t hx;
90b828e6 AJ	129	int e;
90b828e6 AJ	130
bdce812b	131	/* Test for NaN or infinity input. */
90b828e6	132	u.value = xm1;
52e1b618	133	hx = u.parts32.w0;
af1b2fd0 JM	134	if ((hx & 0x7fffffff) >= 0x7fff0000)
af1b2fd0 JM	135	return xm1 + fabsl (xm1);
bdce812b AJ	136
bdce812b AJ	137	/* log1p(+- 0) = +- 0. */
52e1b618 UD	138	if (((hx & 0x7fffffff) == 0)
52e1b618 UD	139	&& (u.parts32.w1 \| u.parts32.w2 \| u.parts32.w3) == 0)
bdce812b AJ	140	return xm1;
bdce812b AJ	141
447885eb DM	142	if ((hx & 0x7fffffff) < 0x3f8e0000)
447885eb DM	143	{
d96164c3	144	math_check_force_underflow (xm1);
447885eb DM	145	if ((int) xm1 == 0)
	146	return xm1;
	147	}
	148
02bbfb41	149	if (xm1 >= L(0x1p113))
1a84c3d6 JM	150	x = xm1;
1a84c3d6 JM	151	else
02bbfb41	152	x = xm1 + 1;
90b828e6	153
bdce812b	154	/* log1p(-1) = -inf */
02bbfb41	155	if (x <= 0)
90b828e6	156	{
02bbfb41 PM	157	if (x == 0)
02bbfb41 PM	158	return (-1 / zero); /* log1p(-1) = -inf */
90b828e6	159	else
52e1b618	160	return (zero / (x - x));
90b828e6 AJ	161	}
	162
	163	/* Separate mantissa from exponent. */
	164
	165	/* Use frexp used so that denormal numbers will be handled properly. */
c5ee217f	166	x = __frexpl (x, &e);
90b828e6 AJ	167
	168	/* Logarithm using log(x) = z + z^3 P(z^2)/Q(z^2),
	169	where z = 2(x-1)/x+1). */
	170	if ((e > 2) \|\| (e < -2))
	171	{
	172	if (x < sqrth)
	173	{ /* 2( 2x-1 )/( 2x+1 ) */
	174	e -= 1;
02bbfb41 PM	175	z = x - L(0.5);
02bbfb41 PM	176	y = L(0.5) * z + L(0.5);
90b828e6 AJ	177	}
	178	else
	179	{ /* 2 (x-1)/(x+1) */
02bbfb41 PM	180	z = x - L(0.5);
	181	z -= L(0.5);
	182	y = L(0.5) * x + L(0.5);
90b828e6 AJ	183	}
	184	x = z / y;
	185	z = x * x;
	186	r = ((((R5 * z
	187	+ R4) * z
	188	+ R3) * z
	189	+ R2) * z
	190	+ R1) * z
	191	+ R0;
	192	s = (((((z
	193	+ S5) * z
	194	+ S4) * z
	195	+ S3) * z
	196	+ S2) * z
	197	+ S1) * z
	198	+ S0;
	199	z = x * (z * r / s);
	200	z = z + e * C2;
	201	z = z + x;
	202	z = z + e * C1;
	203	return (z);
	204	}
	205
	206
	207	/* Logarithm using log(1+x) = x - .5x^2 + x^3 P(x)/Q(x). */
	208
	209	if (x < sqrth)
	210	{
	211	e -= 1;
	212	if (e != 0)
02bbfb41	213	x = 2 * x - 1; /* 2x - 1 */
90b828e6 AJ	214	else
	215	x = xm1;
	216	}
	217	else
	218	{
	219	if (e != 0)
02bbfb41	220	x = x - 1;
90b828e6 AJ	221	else
	222	x = xm1;
	223	}
	224	z = x * x;
	225	r = (((((((((((P12 * x
	226	+ P11) * x
	227	+ P10) * x
	228	+ P9) * x
	229	+ P8) * x
	230	+ P7) * x
	231	+ P6) * x
	232	+ P5) * x
	233	+ P4) * x
	234	+ P3) * x
	235	+ P2) * x
	236	+ P1) * x
	237	+ P0;
	238	s = (((((((((((x
	239	+ Q11) * x
	240	+ Q10) * x
	241	+ Q9) * x
	242	+ Q8) * x
	243	+ Q7) * x
	244	+ Q6) * x
	245	+ Q5) * x
	246	+ Q4) * x
	247	+ Q3) * x
	248	+ Q2) * x
	249	+ Q1) * x
	250	+ Q0;
	251	y = x * (z * r / s);
	252	y = y + e * C2;
02bbfb41	253	z = y - L(0.5) * z;
90b828e6 AJ	254	z = z + x;
	255	z = z + e * C1;
	256	return (z);
	257	}