[thirdparty/glibc.git] / sysdeps / ieee754 / ldbl-128ibm / 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 (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 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);
}
Commit	Line	Data
0ecb606c JJ	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
	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
	61	Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
	62
	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,
135	L102B = -1.14700043360188047862611052755069732318101185E-2L,
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
182	__ieee754_log10l (x)
183	long double x;
184	{
185	long double z;
186	long double y;
187	int e;
188	int64_t hx, lx;
189
190	/* Test for domain */
191	GET_LDOUBLE_WORDS64 (hx, lx, x);
192	if (((hx & 0x7fffffffffffffffLL) \| (lx & 0x7fffffffffffffffLL)) == 0)
193	return (-1.0L / (x - x));
194	if (hx < 0)
195	return (x - x) / (x - x);
196	if (hx >= 0x7ff0000000000000LL)
197	return (x + x);
198
199	/* separate mantissa from exponent */
200
201	/* Note, frexp is used so that denormal numbers
202	* will be handled properly.
203	*/
204	x = __frexpl (x, &e);
205
206
207	/* logarithm using log(x) = z + z**3 P(z)/Q(z),
208	* where z = 2(x-1)/x+1)
209	*/
210	if ((e > 2) \|\| (e < -2))
211	{
212	if (x < SQRTH)
213	{ /* 2( 2x-1 )/( 2x+1 ) */
214	e -= 1;
215	z = x - 0.5L;
216	y = 0.5L * z + 0.5L;
217	}
218	else
219	{ /* 2 (x-1)/(x+1) */
220	z = x - 0.5L;
221	z -= 0.5L;
222	y = 0.5L * x + 0.5L;
223	}
224	x = z / y;
225	z = x * x;
226	y = x * (z * neval (z, R, 5) / deval (z, S, 5));
227	goto done;
228	}
229
230
231	/* logarithm using log(1+x) = x - .5x2 + x3 P(x)/Q(x) */
232
233	if (x < SQRTH)
234	{
235	e -= 1;
236	x = 2.0 * x - 1.0L; /* 2x - 1 */
237	}
238	else
239	{
240	x = x - 1.0L;
241	}
242	z = x * x;
243	y = x * (z * neval (x, P, 12) / deval (x, Q, 11));
244	y = y - 0.5 * z;
245
246	done:
247
248	/* Multiply log of fraction by log10(e)
249	* and base 2 exponent by log10(2).
250	*/
251	z = y * L10EB;
252	z += x * L10EB;
253	z += e * L102B;
254	z += y * L10EA;
255	z += x * L10EA;
256	z += e * L102A;
257	return (z);
258	}