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

/*							log10l.c
 *
 *	Common logarithm, 128-bit long double precision
 *
 *
 *
 * SYNOPSIS:
 *
 * long double x, y, log10l();
 *
 * y = log10l( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns the base 10 logarithm of x.
 *
 * The argument 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(x-1)/x+1),
 *
 *     log(x) = z + z^3 P(z)/Q(z).
 *
 *
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      0.5, 2.0     30000      2.3e-34     4.9e-35
 *    IEEE     exp(+-10000)  30000      1.0e-34     4.1e-35
 *
 * In the tests over the interval exp(+-10000), the logarithms
 * of the random arguments were uniformly distributed over
 * [-10000, +10000].
 *
 */

/*
   Cephes Math Library Release 2.2:  January, 1991
   Copyright 1984, 1991 by Stephen L. Moshier
   Adapted for glibc November, 2001

    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, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA

 */

#include "math.h"
#include "math_private.h"

/* Coefficients for ln(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
 * 1/sqrt(2) <= x < sqrt(2)
 * Theoretical peak relative error = 5.3e-37,
 * relative peak error spread = 2.3e-14
 */
static const long double P[13] =
{
  1.313572404063446165910279910527789794488E4L,
  7.771154681358524243729929227226708890930E4L,
  2.014652742082537582487669938141683759923E5L,
  3.007007295140399532324943111654767187848E5L,
  2.854829159639697837788887080758954924001E5L,
  1.797628303815655343403735250238293741397E5L,
  7.594356839258970405033155585486712125861E4L,
  2.128857716871515081352991964243375186031E4L,
  3.824952356185897735160588078446136783779E3L,
  4.114517881637811823002128927449878962058E2L,
  2.321125933898420063925789532045674660756E1L,
  4.998469661968096229986658302195402690910E-1L,
  1.538612243596254322971797716843006400388E-6L
};
static const long double Q[12] =
{
  3.940717212190338497730839731583397586124E4L,
  2.626900195321832660448791748036714883242E5L,
  7.777690340007566932935753241556479363645E5L,
  1.347518538384329112529391120390701166528E6L,
  1.514882452993549494932585972882995548426E6L,
  1.158019977462989115839826904108208787040E6L,
  6.132189329546557743179177159925690841200E5L,
  2.248234257620569139969141618556349415120E5L,
  5.605842085972455027590989944010492125825E4L,
  9.147150349299596453976674231612674085381E3L,
  9.104928120962988414618126155557301584078E2L,
  4.839208193348159620282142911143429644326E1L
/* 1.000000000000000000000000000000000000000E0L, */
};

/* 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 R[6] =
{
  1.418134209872192732479751274970992665513E5L,
 -8.977257995689735303686582344659576526998E4L,
  2.048819892795278657810231591630928516206E4L,
 -2.024301798136027039250415126250455056397E3L,
  8.057002716646055371965756206836056074715E1L,
 -8.828896441624934385266096344596648080902E-1L
};
static const long double S[6] =
{
  1.701761051846631278975701529965589676574E6L,
 -1.332535117259762928288745111081235577029E6L,
  4.001557694070773974936904547424676279307E5L,
 -5.748542087379434595104154610899551484314E4L,
  3.998526750980007367835804959888064681098E3L,
 -1.186359407982897997337150403816839480438E2L
/* 1.000000000000000000000000000000000000000E0L, */
};

static const long double
/* log10(2) */
L102A = 0.3125L,
L102B = -1.14700043360188047862611052755069732318101185E-2L,
/* log10(e) */
L10EA = 0.5L,
L10EB = -6.570551809674817234887108108339491770560299E-2L,
/* sqrt(2)/2 */
SQRTH = 7.071067811865475244008443621048490392848359E-1L;


/* Evaluate P[n] x^n  +  P[n-1] x^(n-1)  +  ...  +  P[0] */

static long double
neval (long double x, const long double *p, int n)
{
  long double y;

  p += n;
  y = *p--;
  do
    {
      y = y * x + *p--;
    }
  while (--n > 0);
  return y;
}


/* Evaluate x^n+1  +  P[n] x^(n)  +  P[n-1] x^(n-1)  +  ...  +  P[0] */

static long double
deval (long double x, const long double *p, int n)
{
  long double y;

  p += n;
  y = x + *p--;
  do
    {
      y = y * x + *p--;
    }
  while (--n > 0);
  return y;
}


long double
__ieee754_log10l (long double x)
{
  long double z;
  long double y;
  int e;
  int64_t hx, lx;

/* Test for domain */
  GET_LDOUBLE_WORDS64 (hx, lx, x);
  if (((hx & 0x7fffffffffffffffLL) | lx) == 0)
    return (-1.0L / (x - x));
  if (hx < 0)
    return (x - x) / (x - x);
  if (hx >= 0x7fff000000000000LL)
    return (x + x);

/* separate mantissa from exponent */

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


/* logarithm using log(x) = z + z**3 P(z)/Q(z),
 * 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;
      y = x * (z * neval (z, R, 5) / deval (z, S, 5));
      goto done;
    }


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

  if (x < SQRTH)
    {
      e -= 1;
      x = 2.0 * x - 1.0L;	/*  2x - 1  */
    }
  else
    {
      x = x - 1.0L;
    }
  z = x * x;
  y = x * (z * neval (x, P, 12) / deval (x, Q, 11));
  y = y - 0.5 * z;

