Python code coverage for Modules/_math.c

#	count	content
1	n/a	/* Definitions of some C99 math library functions, for those platforms
2	n/a	that don't implement these functions already. */
3	n/a
4	n/a	#include "Python.h"
5	n/a	#include <float.h>
6	n/a	#include "_math.h"
7	n/a
8	n/a	/* The following copyright notice applies to the original
9	n/a	implementations of acosh, asinh and atanh. */
10	n/a
11	n/a	/*
12	n/a	* ====================================================
13	n/a	* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
14	n/a	*
15	n/a	* Developed at SunPro, a Sun Microsystems, Inc. business.
16	n/a	* Permission to use, copy, modify, and distribute this
17	n/a	* software is freely granted, provided that this notice
18	n/a	* is preserved.
19	n/a	* ====================================================
20	n/a	*/
21	n/a
22	n/a	#if !defined(HAVE_ACOSH) \|\| !defined(HAVE_ASINH)
23	n/a	static const double ln2 = 6.93147180559945286227E-01;
24	n/a	static const double two_pow_p28 = 268435456.0; /* 2*28 /
25	n/a	#endif
26	n/a	#if !defined(HAVE_ASINH) \|\| !defined(HAVE_ATANH)
27	n/a	static const double two_pow_m28 = 3.7252902984619141E-09; /* 2*-28 /
28	n/a	#endif
29	n/a	#if !defined(HAVE_ATANH) && !defined(Py_NAN)
30	n/a	static const double zero = 0.0;
31	n/a	#endif
32	n/a
33	n/a
34	n/a	#ifndef HAVE_ACOSH
35	n/a	/* acosh(x)
36	n/a	* Method :
37	n/a	* Based on
38	n/a	* acosh(x) = log [ x + sqrt(x*x-1) ]
39	n/a	* we have
40	n/a	* acosh(x) := log(x)+ln2, if x is large; else
41	n/a	* acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else
42	n/a	* acosh(x) := log1p(t+sqrt(2.0t+tt)); where t=x-1.
43	n/a	*
44	n/a	* Special cases:
45	n/a	* acosh(x) is NaN with signal if x<1.
46	n/a	* acosh(NaN) is NaN without signal.
47	n/a	*/
48	n/a
49	n/a	double
50	n/a	_Py_acosh(double x)
51	n/a	{
52	n/a	if (Py_IS_NAN(x)) {
53	n/a	return x+x;
54	n/a	}
55	n/a	if (x < 1.) { /* x < 1; return a signaling NaN */
56	n/a	errno = EDOM;
57	n/a	#ifdef Py_NAN
58	n/a	return Py_NAN;
59	n/a	#else
60	n/a	return (x-x)/(x-x);
61	n/a	#endif
62	n/a	}
63	n/a	else if (x >= two_pow_p28) { /* x > 2*28 /
64	n/a	if (Py_IS_INFINITY(x)) {
65	n/a	return x+x;
66	n/a	}
67	n/a	else {
68	n/a	return log(x) + ln2; /* acosh(huge)=log(2x) */
69	n/a	}
70	n/a	}
71	n/a	else if (x == 1.) {
72	n/a	return 0.0; /* acosh(1) = 0 */
73	n/a	}
74	n/a	else if (x > 2.) { /* 2 < x < 2*28 /
75	n/a	double t = x * x;
76	n/a	return log(2.0 * x - 1.0 / (x + sqrt(t - 1.0)));
77	n/a	}
78	n/a	else { /* 1 < x <= 2 */
79	n/a	double t = x - 1.0;
80	n/a	return m_log1p(t + sqrt(2.0 * t + t * t));
81	n/a	}
82	n/a	}
83	n/a	#endif /* HAVE_ACOSH */
84	n/a
85	n/a
86	n/a	#ifndef HAVE_ASINH
87	n/a	/* asinh(x)
88	n/a	* Method :
89	n/a	* Based on
90	n/a	* asinh(x) = sign(x) * log [ \|x\| + sqrt(x*x+1) ]
91	n/a	* we have
92	n/a	* asinh(x) := x if 1+x*x=1,
93	n/a	* := sign(x)*(log(x)+ln2)) for large \|x\|, else
94	n/a	* := sign(x)log(2\|x\|+1/(\|x\|+sqrt(xx+1))) if\|x\|>2, else
95	n/a	* := sign(x)*log1p(\|x\| + x^2/(1 + sqrt(1+x^2)))
96	n/a	*/
97	n/a
98	n/a	double
99	n/a	_Py_asinh(double x)
100	n/a	{
101	n/a	double w;
102	n/a	double absx = fabs(x);
103	n/a
104	n/a	if (Py_IS_NAN(x) \|\| Py_IS_INFINITY(x)) {
105	n/a	return x+x;
106	n/a	}
107	n/a	if (absx < two_pow_m28) { /* \|x\| < 2*-28 /
108	n/a	return x; /* return x inexact except 0 */
109	n/a	}
110	n/a	if (absx > two_pow_p28) { /* \|x\| > 2*28 /
111	n/a	w = log(absx) + ln2;
112	n/a	}
113	n/a	else if (absx > 2.0) { /* 2 < \|x\| < 2*28 /
114	n/a	w = log(2.0 * absx + 1.0 / (sqrt(x * x + 1.0) + absx));
115	n/a	}
116	n/a	else { /* 2*-28 <= \|x\| < 2= /
117	n/a	double t = x*x;
118	n/a	w = m_log1p(absx + t / (1.0 + sqrt(1.0 + t)));
119	n/a	}
120	n/a	return copysign(w, x);
121	n/a
122	n/a	}
123	n/a	#endif /* HAVE_ASINH */
124	n/a
125	n/a
126	n/a	#ifndef HAVE_ATANH
127	n/a	/* atanh(x)
128	n/a	* Method :
129	n/a	* 1.Reduced x to positive by atanh(-x) = -atanh(x)
130	n/a	* 2.For x>=0.5
131	n/a	* 1 2x x
132	n/a	* atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * -------)
133	n/a	* 2 1 - x 1 - x
134	n/a	*
135	n/a	* For x<0.5
136	n/a	* atanh(x) = 0.5log1p(2x+2xx/(1-x))
137	n/a	*
138	n/a	* Special cases:
139	n/a	* atanh(x) is NaN if \|x\| >= 1 with signal;
140	n/a	* atanh(NaN) is that NaN with no signal;
141	n/a	*
142	n/a	*/
143	n/a
144	n/a	double
145	n/a	_Py_atanh(double x)
146	n/a	{
147	n/a	double absx;
148	n/a	double t;
149	n/a
150	n/a	if (Py_IS_NAN(x)) {
151	n/a	return x+x;
152	n/a	}
153	n/a	absx = fabs(x);
154	n/a	if (absx >= 1.) { /* \|x\| >= 1 */
155	n/a	errno = EDOM;
156	n/a	#ifdef Py_NAN
157	n/a	return Py_NAN;
158	n/a	#else
159	n/a	return x / zero;
160	n/a	#endif
161	n/a	}
162	n/a	if (absx < two_pow_m28) { /* \|x\| < 2*-28 /
163	n/a	return x;
164	n/a	}
165	n/a	if (absx < 0.5) { /* \|x\| < 0.5 */
166	n/a	t = absx+absx;
167	n/a	t = 0.5 * m_log1p(t + t*absx / (1.0 - absx));
168	n/a	}
169	n/a	else { /* 0.5 <= \|x\| <= 1.0 */
170	n/a	t = 0.5 * m_log1p((absx + absx) / (1.0 - absx));
171	n/a	}
172	n/a	return copysign(t, x);
173	n/a	}
174	n/a	#endif /* HAVE_ATANH */
175	n/a
176	n/a
177	n/a	#ifndef HAVE_EXPM1
178	n/a	/* Mathematically, expm1(x) = exp(x) - 1. The expm1 function is designed
179	n/a	to avoid the significant loss of precision that arises from direct
180	n/a	evaluation of the expression exp(x) - 1, for x near 0. */
181	n/a
182	n/a	double
183	n/a	_Py_expm1(double x)
184	n/a	{
185	n/a	/* For abs(x) >= log(2), it's safe to evaluate exp(x) - 1 directly; this
186	n/a	also works fine for infinities and nans.
187	n/a
188	n/a	For smaller x, we can use a method due to Kahan that achieves close to
189	n/a	full accuracy.
190	n/a	*/
191	n/a
192	n/a	if (fabs(x) < 0.7) {
193	n/a	double u;
194	n/a	u = exp(x);
195	n/a	if (u == 1.0)
196	n/a	return x;
197	n/a	else
198	n/a	return (u - 1.0) * x / log(u);
199	n/a	}
200	n/a	else
201	n/a	return exp(x) - 1.0;
202	n/a	}
203	n/a	#endif /* HAVE_EXPM1 */
204	n/a
205	n/a
206	n/a	/* log1p(x) = log(1+x). The log1p function is designed to avoid the
207	n/a	significant loss of precision that arises from direct evaluation when x is
208	n/a	small. */
209	n/a
210	n/a	double
211	n/a	_Py_log1p(double x)
212	n/a	{
213	n/a	#ifdef HAVE_LOG1P
214	n/a	/* Some platforms supply a log1p function but don't respect the sign of
215	n/a	zero: log1p(-0.0) gives 0.0 instead of the correct result of -0.0.
216	n/a
217	n/a	To save fiddling with configure tests and platform checks, we handle the
218	n/a	special case of zero input directly on all platforms.
219	n/a	*/
220	n/a	if (x == 0.0) {
221	n/a	return x;
222	n/a	}
223	n/a	else {
224	n/a	return log1p(x);
225	n/a	}
226	n/a	#else
227	n/a	/* For x small, we use the following approach. Let y be the nearest float
228	n/a	to 1+x, then
229	n/a
230	n/a	1+x = y * (1 - (y-1-x)/y)
231	n/a
232	n/a	so log(1+x) = log(y) + log(1-(y-1-x)/y). Since (y-1-x)/y is tiny, the
233	n/a	second term is well approximated by (y-1-x)/y. If abs(x) >=
234	n/a	DBL_EPSILON/2 or the rounding-mode is some form of round-to-nearest
235	n/a	then y-1-x will be exactly representable, and is computed exactly by
236	n/a	(y-1)-x.
237	n/a
238	n/a	If abs(x) < DBL_EPSILON/2 and the rounding mode is not known to be
239	n/a	round-to-nearest then this method is slightly dangerous: 1+x could be
240	n/a	rounded up to 1+DBL_EPSILON instead of down to 1, and in that case
241	n/a	y-1-x will not be exactly representable any more and the result can be
242	n/a	off by many ulps. But this is easily fixed: for a floating-point
243	n/a	number \|x\| < DBL_EPSILON/2., the closest floating-point number to
244	n/a	log(1+x) is exactly x.
245	n/a	*/
246	n/a
247	n/a	double y;
248	n/a	if (fabs(x) < DBL_EPSILON / 2.) {
249	n/a	return x;
250	n/a	}
251	n/a	else if (-0.5 <= x && x <= 1.) {
252	n/a	/* WARNING: it's possible that an overeager compiler
253	n/a	will incorrectly optimize the following two lines
254	n/a	to the equivalent of "return log(1.+x)". If this
255	n/a	happens, then results from log1p will be inaccurate
256	n/a	for small x. */
257	n/a	y = 1.+x;
258	n/a	return log(y) - ((y - 1.) - x) / y;
259	n/a	}
260	n/a	else {
261	n/a	/* NaNs and infinities should end up here */
262	n/a	return log(1.+x);
263	n/a	}
264	n/a	#endif /* ifdef HAVE_LOG1P */
265	n/a	}
266	n/a