[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
    <http://www.gnu.org/licenses/>.  */


#include <math.h>
#include <math_private.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 long double
  P12 = 1.538612243596254322971797716843006400388E-6L,
  P11 = 4.998469661968096229986658302195402690910E-1L,
  P10 = 2.321125933898420063925789532045674660756E1L,
  P9 = 4.114517881637811823002128927449878962058E2L,
  P8 = 3.824952356185897735160588078446136783779E3L,
  P7 = 2.128857716871515081352991964243375186031E4L,
  P6 = 7.594356839258970405033155585486712125861E4L,
  P5 = 1.797628303815655343403735250238293741397E5L,
  P4 = 2.854829159639697837788887080758954924001E5L,
  P3 = 3.007007295140399532324943111654767187848E5L,
  P2 = 2.014652742082537582487669938141683759923E5L,
  P1 = 7.771154681358524243729929227226708890930E4L,
  P0 = 1.313572404063446165910279910527789794488E4L,
  /* Q12 = 1.000000000000000000000000000000000000000E0L, */
  Q11 = 4.839208193348159620282142911143429644326E1L,
  Q10 = 9.104928120962988414618126155557301584078E2L,
  Q9 = 9.147150349299596453976674231612674085381E3L,
  Q8 = 5.605842085972455027590989944010492125825E4L,
  Q7 = 2.248234257620569139969141618556349415120E5L,
  Q6 = 6.132189329546557743179177159925690841200E5L,
  Q5 = 1.158019977462989115839826904108208787040E6L,
  Q4 = 1.514882452993549494932585972882995548426E6L,
  Q3 = 1.347518538384329112529391120390701166528E6L,
  Q2 = 7.777690340007566932935753241556479363645E5L,
  Q1 = 2.626900195321832660448791748036714883242E5L,
  Q0 = 3.940717212190338497730839731583397586124E4L;

/* 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 long double
  R5 = -8.828896441624934385266096344596648080902E-1L,
  R4 = 8.057002716646055371965756206836056074715E1L,
  R3 = -2.024301798136027039250415126250455056397E3L,
  R2 = 2.048819892795278657810231591630928516206E4L,
  R1 = -8.977257995689735303686582344659576526998E4L,
  R0 = 1.418134209872192732479751274970992665513E5L,
  /* S6 = 1.000000000000000000000000000000000000000E0L, */
  S5 = -1.186359407982897997337150403816839480438E2L,
  S4 = 3.998526750980007367835804959888064681098E3L,
  S3 = -5.748542087379434595104154610899551484314E4L,
  S2 = 4.001557694070773974936904547424676279307E5L,
  S1 = -1.332535117259762928288745111081235577029E6L,
  S0 = 1.701761051846631278975701529965589676574E6L;

/* C1 + C2 = ln 2 */
static const long double C1 = 6.93145751953125E-1L;
static const long double C2 = 1.428606820309417232121458176568075500134E-6L;

static const long double sqrth = 0.7071067811865475244008443621048490392848L;
/* ln (2^16384 * (1 - 2^-113)) */
static const long double maxlog = 1.1356523406294143949491931077970764891253E4L;
static const long double zero = 0.0L;

long double
__log1pl (long double xm1)
{
  long double 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 >= 0x7fff0000)
    return xm1;

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

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

  x = xm1 + 1.0L;

  /* log1p(-1) = -inf */
  if (x <= 0.0L)
    {
      if (x == 0.0L)
	return (-1.0L / 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 - 0.5L;
	  y = 0.5L * z + 0.5L;
	}
      else
	{			/*  2 (x-1)/(x+1)   */
	  z = x - 0.5L;
	  z -= 0.5L;
	  y = 0.5L * x + 0.5L;
	}
      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.0L * x - 1.0L;	/*  2x - 1  */
      else
	x = xm1;
    }
  else
    {
      if (e != 0)
	x = x - 1.0L;
      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 - 0.5L * z;
  z = z + x;
  z = z + e * C1;
  return (z);
}

weak_alias (__log1pl, log1pl)
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 PE	52	License along with this library; if not, see
59ba27a6 PE	53	<http://www.gnu.org/licenses/>. */
cc7375ce	54
90b828e6	55
1ed0291c RH	56	#include <math.h>
1ed0291c RH	57	#include <math_private.h>
90b828e6 AJ	58
	59	/* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x)
	60	* 1/sqrt(2) <= 1+x < sqrt(2)
	61	* Theoretical peak relative error = 5.3e-37,
	62	* relative peak error spread = 2.3e-14
	63	*/
1f5649f8	64	static const long double
90b828e6 AJ	65	P12 = 1.538612243596254322971797716843006400388E-6L,
	66	P11 = 4.998469661968096229986658302195402690910E-1L,
	67	P10 = 2.321125933898420063925789532045674660756E1L,
	68	P9 = 4.114517881637811823002128927449878962058E2L,
	69	P8 = 3.824952356185897735160588078446136783779E3L,
	70	P7 = 2.128857716871515081352991964243375186031E4L,
	71	P6 = 7.594356839258970405033155585486712125861E4L,
	72	P5 = 1.797628303815655343403735250238293741397E5L,
	73	P4 = 2.854829159639697837788887080758954924001E5L,
	74	P3 = 3.007007295140399532324943111654767187848E5L,
	75	P2 = 2.014652742082537582487669938141683759923E5L,
	76	P1 = 7.771154681358524243729929227226708890930E4L,
	77	P0 = 1.313572404063446165910279910527789794488E4L,
	78	/* Q12 = 1.000000000000000000000000000000000000000E0L, */
	79	Q11 = 4.839208193348159620282142911143429644326E1L,
	80	Q10 = 9.104928120962988414618126155557301584078E2L,
	81	Q9 = 9.147150349299596453976674231612674085381E3L,
	82	Q8 = 5.605842085972455027590989944010492125825E4L,
	83	Q7 = 2.248234257620569139969141618556349415120E5L,
	84	Q6 = 6.132189329546557743179177159925690841200E5L,
	85	Q5 = 1.158019977462989115839826904108208787040E6L,
	86	Q4 = 1.514882452993549494932585972882995548426E6L,
	87	Q3 = 1.347518538384329112529391120390701166528E6L,
	88	Q2 = 7.777690340007566932935753241556479363645E5L,
	89	Q1 = 2.626900195321832660448791748036714883242E5L,
	90	Q0 = 3.940717212190338497730839731583397586124E4L;
	91
	92	/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
	93	* where z = 2(x-1)/(x+1)
	94	* 1/sqrt(2) <= x < sqrt(2)
	95	* Theoretical peak relative error = 1.1e-35,
	96	* relative peak error spread 1.1e-9
	97	*/
1f5649f8	98	static const long double
90b828e6 AJ	99	R5 = -8.828896441624934385266096344596648080902E-1L,
	100	R4 = 8.057002716646055371965756206836056074715E1L,
	101	R3 = -2.024301798136027039250415126250455056397E3L,
	102	R2 = 2.048819892795278657810231591630928516206E4L,
	103	R1 = -8.977257995689735303686582344659576526998E4L,
	104	R0 = 1.418134209872192732479751274970992665513E5L,
	105	/* S6 = 1.000000000000000000000000000000000000000E0L, */
	106	S5 = -1.186359407982897997337150403816839480438E2L,
	107	S4 = 3.998526750980007367835804959888064681098E3L,
	108	S3 = -5.748542087379434595104154610899551484314E4L,
	109	S2 = 4.001557694070773974936904547424676279307E5L,
	110	S1 = -1.332535117259762928288745111081235577029E6L,
	111	S0 = 1.701761051846631278975701529965589676574E6L;
	112
	113	/* C1 + C2 = ln 2 */
1f5649f8 UD	114	static const long double C1 = 6.93145751953125E-1L;
1f5649f8 UD	115	static const long double C2 = 1.428606820309417232121458176568075500134E-6L;
90b828e6	116
1f5649f8	117	static const long double sqrth = 0.7071067811865475244008443621048490392848L;
90b828e6	118	/* ln (2^16384 * (1 - 2^-113)) */
1f5649f8	119	static const long double maxlog = 1.1356523406294143949491931077970764891253E4L;
1f5649f8	120	static const long double zero = 0.0L;
90b828e6	121
90b828e6 AJ	122	long double
	123	__log1pl (long double xm1)
	124	{
	125	long double x, y, z, r, s;
	126	ieee854_long_double_shape_type u;
52e1b618	127	int32_t hx;
90b828e6 AJ	128	int e;
90b828e6 AJ	129
bdce812b	130	/* Test for NaN or infinity input. */
90b828e6	131	u.value = xm1;
52e1b618 UD	132	hx = u.parts32.w0;
52e1b618 UD	133	if (hx >= 0x7fff0000)
bdce812b AJ	134	return xm1;
	135
	136	/* log1p(+- 0) = +- 0. */
52e1b618 UD	137	if (((hx & 0x7fffffff) == 0)
52e1b618 UD	138	&& (u.parts32.w1 \| u.parts32.w2 \| u.parts32.w3) == 0)
bdce812b AJ	139	return xm1;
bdce812b AJ	140
447885eb DM	141	if ((hx & 0x7fffffff) < 0x3f8e0000)
	142	{
	143	if ((int) xm1 == 0)
	144	return xm1;
	145	}
	146
bdce812b	147	x = xm1 + 1.0L;
90b828e6	148
bdce812b	149	/* log1p(-1) = -inf */
90b828e6 AJ	150	if (x <= 0.0L)
	151	{
	152	if (x == 0.0L)
2ca180e9	153	return (-1.0L / zero); /* log1p(-1) = -inf */
90b828e6	154	else
52e1b618	155	return (zero / (x - x));
90b828e6 AJ	156	}
	157
	158	/* Separate mantissa from exponent. */
	159
	160	/* Use frexp used so that denormal numbers will be handled properly. */
c5ee217f	161	x = __frexpl (x, &e);
90b828e6 AJ	162
	163	/* Logarithm using log(x) = z + z^3 P(z^2)/Q(z^2),
	164	where z = 2(x-1)/x+1). */
	165	if ((e > 2) \|\| (e < -2))
	166	{
	167	if (x < sqrth)
	168	{ /* 2( 2x-1 )/( 2x+1 ) */
	169	e -= 1;
	170	z = x - 0.5L;
	171	y = 0.5L * z + 0.5L;
	172	}
	173	else
	174	{ /* 2 (x-1)/(x+1) */
	175	z = x - 0.5L;
	176	z -= 0.5L;
	177	y = 0.5L * x + 0.5L;
	178	}
	179	x = z / y;
	180	z = x * x;
	181	r = ((((R5 * z
	182	+ R4) * z
	183	+ R3) * z
	184	+ R2) * z
	185	+ R1) * z
	186	+ R0;
	187	s = (((((z
	188	+ S5) * z
	189	+ S4) * z
	190	+ S3) * z
	191	+ S2) * z
	192	+ S1) * z
	193	+ S0;
	194	z = x * (z * r / s);
	195	z = z + e * C2;
	196	z = z + x;
	197	z = z + e * C1;
	198	return (z);
	199	}
	200
	201
	202	/* Logarithm using log(1+x) = x - .5x^2 + x^3 P(x)/Q(x). */
	203
	204	if (x < sqrth)
	205	{
	206	e -= 1;
	207	if (e != 0)
	208	x = 2.0L * x - 1.0L; /* 2x - 1 */
	209	else
	210	x = xm1;
	211	}
	212	else
	213	{
	214	if (e != 0)
	215	x = x - 1.0L;
	216	else
	217	x = xm1;
	218	}
	219	z = x * x;
	220	r = (((((((((((P12 * x
	221	+ P11) * x
	222	+ P10) * x
	223	+ P9) * x
	224	+ P8) * x
	225	+ P7) * x
226	+ P6) * x
227	+ P5) * x
228	+ P4) * x
229	+ P3) * x
230	+ P2) * x
231	+ P1) * x
232	+ P0;
233	s = (((((((((((x
234	+ Q11) * x
235	+ Q10) * x
236	+ Q9) * x
237	+ Q8) * x
238	+ Q7) * x
239	+ Q6) * x
240	+ Q5) * x
241	+ Q4) * x
242	+ Q3) * x
243	+ Q2) * x
244	+ Q1) * x
245	+ Q0;
246	y = x * (z * r / s);
247	y = y + e * C2;
248	z = y - 0.5L * z;
249	z = z + x;
250	z = z + e * C1;
251	return (z);
252	}
253
254	weak_alias (__log1pl, log1pl)