done:

  /* Multiply log of fraction by log10(e)
   * and base 2 exponent by log10(2).
   */
  z = y * L10EB;
  z += x * L10EB;
  z += e * L102B;
  z += y * L10EA;
  z += x * L10EA;
  z += e * L102A;
  return (z);
}
strong_alias (__ieee754_log10l, __log10l_finite)
Commit	Line	Data
76321a25 AJ	1	/* log10l.c
	2	*
	3	* Common logarithm, 128-bit long double precision
	4	*
	5	*
	6	*
	7	* SYNOPSIS:
	8	*
	9	* long double x, y, log10l();
	10	*
	11	* y = log10l( x );
	12	*
	13	*
	14	*
	15	* DESCRIPTION:
	16	*
	17	* Returns the base 10 logarithm of x.
	18	*
	19	* The argument is separated into its exponent and fractional
	20	* parts. If the exponent is between -1 and +1, the logarithm
	21	* of the fraction is approximated by
	22	*
	23	* log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
	24	*
	25	* Otherwise, setting z = 2(x-1)/x+1),
	26	*
	27	* log(x) = z + z^3 P(z)/Q(z).
	28	*
	29	*
	30	*
	31	* ACCURACY:
	32	*
	33	* Relative error:
	34	* arithmetic domain # trials peak rms
	35	* IEEE 0.5, 2.0 30000 2.3e-34 4.9e-35
	36	* IEEE exp(+-10000) 30000 1.0e-34 4.1e-35
	37	*
	38	* In the tests over the interval exp(+-10000), the logarithms
	39	* of the random arguments were uniformly distributed over
	40	* [-10000, +10000].
	41	*
	42	*/
	43
	44	/*
	45	Cephes Math Library Release 2.2: January, 1991
	46	Copyright 1984, 1991 by Stephen L. Moshier
	47	Adapted for glibc November, 2001
cc7375ce RM	48
	49	This library is free software; you can redistribute it and/or
	50	modify it under the terms of the GNU Lesser General Public
	51	License as published by the Free Software Foundation; either
	52	version 2.1 of the License, or (at your option) any later version.
	53
	54	This library is distributed in the hope that it will be useful,
	55	but WITHOUT ANY WARRANTY; without even the implied warranty of
	56	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
	57	Lesser General Public License for more details.
	58
	59	You should have received a copy of the GNU Lesser General Public
	60	License along with this library; if not, write to the Free Software
0ac5ae23	61	Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
cc7375ce	62
76321a25 AJ	63	*/
	64
	65	#include "math.h"
	66	#include "math_private.h"
	67
	68	/* Coefficients for ln(1+x) = x - x2/2 + x3 P(x)/Q(x)
	69	* 1/sqrt(2) <= x < sqrt(2)
	70	* Theoretical peak relative error = 5.3e-37,
	71	* relative peak error spread = 2.3e-14
	72	*/
	73	static const long double P[13] =
	74	{
	75	1.313572404063446165910279910527789794488E4L,
	76	7.771154681358524243729929227226708890930E4L,
	77	2.014652742082537582487669938141683759923E5L,
	78	3.007007295140399532324943111654767187848E5L,
	79	2.854829159639697837788887080758954924001E5L,
	80	1.797628303815655343403735250238293741397E5L,
	81	7.594356839258970405033155585486712125861E4L,
	82	2.128857716871515081352991964243375186031E4L,
	83	3.824952356185897735160588078446136783779E3L,
	84	4.114517881637811823002128927449878962058E2L,
	85	2.321125933898420063925789532045674660756E1L,
	86	4.998469661968096229986658302195402690910E-1L,
	87	1.538612243596254322971797716843006400388E-6L
	88	};
	89	static const long double Q[12] =
	90	{
	91	3.940717212190338497730839731583397586124E4L,
	92	2.626900195321832660448791748036714883242E5L,
	93	7.777690340007566932935753241556479363645E5L,
	94	1.347518538384329112529391120390701166528E6L,
	95	1.514882452993549494932585972882995548426E6L,
	96	1.158019977462989115839826904108208787040E6L,
	97	6.132189329546557743179177159925690841200E5L,
	98	2.248234257620569139969141618556349415120E5L,
	99	5.605842085972455027590989944010492125825E4L,
	100	9.147150349299596453976674231612674085381E3L,
	101	9.104928120962988414618126155557301584078E2L,
	102	4.839208193348159620282142911143429644326E1L
	103	/* 1.000000000000000000000000000000000000000E0L, */
	104	};
	105
	106	/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
	107	* where z = 2(x-1)/(x+1)
	108	* 1/sqrt(2) <= x < sqrt(2)
	109	* Theoretical peak relative error = 1.1e-35,
	110	* relative peak error spread 1.1e-9
	111	*/
	112	static const long double R[6] =
	113	{
	114	1.418134209872192732479751274970992665513E5L,
	115	-8.977257995689735303686582344659576526998E4L,
	116	2.048819892795278657810231591630928516206E4L,
	117	-2.024301798136027039250415126250455056397E3L,
	118	8.057002716646055371965756206836056074715E1L,
	119	-8.828896441624934385266096344596648080902E-1L
	120	};
	121	static const long double S[6] =
	122	{
	123	1.701761051846631278975701529965589676574E6L,
	124	-1.332535117259762928288745111081235577029E6L,
	125	4.001557694070773974936904547424676279307E5L,
	126	-5.748542087379434595104154610899551484314E4L,
127	3.998526750980007367835804959888064681098E3L,
128	-1.186359407982897997337150403816839480438E2L
129	/* 1.000000000000000000000000000000000000000E0L, */
130	};
131
132	static const long double
133	/* log10(2) */
134	L102A = 0.3125L,
9992bc08	135	L102B = -1.14700043360188047862611052755069732318101185E-2L,
76321a25 AJ	136	/* log10(e) */
	137	L10EA = 0.5L,
	138	L10EB = -6.570551809674817234887108108339491770560299E-2L,
	139	/* sqrt(2)/2 */
	140	SQRTH = 7.071067811865475244008443621048490392848359E-1L;
	141
	142
	143
	144	/* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */
	145
	146	static long double
	147	neval (long double x, const long double *p, int n)
	148	{
	149	long double y;
	150
	151	p += n;
	152	y = *p--;
	153	do
	154	{
	155	y = y * x + *p--;
	156	}
	157	while (--n > 0);
	158	return y;
	159	}
	160
	161
	162	/* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */
	163
	164	static long double
	165	deval (long double x, const long double *p, int n)
	166	{
	167	long double y;
	168
	169	p += n;
	170	y = x + *p--;
	171	do
	172	{
	173	y = y * x + *p--;
	174	}
	175	while (--n > 0);
	176	return y;
	177	}
	178
	179
	180
	181	long double
