[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;

  x = xm1 + 1.0L;

  /* log1p(-1) = -inf */
  if (x <= 0.0L)
    {
      if (x == 0.0L)
	return (-1.0L / (x - x));
      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
cc7375ce RM	39	/* Copyright 2001 by Stephen L. Moshier
	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 AJ	55
	56	#include "math.h"
	57	#include "math_private.h"
	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;
	140
	141	x = xm1 + 1.0L;
90b828e6	142
bdce812b	143	/* log1p(-1) = -inf */
90b828e6 AJ	144	if (x <= 0.0L)
	145	{
	146	if (x == 0.0L)
52e1b618	147	return (-1.0L / (x - x));
90b828e6	148	else
52e1b618	149	return (zero / (x - x));
90b828e6 AJ	150	}
	151
	152	/* Separate mantissa from exponent. */
	153
	154	/* Use frexp used so that denormal numbers will be handled properly. */
c5ee217f	155	x = __frexpl (x, &e);
90b828e6 AJ	156
	157	/* Logarithm using log(x) = z + z^3 P(z^2)/Q(z^2),
	158	where z = 2(x-1)/x+1). */
	159	if ((e > 2) \|\| (e < -2))
	160	{
	161	if (x < sqrth)
	162	{ /* 2( 2x-1 )/( 2x+1 ) */
	163	e -= 1;
	164	z = x - 0.5L;
	165	y = 0.5L * z + 0.5L;
	166	}
	167	else
	168	{ /* 2 (x-1)/(x+1) */
	169	z = x - 0.5L;
	170	z -= 0.5L;
	171	y = 0.5L * x + 0.5L;
	172	}
	173	x = z / y;
	174	z = x * x;
	175	r = ((((R5 * z
	176	+ R4) * z
	177	+ R3) * z
	178	+ R2) * z
	179	+ R1) * z
	180	+ R0;
	181	s = (((((z
	182	+ S5) * z
	183	+ S4) * z
	184	+ S3) * z
	185	+ S2) * z
	186	+ S1) * z
	187	+ S0;
	188	z = x * (z * r / s);
	189	z = z + e * C2;
	190	z = z + x;
	191	z = z + e * C1;
	192	return (z);
	193	}
	194
	195
	196	/* Logarithm using log(1+x) = x - .5x^2 + x^3 P(x)/Q(x). */
	197
	198	if (x < sqrth)
	199	{
	200	e -= 1;
	201	if (e != 0)
	202	x = 2.0L * x - 1.0L; /* 2x - 1 */
	203	else
	204	x = xm1;
	205	}
	206	else
	207	{
	208	if (e != 0)
	209	x = x - 1.0L;
	210	else
	211	x = xm1;
	212	}
	213	z = x * x;
	214	r = (((((((((((P12 * x
	215	+ P11) * x
	216	+ P10) * x
	217	+ P9) * x
	218	+ P8) * x
	219	+ P7) * x
220	+ P6) * x
221	+ P5) * x
222	+ P4) * x
223	+ P3) * x
224	+ P2) * x
225	+ P1) * x
226	+ P0;
227	s = (((((((((((x
228	+ Q11) * x
229	+ Q10) * x
230	+ Q9) * x
231	+ Q8) * x
232	+ Q7) * x
233	+ Q6) * x
234	+ Q5) * x
235	+ Q4) * x
236	+ Q3) * x
237	+ Q2) * x
238	+ Q1) * x
239	+ Q0;
240	y = x * (z * r / s);
241	y = y + e * C2;
242	z = y - 0.5L * z;
243	z = z + x;
244	z = z + e * C1;
245	return (z);
246	}
247
248	weak_alias (__log1pl, log1pl)