[glibc.git] / sysdeps / ieee754 / ldbl-128ibm / e_log2l.c

/*                                                      log2l.c
 *      Base 2 logarithm, 128-bit long double precision
 *
 *
 *
 * SYNOPSIS:
 *
 * long double x, y, log2l();
 *
 * y = log2l( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns the base 2 logarithm of x.
 *
 * The argument is separated into its exponent and fractional
 * parts.  If the exponent is between -1 and +1, the (natural)
 * 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     100,000    2.6e-34     4.9e-35
 *    IEEE     exp(+-10000)  100,000    9.6e-35     4.0e-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, see <http://www.gnu.org/licenses/>.
 */

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