0ac5ae23	182	__ieee754_log10l (long double x)
76321a25 AJ	183	{
	184	long double z;
	185	long double y;
	186	int e;
52e1b618	187	int64_t hx, lx;
76321a25 AJ	188
76321a25 AJ	189	/* Test for domain */
52e1b618 UD	190	GET_LDOUBLE_WORDS64 (hx, lx, x);
	191	if (((hx & 0x7fffffffffffffffLL) \| lx) == 0)
	192	return (-1.0L / (x - x));
	193	if (hx < 0)
	194	return (x - x) / (x - x);
	195	if (hx >= 0x7fff000000000000LL)
76321a25 AJ	196	return (x + x);
	197
	198	/* separate mantissa from exponent */
	199
	200	/* Note, frexp is used so that denormal numbers
	201	* will be handled properly.
	202	*/
	203	x = __frexpl (x, &e);
	204
	205
	206	/* logarithm using log(x) = z + z**3 P(z)/Q(z),
	207	* where z = 2(x-1)/x+1)
	208	*/
	209	if ((e > 2) \|\| (e < -2))
	210	{
	211	if (x < SQRTH)
	212	{ /* 2( 2x-1 )/( 2x+1 ) */
	213	e -= 1;
	214	z = x - 0.5L;
	215	y = 0.5L * z + 0.5L;
	216	}
	217	else
	218	{ /* 2 (x-1)/(x+1) */
	219	z = x - 0.5L;
	220	z -= 0.5L;
	221	y = 0.5L * x + 0.5L;
	222	}
	223	x = z / y;
	224	z = x * x;
	225	y = x * (z * neval (z, R, 5) / deval (z, S, 5));
	226	goto done;
	227	}
	228
	229
	230	/* logarithm using log(1+x) = x - .5x2 + x3 P(x)/Q(x) */
	231
	232	if (x < SQRTH)
	233	{
	234	e -= 1;
	235	x = 2.0 * x - 1.0L; /* 2x - 1 */
	236	}
	237	else
	238	{
	239	x = x - 1.0L;
	240	}
	241	z = x * x;
	242	y = x * (z * neval (x, P, 12) / deval (x, Q, 11));
	243	y = y - 0.5 * z;
	244
	245	done:
	246
	247	/* Multiply log of fraction by log10(e)
	248	* and base 2 exponent by log10(2).
	249	*/
	250	z = y * L10EB;
	251	z += x * L10EB;
	252	z += e * L102B;
	253	z += y * L10EA;
	254	z += x * L10EA;
	255	z += e * L102A;
	256	return (z);
	257	}
0ac5ae23	258	strong_alias (__ieee754_log10l, __log10l_finite)