Python code coverage for Modules/mathmodule.c

#	count	content
1	n/a	/* Math module -- standard C math library functions, pi and e */
2	n/a
3	n/a	/* Here are some comments from Tim Peters, extracted from the
4	n/a	discussion attached to http://bugs.python.org/issue1640. They
5	n/a	describe the general aims of the math module with respect to
6	n/a	special values, IEEE-754 floating-point exceptions, and Python
7	n/a	exceptions.
8	n/a
9	n/a	These are the "spirit of 754" rules:
10	n/a
11	n/a	1. If the mathematical result is a real number, but of magnitude too
12	n/a	large to approximate by a machine float, overflow is signaled and the
13	n/a	result is an infinity (with the appropriate sign).
14	n/a
15	n/a	2. If the mathematical result is a real number, but of magnitude too
16	n/a	small to approximate by a machine float, underflow is signaled and the
17	n/a	result is a zero (with the appropriate sign).
18	n/a
19	n/a	3. At a singularity (a value x such that the limit of f(y) as y
20	n/a	approaches x exists and is an infinity), "divide by zero" is signaled
21	n/a	and the result is an infinity (with the appropriate sign). This is
22	n/a	complicated a little by that the left-side and right-side limits may
23	n/a	not be the same; e.g., 1/x approaches +inf or -inf as x approaches 0
24	n/a	from the positive or negative directions. In that specific case, the
25	n/a	sign of the zero determines the result of 1/0.
26	n/a
27	n/a	4. At a point where a function has no defined result in the extended
28	n/a	reals (i.e., the reals plus an infinity or two), invalid operation is
29	n/a	signaled and a NaN is returned.
30	n/a
31	n/a	And these are what Python has historically /tried/ to do (but not
32	n/a	always successfully, as platform libm behavior varies a lot):
33	n/a
34	n/a	For #1, raise OverflowError.
35	n/a
36	n/a	For #2, return a zero (with the appropriate sign if that happens by
37	n/a	accident ;-)).
38	n/a
39	n/a	For #3 and #4, raise ValueError. It may have made sense to raise
40	n/a	Python's ZeroDivisionError in #3, but historically that's only been
41	n/a	raised for division by zero and mod by zero.
42	n/a
43	n/a	*/
44	n/a
45	n/a	/*
46	n/a	In general, on an IEEE-754 platform the aim is to follow the C99
47	n/a	standard, including Annex 'F', whenever possible. Where the
48	n/a	standard recommends raising the 'divide-by-zero' or 'invalid'
49	n/a	floating-point exceptions, Python should raise a ValueError. Where
50	n/a	the standard recommends raising 'overflow', Python should raise an
51	n/a	OverflowError. In all other circumstances a value should be
52	n/a	returned.
53	n/a	*/
54	n/a
55	n/a	#include "Python.h"
56	n/a	#include "_math.h"
57	n/a
58	n/a	#include "clinic/mathmodule.c.h"
59	n/a
60	n/a	/*[clinic input]
61	n/a	module math
62	n/a	[clinic start generated code]*/
63	n/a	/[clinic end generated code: output=da39a3ee5e6b4b0d input=76bc7002685dd942]/
64	n/a
65	n/a
66	n/a	/*
67	n/a	sin(pi*x), giving accurate results for all finite x (especially x
68	n/a	integral or close to an integer). This is here for use in the
69	n/a	reflection formula for the gamma function. It conforms to IEEE
70	n/a	754-2008 for finite arguments, but not for infinities or nans.
71	n/a	*/
72	n/a
73	n/a	static const double pi = 3.141592653589793238462643383279502884197;
74	n/a	static const double sqrtpi = 1.772453850905516027298167483341145182798;
75	n/a	static const double logpi = 1.144729885849400174143427351353058711647;
76	n/a
77	n/a	static double
78	n/a	sinpi(double x)
79	n/a	{
80	n/a	double y, r;
81	n/a	int n;
82	n/a	/* this function should only ever be called for finite arguments */
83	n/a	assert(Py_IS_FINITE(x));
84	n/a	y = fmod(fabs(x), 2.0);
85	n/a	n = (int)round(2.0*y);
86	n/a	assert(0 <= n && n <= 4);
87	n/a	switch (n) {
88	n/a	case 0:
89	n/a	r = sin(pi*y);
90	n/a	break;
91	n/a	case 1:
92	n/a	r = cos(pi*(y-0.5));
93	n/a	break;
94	n/a	case 2:
95	n/a	/* N.B. -sin(pi(y-1.0)) is not* equivalent: it would give
96	n/a	-0.0 instead of 0.0 when y == 1.0. */
97	n/a	r = sin(pi*(1.0-y));
98	n/a	break;
99	n/a	case 3:
100	n/a	r = -cos(pi*(y-1.5));
101	n/a	break;
102	n/a	case 4:
103	n/a	r = sin(pi*(y-2.0));
104	n/a	break;
105	n/a	default:
106	n/a	assert(0); /* should never get here */
107	n/a	r = -1.23e200; /* silence gcc warning */
108	n/a	}
109	n/a	return copysign(1.0, x)*r;
110	n/a	}
111	n/a
112	n/a	/* Implementation of the real gamma function. In extensive but non-exhaustive
113	n/a	random tests, this function proved accurate to within <= 10 ulps across the
114	n/a	entire float domain. Note that accuracy may depend on the quality of the
115	n/a	system math functions, the pow function in particular. Special cases
116	n/a	follow C99 annex F. The parameters and method are tailored to platforms
117	n/a	whose double format is the IEEE 754 binary64 format.
118	n/a
119	n/a	Method: for x > 0.0 we use the Lanczos approximation with parameters N=13
120	n/a	and g=6.024680040776729583740234375; these parameters are amongst those
121	n/a	used by the Boost library. Following Boost (again), we re-express the
122	n/a	Lanczos sum as a rational function, and compute it that way. The
123	n/a	coefficients below were computed independently using MPFR, and have been
124	n/a	double-checked against the coefficients in the Boost source code.
125	n/a
126	n/a	For x < 0.0 we use the reflection formula.
127	n/a
128	n/a	There's one minor tweak that deserves explanation: Lanczos' formula for
129	n/a	Gamma(x) involves computing pow(x+g-0.5, x-0.5) / exp(x+g-0.5). For many x
130	n/a	values, x+g-0.5 can be represented exactly. However, in cases where it
131	n/a	can't be represented exactly the small error in x+g-0.5 can be magnified
132	n/a	significantly by the pow and exp calls, especially for large x. A cheap
133	n/a	correction is to multiply by (1 + e*g/(x+g-0.5)), where e is the error
134	n/a	involved in the computation of x+g-0.5 (that is, e = computed value of
135	n/a	x+g-0.5 - exact value of x+g-0.5). Here's the proof:
136	n/a
137	n/a	Correction factor
138	n/a	-----------------
139	n/a	Write x+g-0.5 = y-e, where y is exactly representable as an IEEE 754
140	n/a	double, and e is tiny. Then:
141	n/a
142	n/a	pow(x+g-0.5,x-0.5)/exp(x+g-0.5) = pow(y-e, x-0.5)/exp(y-e)
143	n/a	= pow(y, x-0.5)/exp(y) * C,
144	n/a
145	n/a	where the correction_factor C is given by
146	n/a
147	n/a	C = pow(1-e/y, x-0.5) * exp(e)
148	n/a
149	n/a	Since e is tiny, pow(1-e/y, x-0.5) ~ 1-(x-0.5)*e/y, and exp(x) ~ 1+e, so:
150	n/a
151	n/a	C ~ (1-(x-0.5)e/y) (1+e) ~ 1 + e*(y-(x-0.5))/y
152	n/a
153	n/a	But y-(x-0.5) = g+e, and g+e ~ g. So we get C ~ 1 + e*g/y, and
154	n/a
155	n/a	pow(x+g-0.5,x-0.5)/exp(x+g-0.5) ~ pow(y, x-0.5)/exp(y) * (1 + e*g/y),
156	n/a
157	n/a	Note that for accuracy, when computing r*C it's better to do
158	n/a
159	n/a	r + eg/yr;
160	n/a
161	n/a	than
162	n/a
163	n/a	r * (1 + e*g/y);
164	n/a
165	n/a	since the addition in the latter throws away most of the bits of
166	n/a	information in e*g/y.
167	n/a	*/
168	n/a
169	n/a	#define LANCZOS_N 13
170	n/a	static const double lanczos_g = 6.024680040776729583740234375;
171	n/a	static const double lanczos_g_minus_half = 5.524680040776729583740234375;
172	n/a	static const double lanczos_num_coeffs[LANCZOS_N] = {
173	n/a	23531376880.410759688572007674451636754734846804940,
174	n/a	42919803642.649098768957899047001988850926355848959,
175	n/a	35711959237.355668049440185451547166705960488635843,
176	n/a	17921034426.037209699919755754458931112671403265390,
177	n/a	6039542586.3520280050642916443072979210699388420708,
178	n/a	1439720407.3117216736632230727949123939715485786772,
179	n/a	248874557.86205415651146038641322942321632125127801,
180	n/a	31426415.585400194380614231628318205362874684987640,
181	n/a	2876370.6289353724412254090516208496135991145378768,
182	n/a	186056.26539522349504029498971604569928220784236328,
183	n/a	8071.6720023658162106380029022722506138218516325024,
184	n/a	210.82427775157934587250973392071336271166969580291,
185	n/a	2.5066282746310002701649081771338373386264310793408
186	n/a	};
187	n/a
188	n/a	/* denominator is x(x+1)...(x+LANCZOS_N-2) /
189	n/a	static const double lanczos_den_coeffs[LANCZOS_N] = {
190	n/a	0.0, 39916800.0, 120543840.0, 150917976.0, 105258076.0, 45995730.0,
191	n/a	13339535.0, 2637558.0, 357423.0, 32670.0, 1925.0, 66.0, 1.0};
192	n/a
193	n/a	/* gamma values for small positive integers, 1 though NGAMMA_INTEGRAL */
194	n/a	#define NGAMMA_INTEGRAL 23
195	n/a	static const double gamma_integral[NGAMMA_INTEGRAL] = {
196	n/a	1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0, 362880.0,
197	n/a	3628800.0, 39916800.0, 479001600.0, 6227020800.0, 87178291200.0,
198	n/a	1307674368000.0, 20922789888000.0, 355687428096000.0,
199	n/a	6402373705728000.0, 121645100408832000.0, 2432902008176640000.0,
200	n/a	51090942171709440000.0, 1124000727777607680000.0,
201	n/a	};
202	n/a
203	n/a	/* Lanczos' sum L_g(x), for positive x */
204	n/a
205	n/a	static double
206	n/a	lanczos_sum(double x)
207	n/a	{
208	n/a	double num = 0.0, den = 0.0;
209	n/a	int i;
210	n/a	assert(x > 0.0);
211	n/a	/* evaluate the rational function lanczos_sum(x). For large
212	n/a	x, the obvious algorithm risks overflow, so we instead
213	n/a	rescale the denominator and numerator of the rational
214	n/a	function by x**(1-LANCZOS_N) and treat this as a
215	n/a	rational function in 1/x. This also reduces the error for
216	n/a	larger x values. The choice of cutoff point (5.0 below) is
217	n/a	somewhat arbitrary; in tests, smaller cutoff values than
218	n/a	this resulted in lower accuracy. */
219	n/a	if (x < 5.0) {
220	n/a	for (i = LANCZOS_N; --i >= 0; ) {
221	n/a	num = num * x + lanczos_num_coeffs[i];
222	n/a	den = den * x + lanczos_den_coeffs[i];
223	n/a	}
224	n/a	}
225	n/a	else {
226	n/a	for (i = 0; i < LANCZOS_N; i++) {
227	n/a	num = num / x + lanczos_num_coeffs[i];
228	n/a	den = den / x + lanczos_den_coeffs[i];
229	n/a	}
230	n/a	}
231	n/a	return num/den;
232	n/a	}
233	n/a
234	n/a	/* Constant for +infinity, generated in the same way as float('inf'). */
235	n/a
236	n/a	static double
237	n/a	m_inf(void)
238	n/a	{
239	n/a	#ifndef PY_NO_SHORT_FLOAT_REPR
240	n/a	return _Py_dg_infinity(0);
241	n/a	#else
242	n/a	return Py_HUGE_VAL;
243	n/a	#endif
244	n/a	}
245	n/a
246	n/a	/* Constant nan value, generated in the same way as float('nan'). */
247	n/a	/* We don't currently assume that Py_NAN is defined everywhere. */
248	n/a
249	n/a	#if !defined(PY_NO_SHORT_FLOAT_REPR) \|\| defined(Py_NAN)
250	n/a
251	n/a	static double
252	n/a	m_nan(void)
253	n/a	{
254	n/a	#ifndef PY_NO_SHORT_FLOAT_REPR
255	n/a	return _Py_dg_stdnan(0);
256	n/a	#else
257	n/a	return Py_NAN;
258	n/a	#endif
259	n/a	}
260	n/a
261	n/a	#endif
262	n/a
263	n/a	static double
264	n/a	m_tgamma(double x)
265	n/a	{
266	n/a	double absx, r, y, z, sqrtpow;
267	n/a
268	n/a	/* special cases */
269	n/a	if (!Py_IS_FINITE(x)) {
270	n/a	if (Py_IS_NAN(x) \|\| x > 0.0)
271	n/a	return x; /* tgamma(nan) = nan, tgamma(inf) = inf */
272	n/a	else {
273	n/a	errno = EDOM;
274	n/a	return Py_NAN; /* tgamma(-inf) = nan, invalid */
275	n/a	}
276	n/a	}
277	n/a	if (x == 0.0) {
278	n/a	errno = EDOM;
279	n/a	/* tgamma(+-0.0) = +-inf, divide-by-zero */
280	n/a	return copysign(Py_HUGE_VAL, x);
281	n/a	}
282	n/a
283	n/a	/* integer arguments */
284	n/a	if (x == floor(x)) {
285	n/a	if (x < 0.0) {
286	n/a	errno = EDOM; /* tgamma(n) = nan, invalid for */
287	n/a	return Py_NAN; /* negative integers n */
288	n/a	}
289	n/a	if (x <= NGAMMA_INTEGRAL)
290	n/a	return gamma_integral[(int)x - 1];
291	n/a	}
292	n/a	absx = fabs(x);
293	n/a
294	n/a	/* tiny arguments: tgamma(x) ~ 1/x for x near 0 */
295	n/a	if (absx < 1e-20) {
296	n/a	r = 1.0/x;
297	n/a	if (Py_IS_INFINITY(r))
298	n/a	errno = ERANGE;
299	n/a	return r;
300	n/a	}
301	n/a
302	n/a	/* large arguments: assuming IEEE 754 doubles, tgamma(x) overflows for
303	n/a	x > 200, and underflows to +-0.0 for x < -200, not a negative
304	n/a	integer. */
305	n/a	if (absx > 200.0) {
306	n/a	if (x < 0.0) {
307	n/a	return 0.0/sinpi(x);
308	n/a	}
309	n/a	else {
310	n/a	errno = ERANGE;
311	n/a	return Py_HUGE_VAL;
312	n/a	}
313	n/a	}
314	n/a
315	n/a	y = absx + lanczos_g_minus_half;
316	n/a	/* compute error in sum */
317	n/a	if (absx > lanczos_g_minus_half) {
318	n/a	/* note: the correction can be foiled by an optimizing
319	n/a	compiler that (incorrectly) thinks that an expression like
320	n/a	a + b - a - b can be optimized to 0.0. This shouldn't
321	n/a	happen in a standards-conforming compiler. */
322	n/a	double q = y - absx;
323	n/a	z = q - lanczos_g_minus_half;
324	n/a	}
325	n/a	else {
326	n/a	double q = y - lanczos_g_minus_half;
327	n/a	z = q - absx;
328	n/a	}
329	n/a	z = z * lanczos_g / y;
330	n/a	if (x < 0.0) {
331	n/a	r = -pi / sinpi(absx) / absx * exp(y) / lanczos_sum(absx);
332	n/a	r -= z * r;
333	n/a	if (absx < 140.0) {
334	n/a	r /= pow(y, absx - 0.5);
335	n/a	}
336	n/a	else {
337	n/a	sqrtpow = pow(y, absx / 2.0 - 0.25);
338	n/a	r /= sqrtpow;
339	n/a	r /= sqrtpow;
340	n/a	}
341	n/a	}
342	n/a	else {
343	n/a	r = lanczos_sum(absx) / exp(y);
344	n/a	r += z * r;
345	n/a	if (absx < 140.0) {
346	n/a	r *= pow(y, absx - 0.5);
347	n/a	}
348	n/a	else {
349	n/a	sqrtpow = pow(y, absx / 2.0 - 0.25);
350	n/a	r *= sqrtpow;
351	n/a	r *= sqrtpow;
352	n/a	}
353	n/a	}
354	n/a	if (Py_IS_INFINITY(r))
355	n/a	errno = ERANGE;
356	n/a	return r;
357	n/a	}
358	n/a
359	n/a	/*
360	n/a	lgamma: natural log of the absolute value of the Gamma function.
361	n/a	For large arguments, Lanczos' formula works extremely well here.
362	n/a	*/
363	n/a
364	n/a	static double
365	n/a	m_lgamma(double x)
366	n/a	{
367	n/a	double r, absx;
368	n/a
369	n/a	/* special cases */
370	n/a	if (!Py_IS_FINITE(x)) {
371	n/a	if (Py_IS_NAN(x))
372	n/a	return x; /* lgamma(nan) = nan */
373	n/a	else
374	n/a	return Py_HUGE_VAL; /* lgamma(+-inf) = +inf */
375	n/a	}
376	n/a
377	n/a	/* integer arguments */
378	n/a	if (x == floor(x) && x <= 2.0) {
379	n/a	if (x <= 0.0) {
380	n/a	errno = EDOM; /* lgamma(n) = inf, divide-by-zero for */
381	n/a	return Py_HUGE_VAL; /* integers n <= 0 */
382	n/a	}
383	n/a	else {
384	n/a	return 0.0; /* lgamma(1) = lgamma(2) = 0.0 */
385	n/a	}
386	n/a	}
387	n/a
388	n/a	absx = fabs(x);
389	n/a	/* tiny arguments: lgamma(x) ~ -log(fabs(x)) for small x */
390	n/a	if (absx < 1e-20)
391	n/a	return -log(absx);
392	n/a
393	n/a	/* Lanczos' formula. We could save a fraction of a ulp in accuracy by
394	n/a	having a second set of numerator coefficients for lanczos_sum that
395	n/a	absorbed the exp(-lanczos_g) term, and throwing out the lanczos_g
396	n/a	subtraction below; it's probably not worth it. */
397	n/a	r = log(lanczos_sum(absx)) - lanczos_g;
398	n/a	r += (absx - 0.5) * (log(absx + lanczos_g - 0.5) - 1);
399	n/a	if (x < 0.0)
400	n/a	/* Use reflection formula to get value for negative x. */
401	n/a	r = logpi - log(fabs(sinpi(absx))) - log(absx) - r;
402	n/a	if (Py_IS_INFINITY(r))
403	n/a	errno = ERANGE;
404	n/a	return r;
405	n/a	}
406	n/a
407	n/a	/*
408	n/a	Implementations of the error function erf(x) and the complementary error
409	n/a	function erfc(x).
410	n/a
411	n/a	Method: we use a series approximation for erf for small x, and a continued
412	n/a	fraction approximation for erfc(x) for larger x;
413	n/a	combined with the relations erf(-x) = -erf(x) and erfc(x) = 1.0 - erf(x),
414	n/a	this gives us erf(x) and erfc(x) for all x.
415	n/a
416	n/a	The series expansion used is:
417	n/a
418	n/a	erf(x) = xexp(-xx)/sqrt(pi) * [
419	n/a	2/1 + 4/3 x2 + 8/15 x4 + 16/105 x**6 + ...]
420	n/a
421	n/a	The coefficient of x(2k-2) here is 4kfactorial(k)/factorial(2k).
422	n/a	This series converges well for smallish x, but slowly for larger x.
423	n/a
424	n/a	The continued fraction expansion used is:
425	n/a
426	n/a	erfc(x) = xexp(-xx)/sqrt(pi) * [1/(0.5 + x2 -) 0.5/(2.5 + x2 - )
427	n/a	3.0/(4.5 + x2 - ) 7.5/(6.5 + x2 - ) ...]
428	n/a
429	n/a	after the first term, the general term has the form:
430	n/a
431	n/a	k(k-0.5)/(2k+0.5 + x**2 - ...).
432	n/a
433	n/a	This expansion converges fast for larger x, but convergence becomes
434	n/a	infinitely slow as x approaches 0.0. The (somewhat naive) continued
435	n/a	fraction evaluation algorithm used below also risks overflow for large x;
436	n/a	but for large x, erfc(x) == 0.0 to within machine precision. (For
437	n/a	example, erfc(30.0) is approximately 2.56e-393).
438	n/a
439	n/a	Parameters: use series expansion for abs(x) < ERF_SERIES_CUTOFF and
440	n/a	continued fraction expansion for ERF_SERIES_CUTOFF <= abs(x) <
441	n/a	ERFC_CONTFRAC_CUTOFF. ERFC_SERIES_TERMS and ERFC_CONTFRAC_TERMS are the
442	n/a	numbers of terms to use for the relevant expansions. */
443	n/a
444	n/a	#define ERF_SERIES_CUTOFF 1.5
445	n/a	#define ERF_SERIES_TERMS 25
446	n/a	#define ERFC_CONTFRAC_CUTOFF 30.0
447	n/a	#define ERFC_CONTFRAC_TERMS 50
448	n/a
449	n/a	/*
450	n/a	Error function, via power series.
451	n/a
452	n/a	Given a finite float x, return an approximation to erf(x).
453	n/a	Converges reasonably fast for small x.
454	n/a	*/
455	n/a
456	n/a	static double
457	n/a	m_erf_series(double x)
458	n/a	{
459	n/a	double x2, acc, fk, result;
460	n/a	int i, saved_errno;
461	n/a
462	n/a	x2 = x * x;
463	n/a	acc = 0.0;
464	n/a	fk = (double)ERF_SERIES_TERMS + 0.5;
465	n/a	for (i = 0; i < ERF_SERIES_TERMS; i++) {
466	n/a	acc = 2.0 + x2 * acc / fk;
467	n/a	fk -= 1.0;
468	n/a	}
469	n/a	/* Make sure the exp call doesn't affect errno;
470	n/a	see m_erfc_contfrac for more. */
471	n/a	saved_errno = errno;
472	n/a	result = acc * x * exp(-x2) / sqrtpi;
473	n/a	errno = saved_errno;
474	n/a	return result;
475	n/a	}
476	n/a
477	n/a	/*
478	n/a	Complementary error function, via continued fraction expansion.
479	n/a
480	n/a	Given a positive float x, return an approximation to erfc(x). Converges
481	n/a	reasonably fast for x large (say, x > 2.0), and should be safe from
482	n/a	overflow if x and nterms are not too large. On an IEEE 754 machine, with x
483	n/a	<= 30.0, we're safe up to nterms = 100. For x >= 30.0, erfc(x) is smaller
484	n/a	than the smallest representable nonzero float. */
485	n/a
486	n/a	static double
487	n/a	m_erfc_contfrac(double x)
488	n/a	{
489	n/a	double x2, a, da, p, p_last, q, q_last, b, result;
490	n/a	int i, saved_errno;
491	n/a
492	n/a	if (x >= ERFC_CONTFRAC_CUTOFF)
493	n/a	return 0.0;
494	n/a
495	n/a	x2 = x*x;
496	n/a	a = 0.0;
497	n/a	da = 0.5;
498	n/a	p = 1.0; p_last = 0.0;
499	n/a	q = da + x2; q_last = 1.0;
500	n/a	for (i = 0; i < ERFC_CONTFRAC_TERMS; i++) {
501	n/a	double temp;
502	n/a	a += da;
503	n/a	da += 2.0;
504	n/a	b = da + x2;
505	n/a	temp = p; p = bp - ap_last; p_last = temp;
506	n/a	temp = q; q = bq - aq_last; q_last = temp;
507	n/a	}
508	n/a	/* Issue #8986: On some platforms, exp sets errno on underflow to zero;
509	n/a	save the current errno value so that we can restore it later. */
510	n/a	saved_errno = errno;
511	n/a	result = p / q * x * exp(-x2) / sqrtpi;
512	n/a	errno = saved_errno;
513	n/a	return result;
514	n/a	}
515	n/a
516	n/a	/* Error function erf(x), for general x */
517	n/a
518	n/a	static double
519	n/a	m_erf(double x)
520	n/a	{
521	n/a	double absx, cf;
522	n/a
523	n/a	if (Py_IS_NAN(x))
524	n/a	return x;
525	n/a	absx = fabs(x);
526	n/a	if (absx < ERF_SERIES_CUTOFF)
527	n/a	return m_erf_series(x);
528	n/a	else {
529	n/a	cf = m_erfc_contfrac(absx);
530	n/a	return x > 0.0 ? 1.0 - cf : cf - 1.0;
531	n/a	}
532	n/a	}
533	n/a
534	n/a	/* Complementary error function erfc(x), for general x. */
535	n/a
536	n/a	static double
537	n/a	m_erfc(double x)
538	n/a	{
539	n/a	double absx, cf;
540	n/a
541	n/a	if (Py_IS_NAN(x))
542	n/a	return x;
543	n/a	absx = fabs(x);
544	n/a	if (absx < ERF_SERIES_CUTOFF)
545	n/a	return 1.0 - m_erf_series(x);
546	n/a	else {
547	n/a	cf = m_erfc_contfrac(absx);
548	n/a	return x > 0.0 ? cf : 2.0 - cf;
549	n/a	}
550	n/a	}
551	n/a
552	n/a	/*
553	n/a	wrapper for atan2 that deals directly with special cases before
554	n/a	delegating to the platform libm for the remaining cases. This
555	n/a	is necessary to get consistent behaviour across platforms.
556	n/a	Windows, FreeBSD and alpha Tru64 are amongst platforms that don't
557	n/a	always follow C99.
558	n/a	*/
559	n/a
560	n/a	static double
561	n/a	m_atan2(double y, double x)
562	n/a	{
563	n/a	if (Py_IS_NAN(x) \|\| Py_IS_NAN(y))
564	n/a	return Py_NAN;
565	n/a	if (Py_IS_INFINITY(y)) {
566	n/a	if (Py_IS_INFINITY(x)) {
567	n/a	if (copysign(1., x) == 1.)
568	n/a	/* atan2(+-inf, +inf) == +-pi/4 */
569	n/a	return copysign(0.25*Py_MATH_PI, y);
570	n/a	else
571	n/a	/* atan2(+-inf, -inf) == +-pi3/4 /
572	n/a	return copysign(0.75*Py_MATH_PI, y);
573	n/a	}
574	n/a	/* atan2(+-inf, x) == +-pi/2 for finite x */
575	n/a	return copysign(0.5*Py_MATH_PI, y);
576	n/a	}
577	n/a	if (Py_IS_INFINITY(x) \|\| y == 0.) {
578	n/a	if (copysign(1., x) == 1.)
579	n/a	/* atan2(+-y, +inf) = atan2(+-0, +x) = +-0. */
580	n/a	return copysign(0., y);
581	n/a	else
582	n/a	/* atan2(+-y, -inf) = atan2(+-0., -x) = +-pi. */
583	n/a	return copysign(Py_MATH_PI, y);
584	n/a	}
585	n/a	return atan2(y, x);
586	n/a	}
587	n/a
588	n/a	/*
589	n/a	Various platforms (Solaris, OpenBSD) do nonstandard things for log(0),
590	n/a	log(-ve), log(NaN). Here are wrappers for log and log10 that deal with
591	n/a	special values directly, passing positive non-special values through to
592	n/a	the system log/log10.
593	n/a	*/
594	n/a
595	n/a	static double
596	n/a	m_log(double x)
597	n/a	{
598	n/a	if (Py_IS_FINITE(x)) {
599	n/a	if (x > 0.0)
600	n/a	return log(x);
601	n/a	errno = EDOM;
602	n/a	if (x == 0.0)
603	n/a	return -Py_HUGE_VAL; /* log(0) = -inf */
604	n/a	else
605	n/a	return Py_NAN; /* log(-ve) = nan */
606	n/a	}
607	n/a	else if (Py_IS_NAN(x))
608	n/a	return x; /* log(nan) = nan */
609	n/a	else if (x > 0.0)
610	n/a	return x; /* log(inf) = inf */
611	n/a	else {
612	n/a	errno = EDOM;
613	n/a	return Py_NAN; /* log(-inf) = nan */
614	n/a	}
615	n/a	}
616	n/a
617	n/a	/*
618	n/a	log2: log to base 2.
619	n/a
620	n/a	Uses an algorithm that should:
621	n/a
622	n/a	(a) produce exact results for powers of 2, and
623	n/a	(b) give a monotonic log2 (for positive finite floats),
624	n/a	assuming that the system log is monotonic.
625	n/a	*/
626	n/a
627	n/a	static double
628	n/a	m_log2(double x)
629	n/a	{
630	n/a	if (!Py_IS_FINITE(x)) {
631	n/a	if (Py_IS_NAN(x))
632	n/a	return x; /* log2(nan) = nan */
633	n/a	else if (x > 0.0)
634	n/a	return x; /* log2(+inf) = +inf */
635	n/a	else {
636	n/a	errno = EDOM;
637	n/a	return Py_NAN; /* log2(-inf) = nan, invalid-operation */
638	n/a	}
639	n/a	}
640	n/a
641	n/a	if (x > 0.0) {
642	n/a	#ifdef HAVE_LOG2
643	n/a	return log2(x);
644	n/a	#else
645	n/a	double m;
646	n/a	int e;
647	n/a	m = frexp(x, &e);
648	n/a	/* We want log2(m * 2**e) == log(m) / log(2) + e. Care is needed when
649	n/a	* x is just greater than 1.0: in that case e is 1, log(m) is negative,
650	n/a	* and we get significant cancellation error from the addition of
651	n/a	* log(m) / log(2) to e. The slight rewrite of the expression below
652	n/a	* avoids this problem.
653	n/a	*/
654	n/a	if (x >= 1.0) {
655	n/a	return log(2.0 * m) / log(2.0) + (e - 1);
656	n/a	}
657	n/a	else {
658	n/a	return log(m) / log(2.0) + e;
659	n/a	}
660	n/a	#endif
661	n/a	}
662	n/a	else if (x == 0.0) {
663	n/a	errno = EDOM;
664	n/a	return -Py_HUGE_VAL; /* log2(0) = -inf, divide-by-zero */
665	n/a	}
666	n/a	else {
667	n/a	errno = EDOM;
668	n/a	return Py_NAN; /* log2(-inf) = nan, invalid-operation */
669	n/a	}
670	n/a	}
671	n/a
672	n/a	static double
673	n/a	m_log10(double x)
674	n/a	{
675	n/a	if (Py_IS_FINITE(x)) {
676	n/a	if (x > 0.0)
677	n/a	return log10(x);
678	n/a	errno = EDOM;
679	n/a	if (x == 0.0)
680	n/a	return -Py_HUGE_VAL; /* log10(0) = -inf */
681	n/a	else
682	n/a	return Py_NAN; /* log10(-ve) = nan */
683	n/a	}
684	n/a	else if (Py_IS_NAN(x))
685	n/a	return x; /* log10(nan) = nan */
686	n/a	else if (x > 0.0)
687	n/a	return x; /* log10(inf) = inf */
688	n/a	else {
689	n/a	errno = EDOM;
690	n/a	return Py_NAN; /* log10(-inf) = nan */
691	n/a	}
692	n/a	}
693	n/a
694	n/a
695	n/a	/*[clinic input]
696	n/a	math.gcd
697	n/a
698	n/a	x as a: object
699	n/a	y as b: object
700	n/a	/
701	n/a
702	n/a	greatest common divisor of x and y
703	n/a	[clinic start generated code]*/
704	n/a
705	n/a	static PyObject *
706	n/a	math_gcd_impl(PyObject module, PyObject a, PyObject *b)
707	n/a	/[clinic end generated code: output=7b2e0c151bd7a5d8 input=c2691e57fb2a98fa]/
708	n/a	{
709	n/a	PyObject *g;
710	n/a
711	n/a	a = PyNumber_Index(a);
712	n/a	if (a == NULL)
713	n/a	return NULL;
714	n/a	b = PyNumber_Index(b);
715	n/a	if (b == NULL) {
716	n/a	Py_DECREF(a);
717	n/a	return NULL;
718	n/a	}
719	n/a	g = _PyLong_GCD(a, b);
720	n/a	Py_DECREF(a);
721	n/a	Py_DECREF(b);
722	n/a	return g;
723	n/a	}
724	n/a
725	n/a
726	n/a	/* Call is_error when errno != 0, and where x is the result libm
727	n/a	* returned. is_error will usually set up an exception and return
728	n/a	* true (1), but may return false (0) without setting up an exception.
729	n/a	*/
730	n/a	static int
731	n/a	is_error(double x)
732	n/a	{
733	n/a	int result = 1; /* presumption of guilt */
734	n/a	assert(errno); /* non-zero errno is a precondition for calling */
735	n/a	if (errno == EDOM)
736	n/a	PyErr_SetString(PyExc_ValueError, "math domain error");
737	n/a
738	n/a	else if (errno == ERANGE) {
739	n/a	/* ANSI C generally requires libm functions to set ERANGE
740	n/a	* on overflow, but also generally allows them to set
741	n/a	* ERANGE on underflow too. There's no consistency about
742	n/a	* the latter across platforms.
743	n/a	* Alas, C99 never requires that errno be set.
744	n/a	* Here we suppress the underflow errors (libm functions
745	n/a	* should return a zero on underflow, and +- HUGE_VAL on
746	n/a	* overflow, so testing the result for zero suffices to
747	n/a	* distinguish the cases).
748	n/a	*
749	n/a	* On some platforms (Ubuntu/ia64) it seems that errno can be
750	n/a	* set to ERANGE for subnormal results that do not underflow
751	n/a	* to zero. So to be safe, we'll ignore ERANGE whenever the
752	n/a	* function result is less than one in absolute value.
753	n/a	*/
754	n/a	if (fabs(x) < 1.0)
755	n/a	result = 0;
756	n/a	else
757	n/a	PyErr_SetString(PyExc_OverflowError,
758	n/a	"math range error");
759	n/a	}
760	n/a	else
761	n/a	/* Unexpected math error */
762	n/a	PyErr_SetFromErrno(PyExc_ValueError);
763	n/a	return result;
764	n/a	}
765	n/a
766	n/a	/*
767	n/a	math_1 is used to wrap a libm function f that takes a double
768	n/a	argument and returns a double.
769	n/a
770	n/a	The error reporting follows these rules, which are designed to do
771	n/a	the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
772	n/a	platforms.
773	n/a
774	n/a	- a NaN result from non-NaN inputs causes ValueError to be raised
775	n/a	- an infinite result from finite inputs causes OverflowError to be
776	n/a	raised if can_overflow is 1, or raises ValueError if can_overflow
777	n/a	is 0.
778	n/a	- if the result is finite and errno == EDOM then ValueError is
779	n/a	raised
780	n/a	- if the result is finite and nonzero and errno == ERANGE then
781	n/a	OverflowError is raised
782	n/a
783	n/a	The last rule is used to catch overflow on platforms which follow
784	n/a	C89 but for which HUGE_VAL is not an infinity.
785	n/a
786	n/a	For the majority of one-argument functions these rules are enough
787	n/a	to ensure that Python's functions behave as specified in 'Annex F'
788	n/a	of the C99 standard, with the 'invalid' and 'divide-by-zero'
789	n/a	floating-point exceptions mapping to Python's ValueError and the
790	n/a	'overflow' floating-point exception mapping to OverflowError.
791	n/a	math_1 only works for functions that don't have singularities and
792	n/a	the possibility of overflow; fortunately, that covers everything we
793	n/a	care about right now.
794	n/a	*/
795	n/a
796	n/a	static PyObject *
797	n/a	math_1_to_whatever(PyObject arg, double (func) (double),
798	n/a	PyObject (from_double_func) (double),
799	n/a	int can_overflow)
800	n/a	{
801	n/a	double x, r;
802	n/a	x = PyFloat_AsDouble(arg);
803	n/a	if (x == -1.0 && PyErr_Occurred())
804	n/a	return NULL;
805	n/a	errno = 0;
806	n/a	PyFPE_START_PROTECT("in math_1", return 0);
807	n/a	r = (*func)(x);
808	n/a	PyFPE_END_PROTECT(r);
809	n/a	if (Py_IS_NAN(r) && !Py_IS_NAN(x)) {
810	n/a	PyErr_SetString(PyExc_ValueError,
811	n/a	"math domain error"); /* invalid arg */
812	n/a	return NULL;
813	n/a	}
814	n/a	if (Py_IS_INFINITY(r) && Py_IS_FINITE(x)) {
815	n/a	if (can_overflow)
816	n/a	PyErr_SetString(PyExc_OverflowError,
817	n/a	"math range error"); /* overflow */
818	n/a	else
819	n/a	PyErr_SetString(PyExc_ValueError,
820	n/a	"math domain error"); /* singularity */
821	n/a	return NULL;
822	n/a	}
823	n/a	if (Py_IS_FINITE(r) && errno && is_error(r))
824	n/a	/* this branch unnecessary on most platforms */
825	n/a	return NULL;
826	n/a
827	n/a	return (*from_double_func)(r);
828	n/a	}
829	n/a
830	n/a	/* variant of math_1, to be used when the function being wrapped is known to
831	n/a	set errno properly (that is, errno = EDOM for invalid or divide-by-zero,
832	n/a	errno = ERANGE for overflow). */
833	n/a
834	n/a	static PyObject *
835	n/a	math_1a(PyObject arg, double (func) (double))
836	n/a	{
837	n/a	double x, r;
838	n/a	x = PyFloat_AsDouble(arg);
839	n/a	if (x == -1.0 && PyErr_Occurred())
840	n/a	return NULL;
841	n/a	errno = 0;
842	n/a	PyFPE_START_PROTECT("in math_1a", return 0);
843	n/a	r = (*func)(x);
844	n/a	PyFPE_END_PROTECT(r);
845	n/a	if (errno && is_error(r))
846	n/a	return NULL;
847	n/a	return PyFloat_FromDouble(r);
848	n/a	}
849	n/a
850	n/a	/*
851	n/a	math_2 is used to wrap a libm function f that takes two double
852	n/a	arguments and returns a double.
853	n/a
854	n/a	The error reporting follows these rules, which are designed to do
855	n/a	the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
856	n/a	platforms.
857	n/a
858	n/a	- a NaN result from non-NaN inputs causes ValueError to be raised
859	n/a	- an infinite result from finite inputs causes OverflowError to be
860	n/a	raised.
861	n/a	- if the result is finite and errno == EDOM then ValueError is
862	n/a	raised
863	n/a	- if the result is finite and nonzero and errno == ERANGE then
864	n/a	OverflowError is raised
865	n/a
866	n/a	The last rule is used to catch overflow on platforms which follow
867	n/a	C89 but for which HUGE_VAL is not an infinity.
868	n/a
869	n/a	For most two-argument functions (copysign, fmod, hypot, atan2)
870	n/a	these rules are enough to ensure that Python's functions behave as
871	n/a	specified in 'Annex F' of the C99 standard, with the 'invalid' and
872	n/a	'divide-by-zero' floating-point exceptions mapping to Python's
873	n/a	ValueError and the 'overflow' floating-point exception mapping to
874	n/a	OverflowError.
875	n/a	*/
876	n/a
877	n/a	static PyObject *
878	n/a	math_1(PyObject arg, double (func) (double), int can_overflow)
879	n/a	{
880	n/a	return math_1_to_whatever(arg, func, PyFloat_FromDouble, can_overflow);
881	n/a	}
882	n/a
883	n/a	static PyObject *
884	n/a	math_1_to_int(PyObject arg, double (func) (double), int can_overflow)
885	n/a	{
886	n/a	return math_1_to_whatever(arg, func, PyLong_FromDouble, can_overflow);
887	n/a	}
888	n/a
889	n/a	static PyObject *
890	n/a	math_2(PyObject args, double (func) (double, double), const char *funcname)
891	n/a	{
892	n/a	PyObject ox, oy;
893	n/a	double x, y, r;
894	n/a	if (! PyArg_UnpackTuple(args, funcname, 2, 2, &ox, &oy))
895	n/a	return NULL;
896	n/a	x = PyFloat_AsDouble(ox);
897	n/a	y = PyFloat_AsDouble(oy);
898	n/a	if ((x == -1.0 \|\| y == -1.0) && PyErr_Occurred())
899	n/a	return NULL;
900	n/a	errno = 0;
901	n/a	PyFPE_START_PROTECT("in math_2", return 0);
902	n/a	r = (*func)(x, y);
903	n/a	PyFPE_END_PROTECT(r);
904	n/a	if (Py_IS_NAN(r)) {
905	n/a	if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
906	n/a	errno = EDOM;
907	n/a	else
908	n/a	errno = 0;
909	n/a	}
910	n/a	else if (Py_IS_INFINITY(r)) {
911	n/a	if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
912	n/a	errno = ERANGE;
913	n/a	else
914	n/a	errno = 0;
915	n/a	}
916	n/a	if (errno && is_error(r))
917	n/a	return NULL;
918	n/a	else
919	n/a	return PyFloat_FromDouble(r);
920	n/a	}
921	n/a
922	n/a	#define FUNC1(funcname, func, can_overflow, docstring) \
923	n/a	static PyObject * math_##funcname(PyObject self, PyObject args) { \
924	n/a	return math_1(args, func, can_overflow); \
925	n/a	}\
926	n/a	PyDoc_STRVAR(math_##funcname##_doc, docstring);
927	n/a
928	n/a	#define FUNC1A(funcname, func, docstring) \
929	n/a	static PyObject * math_##funcname(PyObject self, PyObject args) { \
930	n/a	return math_1a(args, func); \
931	n/a	}\
932	n/a	PyDoc_STRVAR(math_##funcname##_doc, docstring);
933	n/a
934	n/a	#define FUNC2(funcname, func, docstring) \
935	n/a	static PyObject * math_##funcname(PyObject self, PyObject args) { \
936	n/a	return math_2(args, func, #funcname); \
937	n/a	}\
938	n/a	PyDoc_STRVAR(math_##funcname##_doc, docstring);
939	n/a
940	n/a	FUNC1(acos, acos, 0,
941	n/a	"acos($module, x, /)\n--\n\n"
942	n/a	"Return the arc cosine (measured in radians) of x.")
943	n/a	FUNC1(acosh, m_acosh, 0,
944	n/a	"acosh($module, x, /)\n--\n\n"
945	n/a	"Return the inverse hyperbolic cosine of x.")
946	n/a	FUNC1(asin, asin, 0,
947	n/a	"asin($module, x, /)\n--\n\n"
948	n/a	"Return the arc sine (measured in radians) of x.")
949	n/a	FUNC1(asinh, m_asinh, 0,
950	n/a	"asinh($module, x, /)\n--\n\n"
951	n/a	"Return the inverse hyperbolic sine of x.")
952	n/a	FUNC1(atan, atan, 0,
953	n/a	"atan($module, x, /)\n--\n\n"
954	n/a	"Return the arc tangent (measured in radians) of x.")
955	n/a	FUNC2(atan2, m_atan2,
956	n/a	"atan2($module, y, x, /)\n--\n\n"
957	n/a	"Return the arc tangent (measured in radians) of y/x.\n\n"
958	n/a	"Unlike atan(y/x), the signs of both x and y are considered.")
959	n/a	FUNC1(atanh, m_atanh, 0,
960	n/a	"atanh($module, x, /)\n--\n\n"
961	n/a	"Return the inverse hyperbolic tangent of x.")
962	n/a
963	n/a	/*[clinic input]
964	n/a	math.ceil
965	n/a
966	n/a	x as number: object
967	n/a	/
968	n/a
969	n/a	Return the ceiling of x as an Integral.
970	n/a
971	n/a	This is the smallest integer >= x.
972	n/a	[clinic start generated code]*/
973	n/a
974	n/a	static PyObject *
975	n/a	math_ceil(PyObject module, PyObject number)
976	n/a	/[clinic end generated code: output=6c3b8a78bc201c67 input=2725352806399cab]/
977	n/a	{
978	n/a	_Py_IDENTIFIER(__ceil__);
979	n/a	PyObject method, result;
980	n/a
981	n/a	method = _PyObject_LookupSpecial(number, &PyId___ceil__);
982	n/a	if (method == NULL) {
983	n/a	if (PyErr_Occurred())
984	n/a	return NULL;
985	n/a	return math_1_to_int(number, ceil, 0);
986	n/a	}
987	n/a	result = _PyObject_CallNoArg(method);
988	n/a	Py_DECREF(method);
989	n/a	return result;
990	n/a	}
991	n/a
992	n/a	FUNC2(copysign, copysign,
993	n/a	"copysign($module, x, y, /)\n--\n\n"
994	n/a	"Return a float with the magnitude (absolute value) of x but the sign of y.\n\n"
995	n/a	"On platforms that support signed zeros, copysign(1.0, -0.0)\n"
996	n/a	"returns -1.0.\n")
997	n/a	FUNC1(cos, cos, 0,
998	n/a	"cos($module, x, /)\n--\n\n"
999	n/a	"Return the cosine of x (measured in radians).")
1000	n/a	FUNC1(cosh, cosh, 1,
1001	n/a	"cosh($module, x, /)\n--\n\n"
1002	n/a	"Return the hyperbolic cosine of x.")
1003	n/a	FUNC1A(erf, m_erf,
1004	n/a	"erf($module, x, /)\n--\n\n"
1005	n/a	"Error function at x.")
1006	n/a	FUNC1A(erfc, m_erfc,
1007	n/a	"erfc($module, x, /)\n--\n\n"
1008	n/a	"Complementary error function at x.")
1009	n/a	FUNC1(exp, exp, 1,
1010	n/a	"exp($module, x, /)\n--\n\n"
1011	n/a	"Return e raised to the power of x.")
1012	n/a	FUNC1(expm1, m_expm1, 1,
1013	n/a	"expm1($module, x, /)\n--\n\n"
1014	n/a	"Return exp(x)-1.\n\n"
1015	n/a	"This function avoids the loss of precision involved in the direct "
1016	n/a	"evaluation of exp(x)-1 for small x.")
1017	n/a	FUNC1(fabs, fabs, 0,
1018	n/a	"fabs($module, x, /)\n--\n\n"
1019	n/a	"Return the absolute value of the float x.")
1020	n/a
1021	n/a	/*[clinic input]
1022	n/a	math.floor
1023	n/a
1024	n/a	x as number: object
1025	n/a	/
1026	n/a
1027	n/a	Return the floor of x as an Integral.
1028	n/a
1029	n/a	This is the largest integer <= x.
1030	n/a	[clinic start generated code]*/
1031	n/a
1032	n/a	static PyObject *
1033	n/a	math_floor(PyObject module, PyObject number)
1034	n/a	/[clinic end generated code: output=c6a65c4884884b8a input=63af6b5d7ebcc3d6]/
1035	n/a	{
1036	n/a	_Py_IDENTIFIER(__floor__);
1037	n/a	PyObject method, result;
1038	n/a
1039	n/a	method = _PyObject_LookupSpecial(number, &PyId___floor__);
1040	n/a	if (method == NULL) {
1041	n/a	if (PyErr_Occurred())
1042	n/a	return NULL;
1043	n/a	return math_1_to_int(number, floor, 0);
1044	n/a	}
1045	n/a	result = _PyObject_CallNoArg(method);
1046	n/a	Py_DECREF(method);
1047	n/a	return result;
1048	n/a	}
1049	n/a
1050	n/a	FUNC1A(gamma, m_tgamma,
1051	n/a	"gamma($module, x, /)\n--\n\n"
1052	n/a	"Gamma function at x.")
1053	n/a	FUNC1A(lgamma, m_lgamma,
1054	n/a	"lgamma($module, x, /)\n--\n\n"
1055	n/a	"Natural logarithm of absolute value of Gamma function at x.")
1056	n/a	FUNC1(log1p, m_log1p, 0,
1057	n/a	"log1p($module, x, /)\n--\n\n"
1058	n/a	"Return the natural logarithm of 1+x (base e).\n\n"
1059	n/a	"The result is computed in a way which is accurate for x near zero.")
1060	n/a	FUNC1(sin, sin, 0,
1061	n/a	"sin($module, x, /)\n--\n\n"
1062	n/a	"Return the sine of x (measured in radians).")
1063	n/a	FUNC1(sinh, sinh, 1,
1064	n/a	"sinh($module, x, /)\n--\n\n"
1065	n/a	"Return the hyperbolic sine of x.")
1066	n/a	FUNC1(sqrt, sqrt, 0,
1067	n/a	"sqrt($module, x, /)\n--\n\n"
1068	n/a	"Return the square root of x.")
1069	n/a	FUNC1(tan, tan, 0,
1070	n/a	"tan($module, x, /)\n--\n\n"
1071	n/a	"Return the tangent of x (measured in radians).")
1072	n/a	FUNC1(tanh, tanh, 0,
1073	n/a	"tanh($module, x, /)\n--\n\n"
1074	n/a	"Return the hyperbolic tangent of x.")
1075	n/a
1076	n/a	/* Precision summation function as msum() by Raymond Hettinger in
1077	n/a	<http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/393090>,
1078	n/a	enhanced with the exact partials sum and roundoff from Mark
1079	n/a	Dickinson's post at <http://bugs.python.org/file10357/msum4.py>.
1080	n/a	See those links for more details, proofs and other references.
1081	n/a
1082	n/a	Note 1: IEEE 754R floating point semantics are assumed,
1083	n/a	but the current implementation does not re-establish special
1084	n/a	value semantics across iterations (i.e. handling -Inf + Inf).
1085	n/a
1086	n/a	Note 2: No provision is made for intermediate overflow handling;
1087	n/a	therefore, sum([1e+308, 1e-308, 1e+308]) returns 1e+308 while
1088	n/a	sum([1e+308, 1e+308, 1e-308]) raises an OverflowError due to the
1089	n/a	overflow of the first partial sum.
1090	n/a
1091	n/a	Note 3: The intermediate values lo, yr, and hi are declared volatile so
1092	n/a	aggressive compilers won't algebraically reduce lo to always be exactly 0.0.
1093	n/a	Also, the volatile declaration forces the values to be stored in memory as
1094	n/a	regular doubles instead of extended long precision (80-bit) values. This
1095	n/a	prevents double rounding because any addition or subtraction of two doubles
1096	n/a	can be resolved exactly into double-sized hi and lo values. As long as the
1097	n/a	hi value gets forced into a double before yr and lo are computed, the extra
1098	n/a	bits in downstream extended precision operations (x87 for example) will be
1099	n/a	exactly zero and therefore can be losslessly stored back into a double,
1100	n/a	thereby preventing double rounding.
1101	n/a
1102	n/a	Note 4: A similar implementation is in Modules/cmathmodule.c.
1103	n/a	Be sure to update both when making changes.
1104	n/a
1105	n/a	Note 5: The signature of math.fsum() differs from builtins.sum()
1106	n/a	because the start argument doesn't make sense in the context of
1107	n/a	accurate summation. Since the partials table is collapsed before
1108	n/a	returning a result, sum(seq2, start=sum(seq1)) may not equal the
1109	n/a	accurate result returned by sum(itertools.chain(seq1, seq2)).
1110	n/a	*/
1111	n/a
1112	n/a	#define NUM_PARTIALS 32 /* initial partials array size, on stack */
1113	n/a
1114	n/a	/* Extend the partials array p[] by doubling its size. */
1115	n/a	static int /* non-zero on error */
1116	n/a	_fsum_realloc(double **p_ptr, Py_ssize_t n,
1117	n/a	double ps, Py_ssize_t m_ptr)
1118	n/a	{
1119	n/a	void *v = NULL;
1120	n/a	Py_ssize_t m = *m_ptr;
1121	n/a
1122	n/a	m += m; /* double */
1123	n/a	if (n < m && (size_t)m < ((size_t)PY_SSIZE_T_MAX / sizeof(double))) {
1124	n/a	double p = p_ptr;
1125	n/a	if (p == ps) {
1126	n/a	v = PyMem_Malloc(sizeof(double) * m);
1127	n/a	if (v != NULL)
1128	n/a	memcpy(v, ps, sizeof(double) * n);
1129	n/a	}
1130	n/a	else
1131	n/a	v = PyMem_Realloc(p, sizeof(double) * m);
1132	n/a	}
1133	n/a	if (v == NULL) { /* size overflow or no memory */
1134	n/a	PyErr_SetString(PyExc_MemoryError, "math.fsum partials");
1135	n/a	return 1;
1136	n/a	}
1137	n/a	p_ptr = (double) v;
1138	n/a	*m_ptr = m;
1139	n/a	return 0;
1140	n/a	}
1141	n/a
1142	n/a	/* Full precision summation of a sequence of floats.
1143	n/a
1144	n/a	def msum(iterable):
1145	n/a	partials = [] # sorted, non-overlapping partial sums
1146	n/a	for x in iterable:
1147	n/a	i = 0
1148	n/a	for y in partials:
1149	n/a	if abs(x) < abs(y):
1150	n/a	x, y = y, x
1151	n/a	hi = x + y
1152	n/a	lo = y - (hi - x)
1153	n/a	if lo:
1154	n/a	partials[i] = lo
1155	n/a	i += 1
1156	n/a	x = hi
1157	n/a	partials[i:] = [x]
1158	n/a	return sum_exact(partials)
1159	n/a
1160	n/a	Rounded x+y stored in hi with the roundoff stored in lo. Together hi+lo
1161	n/a	are exactly equal to x+y. The inner loop applies hi/lo summation to each
1162	n/a	partial so that the list of partial sums remains exact.
1163	n/a
1164	n/a	Sum_exact() adds the partial sums exactly and correctly rounds the final
1165	n/a	result (using the round-half-to-even rule). The items in partials remain
1166	n/a	non-zero, non-special, non-overlapping and strictly increasing in
1167	n/a	magnitude, but possibly not all having the same sign.
1168	n/a
1169	n/a	Depends on IEEE 754 arithmetic guarantees and half-even rounding.
1170	n/a	*/
1171	n/a
1172	n/a	/*[clinic input]
1173	n/a	math.fsum
1174	n/a
1175	n/a	seq: object
1176	n/a	/
1177	n/a
1178	n/a	Return an accurate floating point sum of values in the iterable seq.
1179	n/a
1180	n/a	Assumes IEEE-754 floating point arithmetic.
1181	n/a	[clinic start generated code]*/
1182	n/a
1183	n/a	static PyObject *
1184	n/a	math_fsum(PyObject module, PyObject seq)
1185	n/a	/[clinic end generated code: output=ba5c672b87fe34fc input=c51b7d8caf6f6e82]/
1186	n/a	{
1187	n/a	PyObject item, iter, *sum = NULL;
1188	n/a	Py_ssize_t i, j, n = 0, m = NUM_PARTIALS;
1189	n/a	double x, y, t, ps[NUM_PARTIALS], *p = ps;
1190	n/a	double xsave, special_sum = 0.0, inf_sum = 0.0;
1191	n/a	volatile double hi, yr, lo;
1192	n/a
1193	n/a	iter = PyObject_GetIter(seq);
1194	n/a	if (iter == NULL)
1195	n/a	return NULL;
1196	n/a
1197	n/a	PyFPE_START_PROTECT("fsum", Py_DECREF(iter); return NULL)
1198	n/a
1199	n/a	for(;;) { /* for x in iterable */
1200	n/a	assert(0 <= n && n <= m);
1201	n/a	assert((m == NUM_PARTIALS && p == ps) \|\|
1202	n/a	(m > NUM_PARTIALS && p != NULL));
1203	n/a
1204	n/a	item = PyIter_Next(iter);
1205	n/a	if (item == NULL) {
1206	n/a	if (PyErr_Occurred())
1207	n/a	goto _fsum_error;
1208	n/a	break;
1209	n/a	}
1210	n/a	x = PyFloat_AsDouble(item);
1211	n/a	Py_DECREF(item);
1212	n/a	if (PyErr_Occurred())
1213	n/a	goto _fsum_error;
1214	n/a
1215	n/a	xsave = x;
1216	n/a	for (i = j = 0; j < n; j++) { /* for y in partials */
1217	n/a	y = p[j];
1218	n/a	if (fabs(x) < fabs(y)) {
1219	n/a	t = x; x = y; y = t;
1220	n/a	}
1221	n/a	hi = x + y;
1222	n/a	yr = hi - x;
1223	n/a	lo = y - yr;
1224	n/a	if (lo != 0.0)
1225	n/a	p[i++] = lo;
1226	n/a	x = hi;
1227	n/a	}
1228	n/a
1229	n/a	n = i; /* ps[i:] = [x] */
1230	n/a	if (x != 0.0) {
1231	n/a	if (! Py_IS_FINITE(x)) {
1232	n/a	/* a nonfinite x could arise either as
1233	n/a	a result of intermediate overflow, or
1234	n/a	as a result of a nan or inf in the
1235	n/a	summands */
1236	n/a	if (Py_IS_FINITE(xsave)) {
1237	n/a	PyErr_SetString(PyExc_OverflowError,
1238	n/a	"intermediate overflow in fsum");
1239	n/a	goto _fsum_error;
1240	n/a	}
1241	n/a	if (Py_IS_INFINITY(xsave))
1242	n/a	inf_sum += xsave;
1243	n/a	special_sum += xsave;
1244	n/a	/* reset partials */
1245	n/a	n = 0;
1246	n/a	}
1247	n/a	else if (n >= m && _fsum_realloc(&p, n, ps, &m))
1248	n/a	goto _fsum_error;
1249	n/a	else
1250	n/a	p[n++] = x;
1251	n/a	}
1252	n/a	}
1253	n/a
1254	n/a	if (special_sum != 0.0) {
1255	n/a	if (Py_IS_NAN(inf_sum))
1256	n/a	PyErr_SetString(PyExc_ValueError,
1257	n/a	"-inf + inf in fsum");
1258	n/a	else
1259	n/a	sum = PyFloat_FromDouble(special_sum);
1260	n/a	goto _fsum_error;
1261	n/a	}
1262	n/a
1263	n/a	hi = 0.0;
1264	n/a	if (n > 0) {
1265	n/a	hi = p[--n];
1266	n/a	/* sum_exact(ps, hi) from the top, stop when the sum becomes
1267	n/a	inexact. */
1268	n/a	while (n > 0) {
1269	n/a	x = hi;
1270	n/a	y = p[--n];
1271	n/a	assert(fabs(y) < fabs(x));
1272	n/a	hi = x + y;
1273	n/a	yr = hi - x;
1274	n/a	lo = y - yr;
1275	n/a	if (lo != 0.0)
1276	n/a	break;
1277	n/a	}
1278	n/a	/* Make half-even rounding work across multiple partials.
1279	n/a	Needed so that sum([1e-16, 1, 1e16]) will round-up the last
1280	n/a	digit to two instead of down to zero (the 1e-16 makes the 1
1281	n/a	slightly closer to two). With a potential 1 ULP rounding
1282	n/a	error fixed-up, math.fsum() can guarantee commutativity. */
1283	n/a	if (n > 0 && ((lo < 0.0 && p[n-1] < 0.0) \|\|
1284	n/a	(lo > 0.0 && p[n-1] > 0.0))) {
1285	n/a	y = lo * 2.0;
1286	n/a	x = hi + y;
1287	n/a	yr = x - hi;
1288	n/a	if (y == yr)
1289	n/a	hi = x;
1290	n/a	}
1291	n/a	}
1292	n/a	sum = PyFloat_FromDouble(hi);
1293	n/a
1294	n/a	_fsum_error:
1295	n/a	PyFPE_END_PROTECT(hi)
1296	n/a	Py_DECREF(iter);
1297	n/a	if (p != ps)
1298	n/a	PyMem_Free(p);
1299	n/a	return sum;
1300	n/a	}
1301	n/a
1302	n/a	#undef NUM_PARTIALS
1303	n/a
1304	n/a
1305	n/a	/* Return the smallest integer k such that n < 2**k, or 0 if n == 0.
1306	n/a	* Equivalent to floor(lg(x))+1. Also equivalent to: bitwidth_of_type -
1307	n/a	* count_leading_zero_bits(x)
1308	n/a	*/
1309	n/a
1310	n/a	/* XXX: This routine does more or less the same thing as
1311	n/a	* bits_in_digit() in Objects/longobject.c. Someday it would be nice to
1312	n/a	* consolidate them. On BSD, there's a library function called fls()
1313	n/a	* that we could use, and GCC provides __builtin_clz().
1314	n/a	*/
1315	n/a
1316	n/a	static unsigned long
1317	n/a	bit_length(unsigned long n)
1318	n/a	{
1319	n/a	unsigned long len = 0;
1320	n/a	while (n != 0) {
1321	n/a	++len;
1322	n/a	n >>= 1;
1323	n/a	}
1324	n/a	return len;
1325	n/a	}
1326	n/a
1327	n/a	static unsigned long
1328	n/a	count_set_bits(unsigned long n)
1329	n/a	{
1330	n/a	unsigned long count = 0;
1331	n/a	while (n != 0) {
1332	n/a	++count;
1333	n/a	n &= n - 1; /* clear least significant bit */
1334	n/a	}
1335	n/a	return count;
1336	n/a	}
1337	n/a
1338	n/a	/* Divide-and-conquer factorial algorithm
1339	n/a	*
1340	n/a	* Based on the formula and pseudo-code provided at:
1341	n/a	* http://www.luschny.de/math/factorial/binarysplitfact.html
1342	n/a	*
1343	n/a	* Faster algorithms exist, but they're more complicated and depend on
1344	n/a	* a fast prime factorization algorithm.
1345	n/a	*
1346	n/a	* Notes on the algorithm
1347	n/a	* ----------------------
1348	n/a	*
1349	n/a	* factorial(n) is written in the form 2*k m, with m odd. k and m are
1350	n/a	* computed separately, and then combined using a left shift.
1351	n/a	*
1352	n/a	* The function factorial_odd_part computes the odd part m (i.e., the greatest
1353	n/a	* odd divisor) of factorial(n), using the formula:
1354	n/a	*
1355	n/a	* factorial_odd_part(n) =
1356	n/a	*
1357	n/a	* product_{i >= 0} product_{0 < j <= n / 2**i, j odd} j
1358	n/a	*
1359	n/a	* Example: factorial_odd_part(20) =
1360	n/a	*
1361	n/a	* (1) *
1362	n/a	* (1) *
1363	n/a	* (1 * 3 * 5) *
1364	n/a	* (1 * 3 * 5 * 7 * 9)
1365	n/a	* (1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19)
1366	n/a	*
1367	n/a	* Here i goes from large to small: the first term corresponds to i=4 (any
1368	n/a	* larger i gives an empty product), and the last term corresponds to i=0.
1369	n/a	* Each term can be computed from the last by multiplying by the extra odd
1370	n/a	* numbers required: e.g., to get from the penultimate term to the last one,
1371	n/a	* we multiply by (11 * 13 * 15 * 17 * 19).
1372	n/a	*
1373	n/a	* To see a hint of why this formula works, here are the same numbers as above
1374	n/a	* but with the even parts (i.e., the appropriate powers of 2) included. For
1375	n/a	* each subterm in the product for i, we multiply that subterm by 2**i:
1376	n/a	*
1377	n/a	* factorial(20) =
1378	n/a	*
1379	n/a	* (16) *
1380	n/a	* (8) *
1381	n/a	* (4 * 12 * 20) *
1382	n/a	* (2 * 6 * 10 * 14 * 18) *
1383	n/a	* (1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19)
1384	n/a	*
1385	n/a	* The factorial_partial_product function computes the product of all odd j in
1386	n/a	* range(start, stop) for given start and stop. It's used to compute the
1387	n/a	* partial products like (11 * 13 * 15 * 17 * 19) in the example above. It
1388	n/a	* operates recursively, repeatedly splitting the range into two roughly equal
1389	n/a	* pieces until the subranges are small enough to be computed using only C
1390	n/a	* integer arithmetic.
1391	n/a	*
1392	n/a	* The two-valuation k (i.e., the exponent of the largest power of 2 dividing
1393	n/a	* the factorial) is computed independently in the main math_factorial
1394	n/a	* function. By standard results, its value is:
1395	n/a	*
1396	n/a	* two_valuation = n//2 + n//4 + n//8 + ....
1397	n/a	*
1398	n/a	* It can be shown (e.g., by complete induction on n) that two_valuation is
1399	n/a	* equal to n - count_set_bits(n), where count_set_bits(n) gives the number of
1400	n/a	* '1'-bits in the binary expansion of n.
1401	n/a	*/
1402	n/a
1403	n/a	/* factorial_partial_product: Compute product(range(start, stop, 2)) using
1404	n/a	* divide and conquer. Assumes start and stop are odd and stop > start.
1405	n/a	* max_bits must be >= bit_length(stop - 2). */
1406	n/a
1407	n/a	static PyObject *
1408	n/a	factorial_partial_product(unsigned long start, unsigned long stop,
1409	n/a	unsigned long max_bits)
1410	n/a	{
1411	n/a	unsigned long midpoint, num_operands;
1412	n/a	PyObject left = NULL, right = NULL, *result = NULL;
1413	n/a
1414	n/a	/* If the return value will fit an unsigned long, then we can
1415	n/a	* multiply in a tight, fast loop where each multiply is O(1).
1416	n/a	* Compute an upper bound on the number of bits required to store
1417	n/a	* the answer.
1418	n/a	*
1419	n/a	* Storing some integer z requires floor(lg(z))+1 bits, which is
1420	n/a	* conveniently the value returned by bit_length(z). The
1421	n/a	* product x*y will require at most
1422	n/a	* bit_length(x) + bit_length(y) bits to store, based
1423	n/a	* on the idea that lg product = lg x + lg y.
1424	n/a	*
1425	n/a	* We know that stop - 2 is the largest number to be multiplied. From
1426	n/a	* there, we have: bit_length(answer) <= num_operands *
1427	n/a	* bit_length(stop - 2)
1428	n/a	*/
1429	n/a
1430	n/a	num_operands = (stop - start) / 2;
1431	n/a	/* The "num_operands <= 8 * SIZEOF_LONG" check guards against the
1432	n/a	* unlikely case of an overflow in num_operands * max_bits. */
1433	n/a	if (num_operands <= 8 * SIZEOF_LONG &&
1434	n/a	num_operands * max_bits <= 8 * SIZEOF_LONG) {
1435	n/a	unsigned long j, total;
1436	n/a	for (total = start, j = start + 2; j < stop; j += 2)
1437	n/a	total *= j;
1438	n/a	return PyLong_FromUnsignedLong(total);
1439	n/a	}
1440	n/a
1441	n/a	/* find midpoint of range(start, stop), rounded up to next odd number. */
1442	n/a	midpoint = (start + num_operands) \| 1;
1443	n/a	left = factorial_partial_product(start, midpoint,
1444	n/a	bit_length(midpoint - 2));
1445	n/a	if (left == NULL)
1446	n/a	goto error;
1447	n/a	right = factorial_partial_product(midpoint, stop, max_bits);
1448	n/a	if (right == NULL)
1449	n/a	goto error;
1450	n/a	result = PyNumber_Multiply(left, right);
1451	n/a
1452	n/a	error:
1453	n/a	Py_XDECREF(left);
1454	n/a	Py_XDECREF(right);
1455	n/a	return result;
1456	n/a	}
1457	n/a
1458	n/a	/* factorial_odd_part: compute the odd part of factorial(n). */
1459	n/a
1460	n/a	static PyObject *
1461	n/a	factorial_odd_part(unsigned long n)
1462	n/a	{
1463	n/a	long i;
1464	n/a	unsigned long v, lower, upper;
1465	n/a	PyObject partial, tmp, inner, outer;
1466	n/a
1467	n/a	inner = PyLong_FromLong(1);
1468	n/a	if (inner == NULL)
1469	n/a	return NULL;
1470	n/a	outer = inner;
1471	n/a	Py_INCREF(outer);
1472	n/a
1473	n/a	upper = 3;
1474	n/a	for (i = bit_length(n) - 2; i >= 0; i--) {
1475	n/a	v = n >> i;
1476	n/a	if (v <= 2)
1477	n/a	continue;
1478	n/a	lower = upper;
1479	n/a	/* (v + 1) \| 1 = least odd integer strictly larger than n / 2*i /
1480	n/a	upper = (v + 1) \| 1;
1481	n/a	/* Here inner is the product of all odd integers j in the range (0,
1482	n/a	n/2**(i+1)]. The factorial_partial_product call below gives the
1483	n/a	product of all odd integers j in the range (n/2(i+1), n/2i]. */
1484	n/a	partial = factorial_partial_product(lower, upper, bit_length(upper-2));
1485	n/a	/* inner = partial /
1486	n/a	if (partial == NULL)
1487	n/a	goto error;
1488	n/a	tmp = PyNumber_Multiply(inner, partial);
1489	n/a	Py_DECREF(partial);
1490	n/a	if (tmp == NULL)
1491	n/a	goto error;
1492	n/a	Py_DECREF(inner);
1493	n/a	inner = tmp;
1494	n/a	/* Now inner is the product of all odd integers j in the range (0,
1495	n/a	n/2*i], giving the inner product in the formula above. /
1496	n/a
1497	n/a	/* outer = inner; /
1498	n/a	tmp = PyNumber_Multiply(outer, inner);
1499	n/a	if (tmp == NULL)
1500	n/a	goto error;
1501	n/a	Py_DECREF(outer);
1502	n/a	outer = tmp;
1503	n/a	}
1504	n/a	Py_DECREF(inner);
1505	n/a	return outer;
1506	n/a
1507	n/a	error:
1508	n/a	Py_DECREF(outer);
1509	n/a	Py_DECREF(inner);
1510	n/a	return NULL;
1511	n/a	}
1512	n/a
1513	n/a
1514	n/a	/* Lookup table for small factorial values */
1515	n/a
1516	n/a	static const unsigned long SmallFactorials[] = {
1517	n/a	1, 1, 2, 6, 24, 120, 720, 5040, 40320,
1518	n/a	362880, 3628800, 39916800, 479001600,
1519	n/a	#if SIZEOF_LONG >= 8
1520	n/a	6227020800, 87178291200, 1307674368000,
1521	n/a	20922789888000, 355687428096000, 6402373705728000,
1522	n/a	121645100408832000, 2432902008176640000
1523	n/a	#endif
1524	n/a	};
1525	n/a
1526	n/a	/*[clinic input]
1527	n/a	math.factorial
1528	n/a
1529	n/a	x as arg: object
1530	n/a	/
1531	n/a
1532	n/a	Find x!.
1533	n/a
1534	n/a	Raise a ValueError if x is negative or non-integral.
1535	n/a	[clinic start generated code]*/
1536	n/a
1537	n/a	static PyObject *
1538	n/a	math_factorial(PyObject module, PyObject arg)
1539	n/a	/[clinic end generated code: output=6686f26fae00e9ca input=6d1c8105c0d91fb4]/
1540	n/a	{
1541	n/a	long x;
1542	n/a	int overflow;
1543	n/a	PyObject result, odd_part, *two_valuation;
1544	n/a
1545	n/a	if (PyFloat_Check(arg)) {
1546	n/a	PyObject *lx;
1547	n/a	double dx = PyFloat_AS_DOUBLE((PyFloatObject *)arg);
1548	n/a	if (!(Py_IS_FINITE(dx) && dx == floor(dx))) {
1549	n/a	PyErr_SetString(PyExc_ValueError,
1550	n/a	"factorial() only accepts integral values");
1551	n/a	return NULL;
1552	n/a	}
1553	n/a	lx = PyLong_FromDouble(dx);
1554	n/a	if (lx == NULL)
1555	n/a	return NULL;
1556	n/a	x = PyLong_AsLongAndOverflow(lx, &overflow);
1557	n/a	Py_DECREF(lx);
1558	n/a	}
1559	n/a	else
1560	n/a	x = PyLong_AsLongAndOverflow(arg, &overflow);
1561	n/a
1562	n/a	if (x == -1 && PyErr_Occurred()) {
1563	n/a	return NULL;
1564	n/a	}
1565	n/a	else if (overflow == 1) {
1566	n/a	PyErr_Format(PyExc_OverflowError,
1567	n/a	"factorial() argument should not exceed %ld",
1568	n/a	LONG_MAX);
1569	n/a	return NULL;
1570	n/a	}
1571	n/a	else if (overflow == -1 \|\| x < 0) {
1572	n/a	PyErr_SetString(PyExc_ValueError,
1573	n/a	"factorial() not defined for negative values");
1574	n/a	return NULL;
1575	n/a	}
1576	n/a
1577	n/a	/* use lookup table if x is small */
1578	n/a	if (x < (long)Py_ARRAY_LENGTH(SmallFactorials))
1579	n/a	return PyLong_FromUnsignedLong(SmallFactorials[x]);
1580	n/a
1581	n/a	/* else express in the form odd_part * 2**two_valuation, and compute as
1582	n/a	odd_part << two_valuation. */
1583	n/a	odd_part = factorial_odd_part(x);
1584	n/a	if (odd_part == NULL)
1585	n/a	return NULL;
1586	n/a	two_valuation = PyLong_FromLong(x - count_set_bits(x));
1587	n/a	if (two_valuation == NULL) {
1588	n/a	Py_DECREF(odd_part);
1589	n/a	return NULL;
1590	n/a	}
1591	n/a	result = PyNumber_Lshift(odd_part, two_valuation);
1592	n/a	Py_DECREF(two_valuation);
1593	n/a	Py_DECREF(odd_part);
1594	n/a	return result;
1595	n/a	}
1596	n/a
1597	n/a
1598	n/a	/*[clinic input]
1599	n/a	math.trunc
1600	n/a
1601	n/a	x: object
1602	n/a	/
1603	n/a
1604	n/a	Truncates the Real x to the nearest Integral toward 0.
1605	n/a
1606	n/a	Uses the __trunc__ magic method.
1607	n/a	[clinic start generated code]*/
1608	n/a
1609	n/a	static PyObject *
1610	n/a	math_trunc(PyObject module, PyObject x)
1611	n/a	/[clinic end generated code: output=34b9697b707e1031 input=2168b34e0a09134d]/
1612	n/a	{
1613	n/a	_Py_IDENTIFIER(__trunc__);
1614	n/a	PyObject trunc, result;
1615	n/a
1616	n/a	if (Py_TYPE(x)->tp_dict == NULL) {
1617	n/a	if (PyType_Ready(Py_TYPE(x)) < 0)
1618	n/a	return NULL;
1619	n/a	}
1620	n/a
1621	n/a	trunc = _PyObject_LookupSpecial(x, &PyId___trunc__);
1622	n/a	if (trunc == NULL) {
1623	n/a	if (!PyErr_Occurred())
1624	n/a	PyErr_Format(PyExc_TypeError,
1625	n/a	"type %.100s doesn't define __trunc__ method",
1626	n/a	Py_TYPE(x)->tp_name);
1627	n/a	return NULL;
1628	n/a	}
1629	n/a	result = _PyObject_CallNoArg(trunc);
1630	n/a	Py_DECREF(trunc);
1631	n/a	return result;
1632	n/a	}
1633	n/a
1634	n/a
1635	n/a	/*[clinic input]
1636	n/a	math.frexp
1637	n/a
1638	n/a	x: double
1639	n/a	/
1640	n/a
1641	n/a	Return the mantissa and exponent of x, as pair (m, e).
1642	n/a
1643	n/a	m is a float and e is an int, such that x = m * 2.**e.
1644	n/a	If x is 0, m and e are both 0. Else 0.5 <= abs(m) < 1.0.
1645	n/a	[clinic start generated code]*/
1646	n/a
1647	n/a	static PyObject *
1648	n/a	math_frexp_impl(PyObject *module, double x)
1649	n/a	/[clinic end generated code: output=03e30d252a15ad4a input=96251c9e208bc6e9]/
1650	n/a	{
1651	n/a	int i;
1652	n/a	/* deal with special cases directly, to sidestep platform
1653	n/a	differences */
1654	n/a	if (Py_IS_NAN(x) \|\| Py_IS_INFINITY(x) \|\| !x) {
1655	n/a	i = 0;
1656	n/a	}
1657	n/a	else {
1658	n/a	PyFPE_START_PROTECT("in math_frexp", return 0);
1659	n/a	x = frexp(x, &i);
1660	n/a	PyFPE_END_PROTECT(x);
1661	n/a	}
1662	n/a	return Py_BuildValue("(di)", x, i);
1663	n/a	}
1664	n/a
1665	n/a
1666	n/a	/*[clinic input]
1667	n/a	math.ldexp
1668	n/a
1669	n/a	x: double
1670	n/a	i: object
1671	n/a	/
1672	n/a
1673	n/a	Return x * (2**i).
1674	n/a
1675	n/a	This is essentially the inverse of frexp().
1676	n/a	[clinic start generated code]*/
1677	n/a
1678	n/a	static PyObject *
1679	n/a	math_ldexp_impl(PyObject module, double x, PyObject i)
1680	n/a	/[clinic end generated code: output=b6892f3c2df9cc6a input=17d5970c1a40a8c1]/
1681	n/a	{
1682	n/a	double r;
1683	n/a	long exp;
1684	n/a	int overflow;
1685	n/a
1686	n/a	if (PyLong_Check(i)) {
1687	n/a	/* on overflow, replace exponent with either LONG_MAX
1688	n/a	or LONG_MIN, depending on the sign. */
1689	n/a	exp = PyLong_AsLongAndOverflow(i, &overflow);
1690	n/a	if (exp == -1 && PyErr_Occurred())
1691	n/a	return NULL;
1692	n/a	if (overflow)
1693	n/a	exp = overflow < 0 ? LONG_MIN : LONG_MAX;
1694	n/a	}
1695	n/a	else {
1696	n/a	PyErr_SetString(PyExc_TypeError,
1697	n/a	"Expected an int as second argument to ldexp.");
1698	n/a	return NULL;
1699	n/a	}
1700	n/a
1701	n/a	if (x == 0. \|\| !Py_IS_FINITE(x)) {
1702	n/a	/* NaNs, zeros and infinities are returned unchanged */
1703	n/a	r = x;
1704	n/a	errno = 0;
1705	n/a	} else if (exp > INT_MAX) {
1706	n/a	/* overflow */
1707	n/a	r = copysign(Py_HUGE_VAL, x);
1708	n/a	errno = ERANGE;
1709	n/a	} else if (exp < INT_MIN) {
1710	n/a	/* underflow to +-0 */
1711	n/a	r = copysign(0., x);
1712	n/a	errno = 0;
1713	n/a	} else {
1714	n/a	errno = 0;
1715	n/a	PyFPE_START_PROTECT("in math_ldexp", return 0);
1716	n/a	r = ldexp(x, (int)exp);
1717	n/a	PyFPE_END_PROTECT(r);
1718	n/a	if (Py_IS_INFINITY(r))
1719	n/a	errno = ERANGE;
1720	n/a	}
1721	n/a
1722	n/a	if (errno && is_error(r))
1723	n/a	return NULL;
1724	n/a	return PyFloat_FromDouble(r);
1725	n/a	}
1726	n/a
1727	n/a
1728	n/a	/*[clinic input]
1729	n/a	math.modf
1730	n/a
1731	n/a	x: double
1732	n/a	/
1733	n/a
1734	n/a	Return the fractional and integer parts of x.
1735	n/a
1736	n/a	Both results carry the sign of x and are floats.
1737	n/a	[clinic start generated code]*/
1738	n/a
1739	n/a	static PyObject *
1740	n/a	math_modf_impl(PyObject *module, double x)
1741	n/a	/[clinic end generated code: output=90cee0260014c3c0 input=b4cfb6786afd9035]/
1742	n/a	{
1743	n/a	double y;
1744	n/a	/* some platforms don't do the right thing for NaNs and
1745	n/a	infinities, so we take care of special cases directly. */
1746	n/a	if (!Py_IS_FINITE(x)) {
1747	n/a	if (Py_IS_INFINITY(x))
1748	n/a	return Py_BuildValue("(dd)", copysign(0., x), x);
1749	n/a	else if (Py_IS_NAN(x))
1750	n/a	return Py_BuildValue("(dd)", x, x);
1751	n/a	}
1752	n/a
1753	n/a	errno = 0;
1754	n/a	PyFPE_START_PROTECT("in math_modf", return 0);
1755	n/a	x = modf(x, &y);
1756	n/a	PyFPE_END_PROTECT(x);
1757	n/a	return Py_BuildValue("(dd)", x, y);
1758	n/a	}
1759	n/a
1760	n/a
1761	n/a	/* A decent logarithm is easy to compute even for huge ints, but libm can't
1762	n/a	do that by itself -- loghelper can. func is log or log10, and name is
1763	n/a	"log" or "log10". Note that overflow of the result isn't possible: an int
1764	n/a	can contain no more than INT_MAX * SHIFT bits, so has value certainly less
1765	n/a	than 2(264 * 216) == 2280, and log2 of that is 280, which is
1766	n/a	small enough to fit in an IEEE single. log and log10 are even smaller.
1767	n/a	However, intermediate overflow is possible for an int if the number of bits
1768	n/a	in that int is larger than PY_SSIZE_T_MAX. */
1769	n/a
1770	n/a	static PyObject*
1771	n/a	loghelper(PyObject* arg, double (func)(double), const char funcname)
1772	n/a	{
1773	n/a	/* If it is int, do it ourselves. */
1774	n/a	if (PyLong_Check(arg)) {
1775	n/a	double x, result;
1776	n/a	Py_ssize_t e;
1777	n/a
1778	n/a	/* Negative or zero inputs give a ValueError. */
1779	n/a	if (Py_SIZE(arg) <= 0) {
1780	n/a	PyErr_SetString(PyExc_ValueError,
1781	n/a	"math domain error");
1782	n/a	return NULL;
1783	n/a	}
1784	n/a
1785	n/a	x = PyLong_AsDouble(arg);
1786	n/a	if (x == -1.0 && PyErr_Occurred()) {
1787	n/a	if (!PyErr_ExceptionMatches(PyExc_OverflowError))
1788	n/a	return NULL;
1789	n/a	/* Here the conversion to double overflowed, but it's possible
1790	n/a	to compute the log anyway. Clear the exception and continue. */
1791	n/a	PyErr_Clear();
1792	n/a	x = _PyLong_Frexp((PyLongObject *)arg, &e);
1793	n/a	if (x == -1.0 && PyErr_Occurred())
1794	n/a	return NULL;
1795	n/a	/* Value is ~= x * 2*e, so the log ~= log(x) + log(2) e. */
1796	n/a	result = func(x) + func(2.0) * e;
1797	n/a	}
1798	n/a	else
1799	n/a	/* Successfully converted x to a double. */
1800	n/a	result = func(x);
1801	n/a	return PyFloat_FromDouble(result);
1802	n/a	}
1803	n/a
1804	n/a	/* Else let libm handle it by itself. */
1805	n/a	return math_1(arg, func, 0);
1806	n/a	}
1807	n/a
1808	n/a
1809	n/a	/*[clinic input]
1810	n/a	math.log
1811	n/a
1812	n/a	x: object
1813	n/a	[
1814	n/a	base: object(c_default="NULL") = math.e
1815	n/a	]
1816	n/a	/
1817	n/a
1818	n/a	Return the logarithm of x to the given base.
1819	n/a
1820	n/a	If the base not specified, returns the natural logarithm (base e) of x.
1821	n/a	[clinic start generated code]*/
1822	n/a
1823	n/a	static PyObject *
1824	n/a	math_log_impl(PyObject module, PyObject x, int group_right_1,
1825	n/a	PyObject *base)
1826	n/a	/[clinic end generated code: output=7b5a39e526b73fc9 input=0f62d5726cbfebbd]/
1827	n/a	{
1828	n/a	PyObject num, den;
1829	n/a	PyObject *ans;
1830	n/a
1831	n/a	num = loghelper(x, m_log, "log");
1832	n/a	if (num == NULL \|\| base == NULL)
1833	n/a	return num;
1834	n/a
1835	n/a	den = loghelper(base, m_log, "log");
1836	n/a	if (den == NULL) {
1837	n/a	Py_DECREF(num);
1838	n/a	return NULL;
1839	n/a	}
1840	n/a
1841	n/a	ans = PyNumber_TrueDivide(num, den);
1842	n/a	Py_DECREF(num);
1843	n/a	Py_DECREF(den);
1844	n/a	return ans;
1845	n/a	}
1846	n/a
1847	n/a
1848	n/a	/*[clinic input]
1849	n/a	math.log2
1850	n/a
1851	n/a	x: object
1852	n/a	/
1853	n/a
1854	n/a	Return the base 2 logarithm of x.
1855	n/a	[clinic start generated code]*/
1856	n/a
1857	n/a	static PyObject *
1858	n/a	math_log2(PyObject module, PyObject x)
1859	n/a	/[clinic end generated code: output=5425899a4d5d6acb input=08321262bae4f39b]/
1860	n/a	{
1861	n/a	return loghelper(x, m_log2, "log2");
1862	n/a	}
1863	n/a
1864	n/a
1865	n/a	/*[clinic input]
1866	n/a	math.log10
1867	n/a
1868	n/a	x: object
1869	n/a	/
1870	n/a
1871	n/a	Return the base 10 logarithm of x.
1872	n/a	[clinic start generated code]*/
1873	n/a
1874	n/a	static PyObject *
1875	n/a	math_log10(PyObject module, PyObject x)
1876	n/a	/[clinic end generated code: output=be72a64617df9c6f input=b2469d02c6469e53]/
1877	n/a	{
1878	n/a	return loghelper(x, m_log10, "log10");
1879	n/a	}
1880	n/a
1881	n/a
1882	n/a	/*[clinic input]
1883	n/a	math.fmod
1884	n/a
1885	n/a	x: double
1886	n/a	y: double
1887	n/a	/
1888	n/a
1889	n/a	Return fmod(x, y), according to platform C.
1890	n/a
1891	n/a	x % y may differ.
1892	n/a	[clinic start generated code]*/
1893	n/a
1894	n/a	static PyObject *
1895	n/a	math_fmod_impl(PyObject *module, double x, double y)
1896	n/a	/[clinic end generated code: output=7559d794343a27b5 input=4f84caa8cfc26a03]/
1897	n/a	{
1898	n/a	double r;
1899	n/a	/* fmod(x, +/-Inf) returns x for finite x. */
1900	n/a	if (Py_IS_INFINITY(y) && Py_IS_FINITE(x))
1901	n/a	return PyFloat_FromDouble(x);
1902	n/a	errno = 0;
1903	n/a	PyFPE_START_PROTECT("in math_fmod", return 0);
1904	n/a	r = fmod(x, y);
1905	n/a	PyFPE_END_PROTECT(r);
1906	n/a	if (Py_IS_NAN(r)) {
1907	n/a	if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
1908	n/a	errno = EDOM;
1909	n/a	else
1910	n/a	errno = 0;
1911	n/a	}
1912	n/a	if (errno && is_error(r))
1913	n/a	return NULL;
1914	n/a	else
1915	n/a	return PyFloat_FromDouble(r);
1916	n/a	}
1917	n/a
1918	n/a
1919	n/a	/*[clinic input]
1920	n/a	math.hypot
1921	n/a
1922	n/a	x: double
1923	n/a	y: double
1924	n/a	/
1925	n/a
1926	n/a	Return the Euclidean distance, sqrt(xx + yy).
1927	n/a	[clinic start generated code]*/
1928	n/a
1929	n/a	static PyObject *
1930	n/a	math_hypot_impl(PyObject *module, double x, double y)
1931	n/a	/[clinic end generated code: output=b7686e5be468ef87 input=7f8eea70406474aa]/
1932	n/a	{
1933	n/a	double r;
1934	n/a	/* hypot(x, +/-Inf) returns Inf, even if x is a NaN. */
1935	n/a	if (Py_IS_INFINITY(x))
1936	n/a	return PyFloat_FromDouble(fabs(x));
1937	n/a	if (Py_IS_INFINITY(y))
1938	n/a	return PyFloat_FromDouble(fabs(y));
1939	n/a	errno = 0;
1940	n/a	PyFPE_START_PROTECT("in math_hypot", return 0);
1941	n/a	r = hypot(x, y);
1942	n/a	PyFPE_END_PROTECT(r);
1943	n/a	if (Py_IS_NAN(r)) {
1944	n/a	if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
1945	n/a	errno = EDOM;
1946	n/a	else
1947	n/a	errno = 0;
1948	n/a	}
1949	n/a	else if (Py_IS_INFINITY(r)) {
1950	n/a	if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
1951	n/a	errno = ERANGE;
1952	n/a	else
1953	n/a	errno = 0;
1954	n/a	}
1955	n/a	if (errno && is_error(r))
1956	n/a	return NULL;
1957	n/a	else
1958	n/a	return PyFloat_FromDouble(r);
1959	n/a	}
1960	n/a
1961	n/a
1962	n/a	/* pow can't use math_2, but needs its own wrapper: the problem is
1963	n/a	that an infinite result can arise either as a result of overflow
1964	n/a	(in which case OverflowError should be raised) or as a result of
1965	n/a	e.g. 0.**-5. (for which ValueError needs to be raised.)
1966	n/a	*/
1967	n/a
1968	n/a	/*[clinic input]
1969	n/a	math.pow
1970	n/a
1971	n/a	x: double
1972	n/a	y: double
1973	n/a	/
1974	n/a
1975	n/a	Return x**y (x to the power of y).
1976	n/a	[clinic start generated code]*/
1977	n/a
1978	n/a	static PyObject *
1979	n/a	math_pow_impl(PyObject *module, double x, double y)
1980	n/a	/[clinic end generated code: output=fff93e65abccd6b0 input=c26f1f6075088bfd]/
1981	n/a	{
1982	n/a	double r;
1983	n/a	int odd_y;
1984	n/a
1985	n/a	/* deal directly with IEEE specials, to cope with problems on various
1986	n/a	platforms whose semantics don't exactly match C99 */
1987	n/a	r = 0.; /* silence compiler warning */
1988	n/a	if (!Py_IS_FINITE(x) \|\| !Py_IS_FINITE(y)) {
1989	n/a	errno = 0;
1990	n/a	if (Py_IS_NAN(x))
1991	n/a	r = y == 0. ? 1. : x; /* NaN*0 = 1 /
1992	n/a	else if (Py_IS_NAN(y))
1993	n/a	r = x == 1. ? 1. : y; /* 1*NaN = 1 /
1994	n/a	else if (Py_IS_INFINITY(x)) {
1995	n/a	odd_y = Py_IS_FINITE(y) && fmod(fabs(y), 2.0) == 1.0;
1996	n/a	if (y > 0.)
1997	n/a	r = odd_y ? x : fabs(x);
1998	n/a	else if (y == 0.)
1999	n/a	r = 1.;
2000	n/a	else /* y < 0. */
2001	n/a	r = odd_y ? copysign(0., x) : 0.;
2002	n/a	}
2003	n/a	else if (Py_IS_INFINITY(y)) {
2004	n/a	if (fabs(x) == 1.0)
2005	n/a	r = 1.;
2006	n/a	else if (y > 0. && fabs(x) > 1.0)
2007	n/a	r = y;
2008	n/a	else if (y < 0. && fabs(x) < 1.0) {
2009	n/a	r = -y; /* result is +inf */
2010	n/a	if (x == 0.) /* 0*-inf: divide-by-zero /
2011	n/a	errno = EDOM;
2012	n/a	}
2013	n/a	else
2014	n/a	r = 0.;
2015	n/a	}
2016	n/a	}
2017	n/a	else {
2018	n/a	/* let libm handle finite*finite /
2019	n/a	errno = 0;
2020	n/a	PyFPE_START_PROTECT("in math_pow", return 0);
2021	n/a	r = pow(x, y);
2022	n/a	PyFPE_END_PROTECT(r);
2023	n/a	/* a NaN result should arise only from (-ve)**(finite
2024	n/a	non-integer); in this case we want to raise ValueError. */
2025	n/a	if (!Py_IS_FINITE(r)) {
2026	n/a	if (Py_IS_NAN(r)) {
2027	n/a	errno = EDOM;
2028	n/a	}
2029	n/a	/*
2030	n/a	an infinite result here arises either from:
2031	n/a	(A) (+/-0.)**negative (-> divide-by-zero)
2032	n/a	(B) overflow of x**y with x and y finite
2033	n/a	*/
2034	n/a	else if (Py_IS_INFINITY(r)) {
2035	n/a	if (x == 0.)
2036	n/a	errno = EDOM;
2037	n/a	else
2038	n/a	errno = ERANGE;
2039	n/a	}
2040	n/a	}
2041	n/a	}
2042	n/a
2043	n/a	if (errno && is_error(r))
2044	n/a	return NULL;
2045	n/a	else
2046	n/a	return PyFloat_FromDouble(r);
2047	n/a	}
2048	n/a
2049	n/a
2050	n/a	static const double degToRad = Py_MATH_PI / 180.0;
2051	n/a	static const double radToDeg = 180.0 / Py_MATH_PI;
2052	n/a
2053	n/a	/*[clinic input]
2054	n/a	math.degrees
2055	n/a
2056	n/a	x: double
2057	n/a	/
2058	n/a
2059	n/a	Convert angle x from radians to degrees.
2060	n/a	[clinic start generated code]*/
2061	n/a
2062	n/a	static PyObject *
2063	n/a	math_degrees_impl(PyObject *module, double x)
2064	n/a	/[clinic end generated code: output=7fea78b294acd12f input=81e016555d6e3660]/
2065	n/a	{
2066	n/a	return PyFloat_FromDouble(x * radToDeg);
2067	n/a	}
2068	n/a
2069	n/a
2070	n/a	/*[clinic input]
2071	n/a	math.radians
2072	n/a
2073	n/a	x: double
2074	n/a	/
2075	n/a
2076	n/a	Convert angle x from degrees to radians.
2077	n/a	[clinic start generated code]*/
2078	n/a
2079	n/a	static PyObject *
2080	n/a	math_radians_impl(PyObject *module, double x)
2081	n/a	/[clinic end generated code: output=34daa47caf9b1590 input=91626fc489fe3d63]/
2082	n/a	{
2083	n/a	return PyFloat_FromDouble(x * degToRad);
2084	n/a	}
2085	n/a
2086	n/a
2087	n/a	/*[clinic input]
2088	n/a	math.isfinite
2089	n/a
2090	n/a	x: double
2091	n/a	/
2092	n/a
2093	n/a	Return True if x is neither an infinity nor a NaN, and False otherwise.
2094	n/a	[clinic start generated code]*/
2095	n/a
2096	n/a	static PyObject *
2097	n/a	math_isfinite_impl(PyObject *module, double x)
2098	n/a	/[clinic end generated code: output=8ba1f396440c9901 input=46967d254812e54a]/
2099	n/a	{
2100	n/a	return PyBool_FromLong((long)Py_IS_FINITE(x));
2101	n/a	}
2102	n/a
2103	n/a
2104	n/a	/*[clinic input]
2105	n/a	math.isnan
2106	n/a
2107	n/a	x: double
2108	n/a	/
2109	n/a
2110	n/a	Return True if x is a NaN (not a number), and False otherwise.
2111	n/a	[clinic start generated code]*/
2112	n/a
2113	n/a	static PyObject *
2114	n/a	math_isnan_impl(PyObject *module, double x)
2115	n/a	/[clinic end generated code: output=f537b4d6df878c3e input=935891e66083f46a]/
2116	n/a	{
2117	n/a	return PyBool_FromLong((long)Py_IS_NAN(x));
2118	n/a	}
2119	n/a
2120	n/a
2121	n/a	/*[clinic input]
2122	n/a	math.isinf
2123	n/a
2124	n/a	x: double
2125	n/a	/
2126	n/a
2127	n/a	Return True if x is a positive or negative infinity, and False otherwise.
2128	n/a	[clinic start generated code]*/
2129	n/a
2130	n/a	static PyObject *
2131	n/a	math_isinf_impl(PyObject *module, double x)
2132	n/a	/[clinic end generated code: output=9f00cbec4de7b06b input=32630e4212cf961f]/
2133	n/a	{
2134	n/a	return PyBool_FromLong((long)Py_IS_INFINITY(x));
2135	n/a	}
2136	n/a
2137	n/a
2138	n/a	/*[clinic input]
2139	n/a	math.isclose -> bool
2140	n/a
2141	n/a	a: double
2142	n/a	b: double
2143	n/a	*
2144	n/a	rel_tol: double = 1e-09
2145	n/a	maximum difference for being considered "close", relative to the
2146	n/a	magnitude of the input values
2147	n/a	abs_tol: double = 0.0
2148	n/a	maximum difference for being considered "close", regardless of the
2149	n/a	magnitude of the input values
2150	n/a
2151	n/a	Determine whether two floating point numbers are close in value.
2152	n/a
2153	n/a	Return True if a is close in value to b, and False otherwise.
2154	n/a
2155	n/a	For the values to be considered close, the difference between them
2156	n/a	must be smaller than at least one of the tolerances.
2157	n/a
2158	n/a	-inf, inf and NaN behave similarly to the IEEE 754 Standard. That
2159	n/a	is, NaN is not close to anything, even itself. inf and -inf are
2160	n/a	only close to themselves.
2161	n/a	[clinic start generated code]*/
2162	n/a
2163	n/a	static int
2164	n/a	math_isclose_impl(PyObject *module, double a, double b, double rel_tol,
2165	n/a	double abs_tol)
2166	n/a	/[clinic end generated code: output=b73070207511952d input=f28671871ea5bfba]/
2167	n/a	{
2168	n/a	double diff = 0.0;
2169	n/a
2170	n/a	/* sanity check on the inputs */
2171	n/a	if (rel_tol < 0.0 \|\| abs_tol < 0.0 ) {
2172	n/a	PyErr_SetString(PyExc_ValueError,
2173	n/a	"tolerances must be non-negative");
2174	n/a	return -1;
2175	n/a	}
2176	n/a
2177	n/a	if ( a == b ) {
2178	n/a	/* short circuit exact equality -- needed to catch two infinities of
2179	n/a	the same sign. And perhaps speeds things up a bit sometimes.
2180	n/a	*/
2181	n/a	return 1;
2182	n/a	}
2183	n/a
2184	n/a	/* This catches the case of two infinities of opposite sign, or
2185	n/a	one infinity and one finite number. Two infinities of opposite
2186	n/a	sign would otherwise have an infinite relative tolerance.
2187	n/a	Two infinities of the same sign are caught by the equality check
2188	n/a	above.
2189	n/a	*/
2190	n/a
2191	n/a	if (Py_IS_INFINITY(a) \|\| Py_IS_INFINITY(b)) {
2192	n/a	return 0;
2193	n/a	}
2194	n/a
2195	n/a	/* now do the regular computation
2196	n/a	this is essentially the "weak" test from the Boost library
2197	n/a	*/
2198	n/a
2199	n/a	diff = fabs(b - a);
2200	n/a
2201	n/a	return (((diff <= fabs(rel_tol * b)) \|\|
2202	n/a	(diff <= fabs(rel_tol * a))) \|\|
2203	n/a	(diff <= abs_tol));
2204	n/a	}
2205	n/a
2206	n/a
2207	n/a	static PyMethodDef math_methods[] = {
2208	n/a	{"acos", math_acos, METH_O, math_acos_doc},
2209	n/a	{"acosh", math_acosh, METH_O, math_acosh_doc},
2210	n/a	{"asin", math_asin, METH_O, math_asin_doc},
2211	n/a	{"asinh", math_asinh, METH_O, math_asinh_doc},
2212	n/a	{"atan", math_atan, METH_O, math_atan_doc},
2213	n/a	{"atan2", math_atan2, METH_VARARGS, math_atan2_doc},
2214	n/a	{"atanh", math_atanh, METH_O, math_atanh_doc},
2215	n/a	MATH_CEIL_METHODDEF
2216	n/a	{"copysign", math_copysign, METH_VARARGS, math_copysign_doc},
2217	n/a	{"cos", math_cos, METH_O, math_cos_doc},
2218	n/a	{"cosh", math_cosh, METH_O, math_cosh_doc},
2219	n/a	MATH_DEGREES_METHODDEF
2220	n/a	{"erf", math_erf, METH_O, math_erf_doc},
2221	n/a	{"erfc", math_erfc, METH_O, math_erfc_doc},
2222	n/a	{"exp", math_exp, METH_O, math_exp_doc},
2223	n/a	{"expm1", math_expm1, METH_O, math_expm1_doc},
2224	n/a	{"fabs", math_fabs, METH_O, math_fabs_doc},
2225	n/a	MATH_FACTORIAL_METHODDEF
2226	n/a	MATH_FLOOR_METHODDEF
2227	n/a	MATH_FMOD_METHODDEF
2228	n/a	MATH_FREXP_METHODDEF
2229	n/a	MATH_FSUM_METHODDEF
2230	n/a	{"gamma", math_gamma, METH_O, math_gamma_doc},
2231	n/a	MATH_GCD_METHODDEF
2232	n/a	MATH_HYPOT_METHODDEF
2233	n/a	MATH_ISCLOSE_METHODDEF
2234	n/a	MATH_ISFINITE_METHODDEF
2235	n/a	MATH_ISINF_METHODDEF
2236	n/a	MATH_ISNAN_METHODDEF
2237	n/a	MATH_LDEXP_METHODDEF
2238	n/a	{"lgamma", math_lgamma, METH_O, math_lgamma_doc},
2239	n/a	MATH_LOG_METHODDEF
2240	n/a	{"log1p", math_log1p, METH_O, math_log1p_doc},
2241	n/a	MATH_LOG10_METHODDEF
2242	n/a	MATH_LOG2_METHODDEF
2243	n/a	MATH_MODF_METHODDEF
2244	n/a	MATH_POW_METHODDEF
2245	n/a	MATH_RADIANS_METHODDEF
2246	n/a	{"sin", math_sin, METH_O, math_sin_doc},
2247	n/a	{"sinh", math_sinh, METH_O, math_sinh_doc},
2248	n/a	{"sqrt", math_sqrt, METH_O, math_sqrt_doc},
2249	n/a	{"tan", math_tan, METH_O, math_tan_doc},
2250	n/a	{"tanh", math_tanh, METH_O, math_tanh_doc},
2251	n/a	MATH_TRUNC_METHODDEF
2252	n/a	{NULL, NULL} /* sentinel */
2253	n/a	};
2254	n/a
2255	n/a
2256	n/a	PyDoc_STRVAR(module_doc,
2257	n/a	"This module is always available. It provides access to the\n"
2258	n/a	"mathematical functions defined by the C standard.");
2259	n/a
2260	n/a
2261	n/a	static struct PyModuleDef mathmodule = {
2262	n/a	PyModuleDef_HEAD_INIT,
2263	n/a	"math",
2264	n/a	module_doc,
2265	n/a	-1,
2266	n/a	math_methods,
2267	n/a	NULL,
2268	n/a	NULL,
2269	n/a	NULL,
2270	n/a	NULL
2271	n/a	};
2272	n/a
2273	n/a	PyMODINIT_FUNC
2274	n/a	PyInit_math(void)
2275	n/a	{
2276	n/a	PyObject *m;
2277	n/a
2278	n/a	m = PyModule_Create(&mathmodule);
2279	n/a	if (m == NULL)
2280	n/a	goto finally;
2281	n/a
2282	n/a	PyModule_AddObject(m, "pi", PyFloat_FromDouble(Py_MATH_PI));
2283	n/a	PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E));
2284	n/a	PyModule_AddObject(m, "tau", PyFloat_FromDouble(Py_MATH_TAU)); /* 2pi */
2285	n/a	PyModule_AddObject(m, "inf", PyFloat_FromDouble(m_inf()));
2286	n/a	#if !defined(PY_NO_SHORT_FLOAT_REPR) \|\| defined(Py_NAN)
2287	n/a	PyModule_AddObject(m, "nan", PyFloat_FromDouble(m_nan()));
2288	n/a	#endif
2289	n/a
2290	n/a	finally:
2291	n/a	return m;
2292	n/a	}