From 2a8c0e62f0c6d2fb4c7e202bd10852d318bc1ec6 Mon Sep 17 00:00:00 2001 From: ananthan-123 <60189117+ananthan-123@users.noreply.github.com> Date: Tue, 18 Feb 2020 17:20:04 +0530 Subject: [PATCH 01/18] Update math.rst --- Doc/library/math.rst | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/Doc/library/math.rst b/Doc/library/math.rst index c4c180037f8788..5f7cf95d5b4b59 100644 --- a/Doc/library/math.rst +++ b/Doc/library/math.rst @@ -136,6 +136,13 @@ Number-theoretic and representation functions .. versionadded:: 3.5 +.. function:: lcm(a, b) + + Return the least common multiple of integers *a* and *b*. The value of + ``lcm(a, b)`` is the smallest nonnegative integer that is a multiple of + both *a* and *b*. If either *a* or *b* is zero then ``lcm(a, b)`` is zero. + + .. function:: isclose(a, b, *, rel_tol=1e-09, abs_tol=0.0) Return ``True`` if the values *a* and *b* are close to each other and From 7b5bbf71ade421740def7391dada6a9b070836f6 Mon Sep 17 00:00:00 2001 From: ananthan-123 <60189117+ananthan-123@users.noreply.github.com> Date: Tue, 18 Feb 2020 17:21:40 +0530 Subject: [PATCH 02/18] Update math.rst --- Doc/library/math.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Doc/library/math.rst b/Doc/library/math.rst index 5f7cf95d5b4b59..b89af90a3f82f7 100644 --- a/Doc/library/math.rst +++ b/Doc/library/math.rst @@ -142,6 +142,8 @@ Number-theoretic and representation functions ``lcm(a, b)`` is the smallest nonnegative integer that is a multiple of both *a* and *b*. If either *a* or *b* is zero then ``lcm(a, b)`` is zero. + .. versionadded:: 3.5 + .. function:: isclose(a, b, *, rel_tol=1e-09, abs_tol=0.0) From 2f64f5467bde4a16d7a08ea525be1c3815508df8 Mon Sep 17 00:00:00 2001 From: ananthan-123 <60189117+ananthan-123@users.noreply.github.com> Date: Tue, 18 Feb 2020 17:25:14 +0530 Subject: [PATCH 03/18] updated whats new --- Doc/whatsnew/3.9.rst | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Doc/whatsnew/3.9.rst b/Doc/whatsnew/3.9.rst index 23f0e4306ee63a..7e495563e25ae8 100644 --- a/Doc/whatsnew/3.9.rst +++ b/Doc/whatsnew/3.9.rst @@ -216,6 +216,9 @@ import attempts. math ---- +Add :func:`math.lcm`: return the least common multiple of a and b. +(Contributed by Ananthakrishnan in :issue:`39479`.) + Add :func:`math.nextafter`: return the next floating-point value after *x* towards *y*. (Contributed by Victor Stinner in :issue:`39288`.) From e0e5438167c81abaece9ef208d3197cf7601eb5c Mon Sep 17 00:00:00 2001 From: ananthan-123 <60189117+ananthan-123@users.noreply.github.com> Date: Tue, 18 Feb 2020 17:37:03 +0530 Subject: [PATCH 04/18] Update test_math.py --- Lib/test/test_math.py | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/Lib/test/test_math.py b/Lib/test/test_math.py index b3301f6a5cf74f..58f58a924ca99d 100644 --- a/Lib/test/test_math.py +++ b/Lib/test/test_math.py @@ -974,6 +974,41 @@ def __index__(self): with self.assertRaises(TypeError): math.isqrt(value) + def test_lcm(self): + lcm = math.lcm + self.assertEqual(lcm(0, 0), 0) + self.assertEqual(lcm(1, 0), 0) + self.assertEqual(lcm(-1, 0), 0) + self.assertEqual(lcm(0, 1), 0) + self.assertEqual(lcm(0, -1), 0) + self.assertEqual(lcm(7, 1), 7) + self.assertEqual(lcm(7, -1), 7) + self.assertEqual(lcm(-23, 15), 345) + self.assertEqual(lcm(120, 84), 840) + self.assertEqual(lcm(84, -120), 840) + self.assertEqual(lcm(1216342683557601535506311712, + 436522681849110124616458784), + 16592536571065866494401400422922201534178938447014944) + x = 434610456570399902378880679233098819019853229470286994367836600566 + y = 1064502245825115327754847244914921553977 + + for c in (652560, + 576559230871654959816130551884856912003141446781646602790216406874): + a = x * c + b = y * c + d = x * y * c + self.assertEqual(lcm(a, b), d) + self.assertEqual(lcm(b, a), d) + self.assertEqual(lcm(-a, b), d) + self.assertEqual(lcm(b, -a), d) + self.assertEqual(lcm(a, -b), d) + self.assertEqual(lcm(-b, a), d) + self.assertEqual(lcm(-a, -b), d) + self.assertEqual(lcm(-b, -a), d) + self.assertEqual(lcm(MyIndexable(120), MyIndexable(84)), 840) + self.assertRaises(TypeError, lcm, 120.0, 84) + self.assertRaises(TypeError, lcm, 120, 84.0) + def testLdexp(self): self.assertRaises(TypeError, math.ldexp) self.ftest('ldexp(0,1)', math.ldexp(0,1), 0) From f3aa411bd12e7e42c6afd6e6a88c9cfc7a52ee4d Mon Sep 17 00:00:00 2001 From: ananthan-123 <60189117+ananthan-123@users.noreply.github.com> Date: Tue, 18 Feb 2020 17:48:36 +0530 Subject: [PATCH 05/18] Update mathmodule.c --- Modules/mathmodule.c | 54 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 54 insertions(+) diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c index 309f2291595401..379e16643743ce 100644 --- a/Modules/mathmodule.c +++ b/Modules/mathmodule.c @@ -2017,6 +2017,59 @@ math_factorial(PyObject *module, PyObject *arg) } +/*[clinic input] +math.lcm + x as a: object + y as b: object + / +least common multiple of x and y +[clinic start generated code]*/ + +static PyObject * +math_lcm_impl(PyObject *module, PyObject *a, PyObject *b) +/*[clinic end generated code: output=6f83fb6d671074ba input=bd41b785dc2a4ff1]*/ +{ + PyObject *g, *m, *f, *ab; + + a = PyNumber_Index(a); + if (a == NULL) { + return NULL; + } + b = PyNumber_Index(b); + if (b == NULL) { + Py_DECREF(a); + return NULL; + } + if (_PyLong_Sign(a) == 0 || _PyLong_Sign(b) == 0) { + Py_DECREF(a); + Py_DECREF(b); + return PyLong_FromLong(0); + } + g = _PyLong_GCD(a, b); + if (g == NULL) { + Py_DECREF(a); + Py_DECREF(b); + return NULL; + } + f = PyNumber_FloorDivide(a, g); + Py_DECREF(g); + Py_DECREF(a); + if (f == NULL) { + Py_DECREF(b); + return NULL; + } + m = PyNumber_Multiply(f, b); + Py_DECREF(f); + Py_DECREF(b); + if (m == NULL) { + return NULL; + } + ab = PyNumber_Absolute(m); + Py_DECREF(m); + return ab; +} + + /*[clinic input] math.trunc @@ -3362,6 +3415,7 @@ static PyMethodDef math_methods[] = { MATH_ISINF_METHODDEF MATH_ISNAN_METHODDEF MATH_ISQRT_METHODDEF + MATH_LCM_METHODDEF MATH_LDEXP_METHODDEF {"lgamma", math_lgamma, METH_O, math_lgamma_doc}, MATH_LOG_METHODDEF From ff27d2fb7022525e08ecd7621dd0281bf267a4b8 Mon Sep 17 00:00:00 2001 From: ananthan-123 <60189117+ananthan-123@users.noreply.github.com> Date: Tue, 18 Feb 2020 17:50:52 +0530 Subject: [PATCH 06/18] Update mathmodule.c.h --- Modules/clinic/mathmodule.c.h | 32 +++++++++++++++++++++++++++++++- 1 file changed, 31 insertions(+), 1 deletion(-) diff --git a/Modules/clinic/mathmodule.c.h b/Modules/clinic/mathmodule.c.h index f95d291d41f804..df45a1a0c5e855 100644 --- a/Modules/clinic/mathmodule.c.h +++ b/Modules/clinic/mathmodule.c.h @@ -85,6 +85,36 @@ PyDoc_STRVAR(math_factorial__doc__, #define MATH_FACTORIAL_METHODDEF \ {"factorial", (PyCFunction)math_factorial, METH_O, math_factorial__doc__}, +PyDoc_STRVAR(math_lcm__doc__, +"lcm($module, x, y, /)\n" +"--\n" +"\n" +"least common multiple of x and y"); + +#define MATH_LCM_METHODDEF \ + {"lcm", (PyCFunction)(void(*)(void))math_lcm, METH_FASTCALL, math_lcm__doc__}, + +static PyObject * +math_lcm_impl(PyObject *module, PyObject *a, PyObject *b); + +static PyObject * +math_lcm(PyObject *module, PyObject *const *args, Py_ssize_t nargs) +{ + PyObject *return_value = NULL; + PyObject *a; + PyObject *b; + + if (!_PyArg_CheckPositional("lcm", nargs, 2, 2)) { + goto exit; + } + a = args[0]; + b = args[1]; + return_value = math_lcm_impl(module, a, b); + +exit: + return return_value; +} + PyDoc_STRVAR(math_trunc__doc__, "trunc($module, x, /)\n" "--\n" @@ -895,4 +925,4 @@ math_ulp(PyObject *module, PyObject *arg) exit: return return_value; } -/*[clinic end generated code: output=9b51d215dbcac060 input=a9049054013a1b77]*/ +/*[clinic end generated code: output=f8daa185c043a7b7 input=a9049054013a1b77]*/ From b13f01fd2a9809a73de6d84ac6dd502b853b8c81 Mon Sep 17 00:00:00 2001 From: ananthan-123 <60189117+ananthan-123@users.noreply.github.com> Date: Tue, 18 Feb 2020 17:56:52 +0530 Subject: [PATCH 07/18] Update ACKS --- Misc/ACKS | 1 + 1 file changed, 1 insertion(+) diff --git a/Misc/ACKS b/Misc/ACKS index e8ce30310bfd8e..1c7a73dfe72cc9 100644 --- a/Misc/ACKS +++ b/Misc/ACKS @@ -12,6 +12,7 @@ PS: In the standard Python distribution, this file is encoded in UTF-8 and the list is in rough alphabetical order by last names. Aahz +Ananthakrishnan Edison Abahurire Michael Abbott Rajiv Abraham From 508309db3f6e944952da39f14bcd96068d11c5cb Mon Sep 17 00:00:00 2001 From: "blurb-it[bot]" <43283697+blurb-it[bot]@users.noreply.github.com> Date: Tue, 18 Feb 2020 12:37:17 +0000 Subject: [PATCH 08/18] =?UTF-8?q?=F0=9F=93=9C=F0=9F=A4=96=20Added=20by=20b?= =?UTF-8?q?lurb=5Fit.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../NEWS.d/next/Library/2020-02-18-12-37-16.bpo-39479.j3UcCq.rst | 1 + 1 file changed, 1 insertion(+) create mode 100644 Misc/NEWS.d/next/Library/2020-02-18-12-37-16.bpo-39479.j3UcCq.rst diff --git a/Misc/NEWS.d/next/Library/2020-02-18-12-37-16.bpo-39479.j3UcCq.rst b/Misc/NEWS.d/next/Library/2020-02-18-12-37-16.bpo-39479.j3UcCq.rst new file mode 100644 index 00000000000000..7bbdcb3aa1cb7e --- /dev/null +++ b/Misc/NEWS.d/next/Library/2020-02-18-12-37-16.bpo-39479.j3UcCq.rst @@ -0,0 +1 @@ +Add func:math.lcm function: Least Common Multiple. \ No newline at end of file From 49f1b4b8d4fd38568d068dc1dd427b3535fbfd57 Mon Sep 17 00:00:00 2001 From: ananthan-123 <60189117+ananthan-123@users.noreply.github.com> Date: Tue, 18 Feb 2020 18:54:35 +0530 Subject: [PATCH 09/18] Update 3.9.rst --- Doc/whatsnew/3.9.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Doc/whatsnew/3.9.rst b/Doc/whatsnew/3.9.rst index 7e495563e25ae8..c8f650972d298f 100644 --- a/Doc/whatsnew/3.9.rst +++ b/Doc/whatsnew/3.9.rst @@ -216,7 +216,7 @@ import attempts. math ---- -Add :func:`math.lcm`: return the least common multiple of a and b. +Add :func:`math.lcm`: return the least common multiple of *a* and *b*. (Contributed by Ananthakrishnan in :issue:`39479`.) Add :func:`math.nextafter`: return the next floating-point value after *x* From cd06d73a65768427061c0a4a059f031b8eb1b21d Mon Sep 17 00:00:00 2001 From: ananthan-123 <60189117+ananthan-123@users.noreply.github.com> Date: Tue, 18 Feb 2020 19:56:49 +0530 Subject: [PATCH 10/18] Update 2020-02-18-12-37-16.bpo-39479.j3UcCq.rst --- .../next/Library/2020-02-18-12-37-16.bpo-39479.j3UcCq.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Misc/NEWS.d/next/Library/2020-02-18-12-37-16.bpo-39479.j3UcCq.rst b/Misc/NEWS.d/next/Library/2020-02-18-12-37-16.bpo-39479.j3UcCq.rst index 7bbdcb3aa1cb7e..1620adeff7c361 100644 --- a/Misc/NEWS.d/next/Library/2020-02-18-12-37-16.bpo-39479.j3UcCq.rst +++ b/Misc/NEWS.d/next/Library/2020-02-18-12-37-16.bpo-39479.j3UcCq.rst @@ -1 +1 @@ -Add func:math.lcm function: Least Common Multiple. \ No newline at end of file +Add :func:`math.lcm` function: Least Common Multiple. From 62fd3e7211f6d2c98680ba0267d37bef2760fb2a Mon Sep 17 00:00:00 2001 From: ananthan-123 Date: Wed, 19 Feb 2020 00:30:21 +0530 Subject: [PATCH 11/18] Update math.rst --- Doc/library/math.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Doc/library/math.rst b/Doc/library/math.rst index b89af90a3f82f7..d8ac35256d1b4b 100644 --- a/Doc/library/math.rst +++ b/Doc/library/math.rst @@ -142,7 +142,7 @@ Number-theoretic and representation functions ``lcm(a, b)`` is the smallest nonnegative integer that is a multiple of both *a* and *b*. If either *a* or *b* is zero then ``lcm(a, b)`` is zero. - .. versionadded:: 3.5 + .. versionadded:: 3.9 .. function:: isclose(a, b, *, rel_tol=1e-09, abs_tol=0.0) From 5f2f7f693b32f5ef939c778b3ebe013aea54d626 Mon Sep 17 00:00:00 2001 From: ananthan-123 Date: Wed, 19 Feb 2020 08:22:02 +0530 Subject: [PATCH 12/18] Update 2020-02-18-12-37-16.bpo-39479.j3UcCq.rst --- .../next/Library/2020-02-18-12-37-16.bpo-39479.j3UcCq.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Misc/NEWS.d/next/Library/2020-02-18-12-37-16.bpo-39479.j3UcCq.rst b/Misc/NEWS.d/next/Library/2020-02-18-12-37-16.bpo-39479.j3UcCq.rst index 1620adeff7c361..6f16623a8f5cd1 100644 --- a/Misc/NEWS.d/next/Library/2020-02-18-12-37-16.bpo-39479.j3UcCq.rst +++ b/Misc/NEWS.d/next/Library/2020-02-18-12-37-16.bpo-39479.j3UcCq.rst @@ -1 +1 @@ -Add :func:`math.lcm` function: Least Common Multiple. +Add :func:`math.lcm` function: least common multiple. From d3f527842d21105d262235a652f6c42fc3c48ce9 Mon Sep 17 00:00:00 2001 From: ananthan-123 Date: Wed, 19 Feb 2020 15:22:58 +0530 Subject: [PATCH 13/18] Update test_math.py --- Lib/test/test_math.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Lib/test/test_math.py b/Lib/test/test_math.py index 58f58a924ca99d..ad8273d50cc39a 100644 --- a/Lib/test/test_math.py +++ b/Lib/test/test_math.py @@ -989,11 +989,11 @@ def test_lcm(self): self.assertEqual(lcm(1216342683557601535506311712, 436522681849110124616458784), 16592536571065866494401400422922201534178938447014944) - x = 434610456570399902378880679233098819019853229470286994367836600566 - y = 1064502245825115327754847244914921553977 + x = 43461045657039990237 + y = 10645022458251153277 for c in (652560, - 576559230871654959816130551884856912003141446781646602790216406874): + 57655923087165495981): a = x * c b = y * c d = x * y * c From 5e902a148231da156219b61a6f3b0a437a67472b Mon Sep 17 00:00:00 2001 From: ananthan-123 Date: Wed, 19 Feb 2020 19:21:12 +0530 Subject: [PATCH 14/18] Update ACKS --- Misc/ACKS | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Misc/ACKS b/Misc/ACKS index 1c7a73dfe72cc9..b11ee9cdea67d7 100644 --- a/Misc/ACKS +++ b/Misc/ACKS @@ -12,7 +12,6 @@ PS: In the standard Python distribution, this file is encoded in UTF-8 and the list is in rough alphabetical order by last names. Aahz -Ananthakrishnan Edison Abahurire Michael Abbott Rajiv Abraham @@ -46,6 +45,7 @@ Rose Ames A. Amoroso Mark Anacker Shashwat Anand +Ananthakrishnan Anders Andersen Tycho Andersen John Anderson From cb0e634aeabde0a2d4154a98f485eb214211316b Mon Sep 17 00:00:00 2001 From: ananthan-123 Date: Wed, 19 Feb 2020 20:45:57 +0530 Subject: [PATCH 15/18] Update mathmodule.c.h --- Modules/clinic/mathmodule.c.h | 4051 +++++++++++++++++++++++++++------ 1 file changed, 3301 insertions(+), 750 deletions(-) diff --git a/Modules/clinic/mathmodule.c.h b/Modules/clinic/mathmodule.c.h index df45a1a0c5e855..93158246d389b5 100644 --- a/Modules/clinic/mathmodule.c.h +++ b/Modules/clinic/mathmodule.c.h @@ -1,928 +1,3479 @@ -/*[clinic input] -preserve -[clinic start generated code]*/ - -PyDoc_STRVAR(math_gcd__doc__, -"gcd($module, x, y, /)\n" -"--\n" -"\n" -"greatest common divisor of x and y"); - -#define MATH_GCD_METHODDEF \ - {"gcd", (PyCFunction)(void(*)(void))math_gcd, METH_FASTCALL, math_gcd__doc__}, - -static PyObject * -math_gcd_impl(PyObject *module, PyObject *a, PyObject *b); - -static PyObject * -math_gcd(PyObject *module, PyObject *const *args, Py_ssize_t nargs) -{ - PyObject *return_value = NULL; - PyObject *a; - PyObject *b; +/* Math module -- standard C math library functions, pi and e */ - if (!_PyArg_CheckPositional("gcd", nargs, 2, 2)) { - goto exit; - } - a = args[0]; - b = args[1]; - return_value = math_gcd_impl(module, a, b); +/* Here are some comments from Tim Peters, extracted from the + discussion attached to http://bugs.python.org/issue1640. They + describe the general aims of the math module with respect to + special values, IEEE-754 floating-point exceptions, and Python + exceptions. -exit: - return return_value; -} +These are the "spirit of 754" rules: -PyDoc_STRVAR(math_ceil__doc__, -"ceil($module, x, /)\n" -"--\n" -"\n" -"Return the ceiling of x as an Integral.\n" -"\n" -"This is the smallest integer >= x."); +1. If the mathematical result is a real number, but of magnitude too +large to approximate by a machine float, overflow is signaled and the +result is an infinity (with the appropriate sign). -#define MATH_CEIL_METHODDEF \ - {"ceil", (PyCFunction)math_ceil, METH_O, math_ceil__doc__}, +2. If the mathematical result is a real number, but of magnitude too +small to approximate by a machine float, underflow is signaled and the +result is a zero (with the appropriate sign). -PyDoc_STRVAR(math_floor__doc__, -"floor($module, x, /)\n" -"--\n" -"\n" -"Return the floor of x as an Integral.\n" -"\n" -"This is the largest integer <= x."); +3. At a singularity (a value x such that the limit of f(y) as y +approaches x exists and is an infinity), "divide by zero" is signaled +and the result is an infinity (with the appropriate sign). This is +complicated a little by that the left-side and right-side limits may +not be the same; e.g., 1/x approaches +inf or -inf as x approaches 0 +from the positive or negative directions. In that specific case, the +sign of the zero determines the result of 1/0. -#define MATH_FLOOR_METHODDEF \ - {"floor", (PyCFunction)math_floor, METH_O, math_floor__doc__}, +4. At a point where a function has no defined result in the extended +reals (i.e., the reals plus an infinity or two), invalid operation is +signaled and a NaN is returned. -PyDoc_STRVAR(math_fsum__doc__, -"fsum($module, seq, /)\n" -"--\n" -"\n" -"Return an accurate floating point sum of values in the iterable seq.\n" -"\n" -"Assumes IEEE-754 floating point arithmetic."); +And these are what Python has historically /tried/ to do (but not +always successfully, as platform libm behavior varies a lot): -#define MATH_FSUM_METHODDEF \ - {"fsum", (PyCFunction)math_fsum, METH_O, math_fsum__doc__}, +For #1, raise OverflowError. -PyDoc_STRVAR(math_isqrt__doc__, -"isqrt($module, n, /)\n" -"--\n" -"\n" -"Return the integer part of the square root of the input."); +For #2, return a zero (with the appropriate sign if that happens by +accident ;-)). -#define MATH_ISQRT_METHODDEF \ - {"isqrt", (PyCFunction)math_isqrt, METH_O, math_isqrt__doc__}, +For #3 and #4, raise ValueError. It may have made sense to raise +Python's ZeroDivisionError in #3, but historically that's only been +raised for division by zero and mod by zero. -PyDoc_STRVAR(math_factorial__doc__, -"factorial($module, x, /)\n" -"--\n" -"\n" -"Find x!.\n" -"\n" -"Raise a ValueError if x is negative or non-integral."); +*/ -#define MATH_FACTORIAL_METHODDEF \ - {"factorial", (PyCFunction)math_factorial, METH_O, math_factorial__doc__}, +/* + In general, on an IEEE-754 platform the aim is to follow the C99 + standard, including Annex 'F', whenever possible. Where the + standard recommends raising the 'divide-by-zero' or 'invalid' + floating-point exceptions, Python should raise a ValueError. Where + the standard recommends raising 'overflow', Python should raise an + OverflowError. In all other circumstances a value should be + returned. + */ -PyDoc_STRVAR(math_lcm__doc__, -"lcm($module, x, y, /)\n" -"--\n" -"\n" -"least common multiple of x and y"); +#include "Python.h" +#include "pycore_dtoa.h" +#include "_math.h" -#define MATH_LCM_METHODDEF \ - {"lcm", (PyCFunction)(void(*)(void))math_lcm, METH_FASTCALL, math_lcm__doc__}, +#include "clinic/mathmodule.c.h" -static PyObject * -math_lcm_impl(PyObject *module, PyObject *a, PyObject *b); +/*[clinic input] +module math +[clinic start generated code]*/ +/*[clinic end generated code: output=da39a3ee5e6b4b0d input=76bc7002685dd942]*/ + + +/* + sin(pi*x), giving accurate results for all finite x (especially x + integral or close to an integer). This is here for use in the + reflection formula for the gamma function. It conforms to IEEE + 754-2008 for finite arguments, but not for infinities or nans. +*/ + +static const double pi = 3.141592653589793238462643383279502884197; +static const double logpi = 1.144729885849400174143427351353058711647; +#if !defined(HAVE_ERF) || !defined(HAVE_ERFC) +static const double sqrtpi = 1.772453850905516027298167483341145182798; +#endif /* !defined(HAVE_ERF) || !defined(HAVE_ERFC) */ + + +/* Version of PyFloat_AsDouble() with in-line fast paths + for exact floats and integers. Gives a substantial + speed improvement for extracting float arguments. +*/ + +#define ASSIGN_DOUBLE(target_var, obj, error_label) \ + if (PyFloat_CheckExact(obj)) { \ + target_var = PyFloat_AS_DOUBLE(obj); \ + } \ + else if (PyLong_CheckExact(obj)) { \ + target_var = PyLong_AsDouble(obj); \ + if (target_var == -1.0 && PyErr_Occurred()) { \ + goto error_label; \ + } \ + } \ + else { \ + target_var = PyFloat_AsDouble(obj); \ + if (target_var == -1.0 && PyErr_Occurred()) { \ + goto error_label; \ + } \ + } -static PyObject * -math_lcm(PyObject *module, PyObject *const *args, Py_ssize_t nargs) +static double +m_sinpi(double x) { - PyObject *return_value = NULL; - PyObject *a; - PyObject *b; - - if (!_PyArg_CheckPositional("lcm", nargs, 2, 2)) { - goto exit; + double y, r; + int n; + /* this function should only ever be called for finite arguments */ + assert(Py_IS_FINITE(x)); + y = fmod(fabs(x), 2.0); + n = (int)round(2.0*y); + assert(0 <= n && n <= 4); + switch (n) { + case 0: + r = sin(pi*y); + break; + case 1: + r = cos(pi*(y-0.5)); + break; + case 2: + /* N.B. -sin(pi*(y-1.0)) is *not* equivalent: it would give + -0.0 instead of 0.0 when y == 1.0. */ + r = sin(pi*(1.0-y)); + break; + case 3: + r = -cos(pi*(y-1.5)); + break; + case 4: + r = sin(pi*(y-2.0)); + break; + default: + Py_UNREACHABLE(); } - a = args[0]; - b = args[1]; - return_value = math_lcm_impl(module, a, b); - -exit: - return return_value; + return copysign(1.0, x)*r; } -PyDoc_STRVAR(math_trunc__doc__, -"trunc($module, x, /)\n" -"--\n" -"\n" -"Truncates the Real x to the nearest Integral toward 0.\n" -"\n" -"Uses the __trunc__ magic method."); - -#define MATH_TRUNC_METHODDEF \ - {"trunc", (PyCFunction)math_trunc, METH_O, math_trunc__doc__}, - -PyDoc_STRVAR(math_frexp__doc__, -"frexp($module, x, /)\n" -"--\n" -"\n" -"Return the mantissa and exponent of x, as pair (m, e).\n" -"\n" -"m is a float and e is an int, such that x = m * 2.**e.\n" -"If x is 0, m and e are both 0. Else 0.5 <= abs(m) < 1.0."); - -#define MATH_FREXP_METHODDEF \ - {"frexp", (PyCFunction)math_frexp, METH_O, math_frexp__doc__}, - -static PyObject * -math_frexp_impl(PyObject *module, double x); +/* Implementation of the real gamma function. In extensive but non-exhaustive + random tests, this function proved accurate to within <= 10 ulps across the + entire float domain. Note that accuracy may depend on the quality of the + system math functions, the pow function in particular. Special cases + follow C99 annex F. The parameters and method are tailored to platforms + whose double format is the IEEE 754 binary64 format. + + Method: for x > 0.0 we use the Lanczos approximation with parameters N=13 + and g=6.024680040776729583740234375; these parameters are amongst those + used by the Boost library. Following Boost (again), we re-express the + Lanczos sum as a rational function, and compute it that way. The + coefficients below were computed independently using MPFR, and have been + double-checked against the coefficients in the Boost source code. + + For x < 0.0 we use the reflection formula. + + There's one minor tweak that deserves explanation: Lanczos' formula for + Gamma(x) involves computing pow(x+g-0.5, x-0.5) / exp(x+g-0.5). For many x + values, x+g-0.5 can be represented exactly. However, in cases where it + can't be represented exactly the small error in x+g-0.5 can be magnified + significantly by the pow and exp calls, especially for large x. A cheap + correction is to multiply by (1 + e*g/(x+g-0.5)), where e is the error + involved in the computation of x+g-0.5 (that is, e = computed value of + x+g-0.5 - exact value of x+g-0.5). Here's the proof: + + Correction factor + ----------------- + Write x+g-0.5 = y-e, where y is exactly representable as an IEEE 754 + double, and e is tiny. Then: + + pow(x+g-0.5,x-0.5)/exp(x+g-0.5) = pow(y-e, x-0.5)/exp(y-e) + = pow(y, x-0.5)/exp(y) * C, + + where the correction_factor C is given by + + C = pow(1-e/y, x-0.5) * exp(e) + + Since e is tiny, pow(1-e/y, x-0.5) ~ 1-(x-0.5)*e/y, and exp(x) ~ 1+e, so: + + C ~ (1-(x-0.5)*e/y) * (1+e) ~ 1 + e*(y-(x-0.5))/y + + But y-(x-0.5) = g+e, and g+e ~ g. So we get C ~ 1 + e*g/y, and + + pow(x+g-0.5,x-0.5)/exp(x+g-0.5) ~ pow(y, x-0.5)/exp(y) * (1 + e*g/y), + + Note that for accuracy, when computing r*C it's better to do + + r + e*g/y*r; + + than + + r * (1 + e*g/y); + + since the addition in the latter throws away most of the bits of + information in e*g/y. +*/ + +#define LANCZOS_N 13 +static const double lanczos_g = 6.024680040776729583740234375; +static const double lanczos_g_minus_half = 5.524680040776729583740234375; +static const double lanczos_num_coeffs[LANCZOS_N] = { + 23531376880.410759688572007674451636754734846804940, + 42919803642.649098768957899047001988850926355848959, + 35711959237.355668049440185451547166705960488635843, + 17921034426.037209699919755754458931112671403265390, + 6039542586.3520280050642916443072979210699388420708, + 1439720407.3117216736632230727949123939715485786772, + 248874557.86205415651146038641322942321632125127801, + 31426415.585400194380614231628318205362874684987640, + 2876370.6289353724412254090516208496135991145378768, + 186056.26539522349504029498971604569928220784236328, + 8071.6720023658162106380029022722506138218516325024, + 210.82427775157934587250973392071336271166969580291, + 2.5066282746310002701649081771338373386264310793408 +}; + +/* denominator is x*(x+1)*...*(x+LANCZOS_N-2) */ +static const double lanczos_den_coeffs[LANCZOS_N] = { + 0.0, 39916800.0, 120543840.0, 150917976.0, 105258076.0, 45995730.0, + 13339535.0, 2637558.0, 357423.0, 32670.0, 1925.0, 66.0, 1.0}; + +/* gamma values for small positive integers, 1 though NGAMMA_INTEGRAL */ +#define NGAMMA_INTEGRAL 23 +static const double gamma_integral[NGAMMA_INTEGRAL] = { + 1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0, 362880.0, + 3628800.0, 39916800.0, 479001600.0, 6227020800.0, 87178291200.0, + 1307674368000.0, 20922789888000.0, 355687428096000.0, + 6402373705728000.0, 121645100408832000.0, 2432902008176640000.0, + 51090942171709440000.0, 1124000727777607680000.0, +}; + +/* Lanczos' sum L_g(x), for positive x */ -static PyObject * -math_frexp(PyObject *module, PyObject *arg) +static double +lanczos_sum(double x) { - PyObject *return_value = NULL; - double x; - - if (PyFloat_CheckExact(arg)) { - x = PyFloat_AS_DOUBLE(arg); + double num = 0.0, den = 0.0; + int i; + assert(x > 0.0); + /* evaluate the rational function lanczos_sum(x). For large + x, the obvious algorithm risks overflow, so we instead + rescale the denominator and numerator of the rational + function by x**(1-LANCZOS_N) and treat this as a + rational function in 1/x. This also reduces the error for + larger x values. The choice of cutoff point (5.0 below) is + somewhat arbitrary; in tests, smaller cutoff values than + this resulted in lower accuracy. */ + if (x < 5.0) { + for (i = LANCZOS_N; --i >= 0; ) { + num = num * x + lanczos_num_coeffs[i]; + den = den * x + lanczos_den_coeffs[i]; + } } - else - { - x = PyFloat_AsDouble(arg); - if (x == -1.0 && PyErr_Occurred()) { - goto exit; + else { + for (i = 0; i < LANCZOS_N; i++) { + num = num / x + lanczos_num_coeffs[i]; + den = den / x + lanczos_den_coeffs[i]; } } - return_value = math_frexp_impl(module, x); + return num/den; +} + +/* Constant for +infinity, generated in the same way as float('inf'). */ -exit: - return return_value; +static double +m_inf(void) +{ +#ifndef PY_NO_SHORT_FLOAT_REPR + return _Py_dg_infinity(0); +#else + return Py_HUGE_VAL; +#endif } -PyDoc_STRVAR(math_ldexp__doc__, -"ldexp($module, x, i, /)\n" -"--\n" -"\n" -"Return x * (2**i).\n" -"\n" -"This is essentially the inverse of frexp()."); +/* Constant nan value, generated in the same way as float('nan'). */ +/* We don't currently assume that Py_NAN is defined everywhere. */ -#define MATH_LDEXP_METHODDEF \ - {"ldexp", (PyCFunction)(void(*)(void))math_ldexp, METH_FASTCALL, math_ldexp__doc__}, +#if !defined(PY_NO_SHORT_FLOAT_REPR) || defined(Py_NAN) -static PyObject * -math_ldexp_impl(PyObject *module, double x, PyObject *i); +static double +m_nan(void) +{ +#ifndef PY_NO_SHORT_FLOAT_REPR + return _Py_dg_stdnan(0); +#else + return Py_NAN; +#endif +} -static PyObject * -math_ldexp(PyObject *module, PyObject *const *args, Py_ssize_t nargs) +#endif + +static double +m_tgamma(double x) { - PyObject *return_value = NULL; - double x; - PyObject *i; + double absx, r, y, z, sqrtpow; + + /* special cases */ + if (!Py_IS_FINITE(x)) { + if (Py_IS_NAN(x) || x > 0.0) + return x; /* tgamma(nan) = nan, tgamma(inf) = inf */ + else { + errno = EDOM; + return Py_NAN; /* tgamma(-inf) = nan, invalid */ + } + } + if (x == 0.0) { + errno = EDOM; + /* tgamma(+-0.0) = +-inf, divide-by-zero */ + return copysign(Py_HUGE_VAL, x); + } - if (!_PyArg_CheckPositional("ldexp", nargs, 2, 2)) { - goto exit; + /* integer arguments */ + if (x == floor(x)) { + if (x < 0.0) { + errno = EDOM; /* tgamma(n) = nan, invalid for */ + return Py_NAN; /* negative integers n */ + } + if (x <= NGAMMA_INTEGRAL) + return gamma_integral[(int)x - 1]; } - if (PyFloat_CheckExact(args[0])) { - x = PyFloat_AS_DOUBLE(args[0]); + absx = fabs(x); + + /* tiny arguments: tgamma(x) ~ 1/x for x near 0 */ + if (absx < 1e-20) { + r = 1.0/x; + if (Py_IS_INFINITY(r)) + errno = ERANGE; + return r; } - else - { - x = PyFloat_AsDouble(args[0]); - if (x == -1.0 && PyErr_Occurred()) { - goto exit; + + /* large arguments: assuming IEEE 754 doubles, tgamma(x) overflows for + x > 200, and underflows to +-0.0 for x < -200, not a negative + integer. */ + if (absx > 200.0) { + if (x < 0.0) { + return 0.0/m_sinpi(x); + } + else { + errno = ERANGE; + return Py_HUGE_VAL; } } - i = args[1]; - return_value = math_ldexp_impl(module, x, i); -exit: - return return_value; + y = absx + lanczos_g_minus_half; + /* compute error in sum */ + if (absx > lanczos_g_minus_half) { + /* note: the correction can be foiled by an optimizing + compiler that (incorrectly) thinks that an expression like + a + b - a - b can be optimized to 0.0. This shouldn't + happen in a standards-conforming compiler. */ + double q = y - absx; + z = q - lanczos_g_minus_half; + } + else { + double q = y - lanczos_g_minus_half; + z = q - absx; + } + z = z * lanczos_g / y; + if (x < 0.0) { + r = -pi / m_sinpi(absx) / absx * exp(y) / lanczos_sum(absx); + r -= z * r; + if (absx < 140.0) { + r /= pow(y, absx - 0.5); + } + else { + sqrtpow = pow(y, absx / 2.0 - 0.25); + r /= sqrtpow; + r /= sqrtpow; + } + } + else { + r = lanczos_sum(absx) / exp(y); + r += z * r; + if (absx < 140.0) { + r *= pow(y, absx - 0.5); + } + else { + sqrtpow = pow(y, absx / 2.0 - 0.25); + r *= sqrtpow; + r *= sqrtpow; + } + } + if (Py_IS_INFINITY(r)) + errno = ERANGE; + return r; } -PyDoc_STRVAR(math_modf__doc__, -"modf($module, x, /)\n" -"--\n" -"\n" -"Return the fractional and integer parts of x.\n" -"\n" -"Both results carry the sign of x and are floats."); - -#define MATH_MODF_METHODDEF \ - {"modf", (PyCFunction)math_modf, METH_O, math_modf__doc__}, - -static PyObject * -math_modf_impl(PyObject *module, double x); +/* + lgamma: natural log of the absolute value of the Gamma function. + For large arguments, Lanczos' formula works extremely well here. +*/ -static PyObject * -math_modf(PyObject *module, PyObject *arg) +static double +m_lgamma(double x) { - PyObject *return_value = NULL; - double x; + double r; + double absx; - if (PyFloat_CheckExact(arg)) { - x = PyFloat_AS_DOUBLE(arg); + /* special cases */ + if (!Py_IS_FINITE(x)) { + if (Py_IS_NAN(x)) + return x; /* lgamma(nan) = nan */ + else + return Py_HUGE_VAL; /* lgamma(+-inf) = +inf */ } - else - { - x = PyFloat_AsDouble(arg); - if (x == -1.0 && PyErr_Occurred()) { - goto exit; + + /* integer arguments */ + if (x == floor(x) && x <= 2.0) { + if (x <= 0.0) { + errno = EDOM; /* lgamma(n) = inf, divide-by-zero for */ + return Py_HUGE_VAL; /* integers n <= 0 */ + } + else { + return 0.0; /* lgamma(1) = lgamma(2) = 0.0 */ } } - return_value = math_modf_impl(module, x); -exit: - return return_value; + absx = fabs(x); + /* tiny arguments: lgamma(x) ~ -log(fabs(x)) for small x */ + if (absx < 1e-20) + return -log(absx); + + /* Lanczos' formula. We could save a fraction of a ulp in accuracy by + having a second set of numerator coefficients for lanczos_sum that + absorbed the exp(-lanczos_g) term, and throwing out the lanczos_g + subtraction below; it's probably not worth it. */ + r = log(lanczos_sum(absx)) - lanczos_g; + r += (absx - 0.5) * (log(absx + lanczos_g - 0.5) - 1); + if (x < 0.0) + /* Use reflection formula to get value for negative x. */ + r = logpi - log(fabs(m_sinpi(absx))) - log(absx) - r; + if (Py_IS_INFINITY(r)) + errno = ERANGE; + return r; } -PyDoc_STRVAR(math_log__doc__, -"log(x, [base=math.e])\n" -"Return the logarithm of x to the given base.\n" -"\n" -"If the base not specified, returns the natural logarithm (base e) of x."); +#if !defined(HAVE_ERF) || !defined(HAVE_ERFC) -#define MATH_LOG_METHODDEF \ - {"log", (PyCFunction)math_log, METH_VARARGS, math_log__doc__}, +/* + Implementations of the error function erf(x) and the complementary error + function erfc(x). -static PyObject * -math_log_impl(PyObject *module, PyObject *x, int group_right_1, - PyObject *base); + Method: we use a series approximation for erf for small x, and a continued + fraction approximation for erfc(x) for larger x; + combined with the relations erf(-x) = -erf(x) and erfc(x) = 1.0 - erf(x), + this gives us erf(x) and erfc(x) for all x. -static PyObject * -math_log(PyObject *module, PyObject *args) -{ - PyObject *return_value = NULL; - PyObject *x; - int group_right_1 = 0; - PyObject *base = NULL; + The series expansion used is: - switch (PyTuple_GET_SIZE(args)) { - case 1: - if (!PyArg_ParseTuple(args, "O:log", &x)) { - goto exit; - } - break; - case 2: - if (!PyArg_ParseTuple(args, "OO:log", &x, &base)) { - goto exit; - } - group_right_1 = 1; - break; - default: - PyErr_SetString(PyExc_TypeError, "math.log requires 1 to 2 arguments"); - goto exit; - } - return_value = math_log_impl(module, x, group_right_1, base); + erf(x) = x*exp(-x*x)/sqrt(pi) * [ + 2/1 + 4/3 x**2 + 8/15 x**4 + 16/105 x**6 + ...] -exit: - return return_value; -} + The coefficient of x**(2k-2) here is 4**k*factorial(k)/factorial(2*k). + This series converges well for smallish x, but slowly for larger x. -PyDoc_STRVAR(math_log2__doc__, -"log2($module, x, /)\n" -"--\n" -"\n" -"Return the base 2 logarithm of x."); + The continued fraction expansion used is: -#define MATH_LOG2_METHODDEF \ - {"log2", (PyCFunction)math_log2, METH_O, math_log2__doc__}, + erfc(x) = x*exp(-x*x)/sqrt(pi) * [1/(0.5 + x**2 -) 0.5/(2.5 + x**2 - ) + 3.0/(4.5 + x**2 - ) 7.5/(6.5 + x**2 - ) ...] -PyDoc_STRVAR(math_log10__doc__, -"log10($module, x, /)\n" -"--\n" -"\n" -"Return the base 10 logarithm of x."); + after the first term, the general term has the form: -#define MATH_LOG10_METHODDEF \ - {"log10", (PyCFunction)math_log10, METH_O, math_log10__doc__}, + k*(k-0.5)/(2*k+0.5 + x**2 - ...). -PyDoc_STRVAR(math_fmod__doc__, -"fmod($module, x, y, /)\n" -"--\n" -"\n" -"Return fmod(x, y), according to platform C.\n" -"\n" -"x % y may differ."); + This expansion converges fast for larger x, but convergence becomes + infinitely slow as x approaches 0.0. The (somewhat naive) continued + fraction evaluation algorithm used below also risks overflow for large x; + but for large x, erfc(x) == 0.0 to within machine precision. (For + example, erfc(30.0) is approximately 2.56e-393). -#define MATH_FMOD_METHODDEF \ - {"fmod", (PyCFunction)(void(*)(void))math_fmod, METH_FASTCALL, math_fmod__doc__}, + Parameters: use series expansion for abs(x) < ERF_SERIES_CUTOFF and + continued fraction expansion for ERF_SERIES_CUTOFF <= abs(x) < + ERFC_CONTFRAC_CUTOFF. ERFC_SERIES_TERMS and ERFC_CONTFRAC_TERMS are the + numbers of terms to use for the relevant expansions. */ -static PyObject * -math_fmod_impl(PyObject *module, double x, double y); +#define ERF_SERIES_CUTOFF 1.5 +#define ERF_SERIES_TERMS 25 +#define ERFC_CONTFRAC_CUTOFF 30.0 +#define ERFC_CONTFRAC_TERMS 50 -static PyObject * -math_fmod(PyObject *module, PyObject *const *args, Py_ssize_t nargs) -{ - PyObject *return_value = NULL; - double x; - double y; +/* + Error function, via power series. - if (!_PyArg_CheckPositional("fmod", nargs, 2, 2)) { - goto exit; - } - if (PyFloat_CheckExact(args[0])) { - x = PyFloat_AS_DOUBLE(args[0]); - } - else - { - x = PyFloat_AsDouble(args[0]); - if (x == -1.0 && PyErr_Occurred()) { - goto exit; - } - } - if (PyFloat_CheckExact(args[1])) { - y = PyFloat_AS_DOUBLE(args[1]); - } - else - { - y = PyFloat_AsDouble(args[1]); - if (y == -1.0 && PyErr_Occurred()) { - goto exit; - } - } - return_value = math_fmod_impl(module, x, y); + Given a finite float x, return an approximation to erf(x). + Converges reasonably fast for small x. +*/ -exit: - return return_value; +static double +m_erf_series(double x) +{ + double x2, acc, fk, result; + int i, saved_errno; + + x2 = x * x; + acc = 0.0; + fk = (double)ERF_SERIES_TERMS + 0.5; + for (i = 0; i < ERF_SERIES_TERMS; i++) { + acc = 2.0 + x2 * acc / fk; + fk -= 1.0; + } + /* Make sure the exp call doesn't affect errno; + see m_erfc_contfrac for more. */ + saved_errno = errno; + result = acc * x * exp(-x2) / sqrtpi; + errno = saved_errno; + return result; } -PyDoc_STRVAR(math_dist__doc__, -"dist($module, p, q, /)\n" -"--\n" -"\n" -"Return the Euclidean distance between two points p and q.\n" -"\n" -"The points should be specified as sequences (or iterables) of\n" -"coordinates. Both inputs must have the same dimension.\n" -"\n" -"Roughly equivalent to:\n" -" sqrt(sum((px - qx) ** 2.0 for px, qx in zip(p, q)))"); +/* + Complementary error function, via continued fraction expansion. -#define MATH_DIST_METHODDEF \ - {"dist", (PyCFunction)(void(*)(void))math_dist, METH_FASTCALL, math_dist__doc__}, - -static PyObject * -math_dist_impl(PyObject *module, PyObject *p, PyObject *q); + Given a positive float x, return an approximation to erfc(x). Converges + reasonably fast for x large (say, x > 2.0), and should be safe from + overflow if x and nterms are not too large. On an IEEE 754 machine, with x + <= 30.0, we're safe up to nterms = 100. For x >= 30.0, erfc(x) is smaller + than the smallest representable nonzero float. */ -static PyObject * -math_dist(PyObject *module, PyObject *const *args, Py_ssize_t nargs) +static double +m_erfc_contfrac(double x) { - PyObject *return_value = NULL; - PyObject *p; - PyObject *q; - - if (!_PyArg_CheckPositional("dist", nargs, 2, 2)) { - goto exit; + double x2, a, da, p, p_last, q, q_last, b, result; + int i, saved_errno; + + if (x >= ERFC_CONTFRAC_CUTOFF) + return 0.0; + + x2 = x*x; + a = 0.0; + da = 0.5; + p = 1.0; p_last = 0.0; + q = da + x2; q_last = 1.0; + for (i = 0; i < ERFC_CONTFRAC_TERMS; i++) { + double temp; + a += da; + da += 2.0; + b = da + x2; + temp = p; p = b*p - a*p_last; p_last = temp; + temp = q; q = b*q - a*q_last; q_last = temp; } - p = args[0]; - q = args[1]; - return_value = math_dist_impl(module, p, q); - -exit: - return return_value; + /* Issue #8986: On some platforms, exp sets errno on underflow to zero; + save the current errno value so that we can restore it later. */ + saved_errno = errno; + result = p / q * x * exp(-x2) / sqrtpi; + errno = saved_errno; + return result; } -PyDoc_STRVAR(math_pow__doc__, -"pow($module, x, y, /)\n" -"--\n" -"\n" -"Return x**y (x to the power of y)."); +#endif /* !defined(HAVE_ERF) || !defined(HAVE_ERFC) */ -#define MATH_POW_METHODDEF \ - {"pow", (PyCFunction)(void(*)(void))math_pow, METH_FASTCALL, math_pow__doc__}, +/* Error function erf(x), for general x */ -static PyObject * -math_pow_impl(PyObject *module, double x, double y); - -static PyObject * -math_pow(PyObject *module, PyObject *const *args, Py_ssize_t nargs) +static double +m_erf(double x) { - PyObject *return_value = NULL; - double x; - double y; - - if (!_PyArg_CheckPositional("pow", nargs, 2, 2)) { - goto exit; +#ifdef HAVE_ERF + return erf(x); +#else + double absx, cf; + + if (Py_IS_NAN(x)) + return x; + absx = fabs(x); + if (absx < ERF_SERIES_CUTOFF) + return m_erf_series(x); + else { + cf = m_erfc_contfrac(absx); + return x > 0.0 ? 1.0 - cf : cf - 1.0; } - if (PyFloat_CheckExact(args[0])) { - x = PyFloat_AS_DOUBLE(args[0]); +#endif +} + +/* Complementary error function erfc(x), for general x. */ + +static double +m_erfc(double x) +{ +#ifdef HAVE_ERFC + return erfc(x); +#else + double absx, cf; + + if (Py_IS_NAN(x)) + return x; + absx = fabs(x); + if (absx < ERF_SERIES_CUTOFF) + return 1.0 - m_erf_series(x); + else { + cf = m_erfc_contfrac(absx); + return x > 0.0 ? cf : 2.0 - cf; } - else - { - x = PyFloat_AsDouble(args[0]); - if (x == -1.0 && PyErr_Occurred()) { - goto exit; +#endif +} + +/* + wrapper for atan2 that deals directly with special cases before + delegating to the platform libm for the remaining cases. This + is necessary to get consistent behaviour across platforms. + Windows, FreeBSD and alpha Tru64 are amongst platforms that don't + always follow C99. +*/ + +static double +m_atan2(double y, double x) +{ + if (Py_IS_NAN(x) || Py_IS_NAN(y)) + return Py_NAN; + if (Py_IS_INFINITY(y)) { + if (Py_IS_INFINITY(x)) { + if (copysign(1., x) == 1.) + /* atan2(+-inf, +inf) == +-pi/4 */ + return copysign(0.25*Py_MATH_PI, y); + else + /* atan2(+-inf, -inf) == +-pi*3/4 */ + return copysign(0.75*Py_MATH_PI, y); } + /* atan2(+-inf, x) == +-pi/2 for finite x */ + return copysign(0.5*Py_MATH_PI, y); } - if (PyFloat_CheckExact(args[1])) { - y = PyFloat_AS_DOUBLE(args[1]); - } - else - { - y = PyFloat_AsDouble(args[1]); - if (y == -1.0 && PyErr_Occurred()) { - goto exit; - } + if (Py_IS_INFINITY(x) || y == 0.) { + if (copysign(1., x) == 1.) + /* atan2(+-y, +inf) = atan2(+-0, +x) = +-0. */ + return copysign(0., y); + else + /* atan2(+-y, -inf) = atan2(+-0., -x) = +-pi. */ + return copysign(Py_MATH_PI, y); } - return_value = math_pow_impl(module, x, y); - -exit: - return return_value; + return atan2(y, x); } -PyDoc_STRVAR(math_degrees__doc__, -"degrees($module, x, /)\n" -"--\n" -"\n" -"Convert angle x from radians to degrees."); - -#define MATH_DEGREES_METHODDEF \ - {"degrees", (PyCFunction)math_degrees, METH_O, math_degrees__doc__}, -static PyObject * -math_degrees_impl(PyObject *module, double x); +/* IEEE 754-style remainder operation: x - n*y where n*y is the nearest + multiple of y to x, taking n even in the case of a tie. Assuming an IEEE 754 + binary floating-point format, the result is always exact. */ -static PyObject * -math_degrees(PyObject *module, PyObject *arg) +static double +m_remainder(double x, double y) { - PyObject *return_value = NULL; - double x; + /* Deal with most common case first. */ + if (Py_IS_FINITE(x) && Py_IS_FINITE(y)) { + double absx, absy, c, m, r; - if (PyFloat_CheckExact(arg)) { - x = PyFloat_AS_DOUBLE(arg); - } - else - { - x = PyFloat_AsDouble(arg); - if (x == -1.0 && PyErr_Occurred()) { - goto exit; + if (y == 0.0) { + return Py_NAN; } - } - return_value = math_degrees_impl(module, x); -exit: - return return_value; -} + absx = fabs(x); + absy = fabs(y); + m = fmod(absx, absy); + + /* + Warning: some subtlety here. What we *want* to know at this point is + whether the remainder m is less than, equal to, or greater than half + of absy. However, we can't do that comparison directly because we + can't be sure that 0.5*absy is representable (the multiplication + might incur precision loss due to underflow). So instead we compare + m with the complement c = absy - m: m < 0.5*absy if and only if m < + c, and so on. The catch is that absy - m might also not be + representable, but it turns out that it doesn't matter: + + - if m > 0.5*absy then absy - m is exactly representable, by + Sterbenz's lemma, so m > c + - if m == 0.5*absy then again absy - m is exactly representable + and m == c + - if m < 0.5*absy then either (i) 0.5*absy is exactly representable, + in which case 0.5*absy < absy - m, so 0.5*absy <= c and hence m < + c, or (ii) absy is tiny, either subnormal or in the lowest normal + binade. Then absy - m is exactly representable and again m < c. + */ + + c = absy - m; + if (m < c) { + r = m; + } + else if (m > c) { + r = -c; + } + else { + /* + Here absx is exactly halfway between two multiples of absy, + and we need to choose the even multiple. x now has the form -PyDoc_STRVAR(math_radians__doc__, -"radians($module, x, /)\n" -"--\n" -"\n" -"Convert angle x from degrees to radians."); + absx = n * absy + m -#define MATH_RADIANS_METHODDEF \ - {"radians", (PyCFunction)math_radians, METH_O, math_radians__doc__}, + for some integer n (recalling that m = 0.5*absy at this point). + If n is even we want to return m; if n is odd, we need to + return -m. -static PyObject * -math_radians_impl(PyObject *module, double x); + So -static PyObject * -math_radians(PyObject *module, PyObject *arg) -{ - PyObject *return_value = NULL; - double x; + 0.5 * (absx - m) = (n/2) * absy - if (PyFloat_CheckExact(arg)) { - x = PyFloat_AS_DOUBLE(arg); - } - else - { - x = PyFloat_AsDouble(arg); - if (x == -1.0 && PyErr_Occurred()) { - goto exit; + and now reducing modulo absy gives us: + + | m, if n is odd + fmod(0.5 * (absx - m), absy) = | + | 0, if n is even + + Now m - 2.0 * fmod(...) gives the desired result: m + if n is even, -m if m is odd. + + Note that all steps in fmod(0.5 * (absx - m), absy) + will be computed exactly, with no rounding error + introduced. + */ + assert(m == c); + r = m - 2.0 * fmod(0.5 * (absx - m), absy); } + return copysign(1.0, x) * r; } - return_value = math_radians_impl(module, x); -exit: - return return_value; + /* Special values. */ + if (Py_IS_NAN(x)) { + return x; + } + if (Py_IS_NAN(y)) { + return y; + } + if (Py_IS_INFINITY(x)) { + return Py_NAN; + } + assert(Py_IS_INFINITY(y)); + return x; } -PyDoc_STRVAR(math_isfinite__doc__, -"isfinite($module, x, /)\n" -"--\n" -"\n" -"Return True if x is neither an infinity nor a NaN, and False otherwise."); -#define MATH_ISFINITE_METHODDEF \ - {"isfinite", (PyCFunction)math_isfinite, METH_O, math_isfinite__doc__}, - -static PyObject * -math_isfinite_impl(PyObject *module, double x); +/* + Various platforms (Solaris, OpenBSD) do nonstandard things for log(0), + log(-ve), log(NaN). Here are wrappers for log and log10 that deal with + special values directly, passing positive non-special values through to + the system log/log10. + */ -static PyObject * -math_isfinite(PyObject *module, PyObject *arg) +static double +m_log(double x) { - PyObject *return_value = NULL; - double x; - - if (PyFloat_CheckExact(arg)) { - x = PyFloat_AS_DOUBLE(arg); + if (Py_IS_FINITE(x)) { + if (x > 0.0) + return log(x); + errno = EDOM; + if (x == 0.0) + return -Py_HUGE_VAL; /* log(0) = -inf */ + else + return Py_NAN; /* log(-ve) = nan */ } - else - { - x = PyFloat_AsDouble(arg); - if (x == -1.0 && PyErr_Occurred()) { - goto exit; - } + else if (Py_IS_NAN(x)) + return x; /* log(nan) = nan */ + else if (x > 0.0) + return x; /* log(inf) = inf */ + else { + errno = EDOM; + return Py_NAN; /* log(-inf) = nan */ } - return_value = math_isfinite_impl(module, x); - -exit: - return return_value; } -PyDoc_STRVAR(math_isnan__doc__, -"isnan($module, x, /)\n" -"--\n" -"\n" -"Return True if x is a NaN (not a number), and False otherwise."); +/* + log2: log to base 2. -#define MATH_ISNAN_METHODDEF \ - {"isnan", (PyCFunction)math_isnan, METH_O, math_isnan__doc__}, + Uses an algorithm that should: -static PyObject * -math_isnan_impl(PyObject *module, double x); + (a) produce exact results for powers of 2, and + (b) give a monotonic log2 (for positive finite floats), + assuming that the system log is monotonic. +*/ -static PyObject * -math_isnan(PyObject *module, PyObject *arg) +static double +m_log2(double x) { - PyObject *return_value = NULL; - double x; - - if (PyFloat_CheckExact(arg)) { - x = PyFloat_AS_DOUBLE(arg); + if (!Py_IS_FINITE(x)) { + if (Py_IS_NAN(x)) + return x; /* log2(nan) = nan */ + else if (x > 0.0) + return x; /* log2(+inf) = +inf */ + else { + errno = EDOM; + return Py_NAN; /* log2(-inf) = nan, invalid-operation */ + } } - else - { - x = PyFloat_AsDouble(arg); - if (x == -1.0 && PyErr_Occurred()) { - goto exit; + + if (x > 0.0) { +#ifdef HAVE_LOG2 + return log2(x); +#else + double m; + int e; + m = frexp(x, &e); + /* We want log2(m * 2**e) == log(m) / log(2) + e. Care is needed when + * x is just greater than 1.0: in that case e is 1, log(m) is negative, + * and we get significant cancellation error from the addition of + * log(m) / log(2) to e. The slight rewrite of the expression below + * avoids this problem. + */ + if (x >= 1.0) { + return log(2.0 * m) / log(2.0) + (e - 1); } + else { + return log(m) / log(2.0) + e; + } +#endif + } + else if (x == 0.0) { + errno = EDOM; + return -Py_HUGE_VAL; /* log2(0) = -inf, divide-by-zero */ } - return_value = math_isnan_impl(module, x); + else { + errno = EDOM; + return Py_NAN; /* log2(-inf) = nan, invalid-operation */ + } +} -exit: - return return_value; +static double +m_log10(double x) +{ + if (Py_IS_FINITE(x)) { + if (x > 0.0) + return log10(x); + errno = EDOM; + if (x == 0.0) + return -Py_HUGE_VAL; /* log10(0) = -inf */ + else + return Py_NAN; /* log10(-ve) = nan */ + } + else if (Py_IS_NAN(x)) + return x; /* log10(nan) = nan */ + else if (x > 0.0) + return x; /* log10(inf) = inf */ + else { + errno = EDOM; + return Py_NAN; /* log10(-inf) = nan */ + } } -PyDoc_STRVAR(math_isinf__doc__, -"isinf($module, x, /)\n" -"--\n" -"\n" -"Return True if x is a positive or negative infinity, and False otherwise."); -#define MATH_ISINF_METHODDEF \ - {"isinf", (PyCFunction)math_isinf, METH_O, math_isinf__doc__}, +/*[clinic input] +math.gcd + + x as a: object + y as b: object + / -static PyObject * -math_isinf_impl(PyObject *module, double x); +greatest common divisor of x and y +[clinic start generated code]*/ static PyObject * -math_isinf(PyObject *module, PyObject *arg) +math_gcd_impl(PyObject *module, PyObject *a, PyObject *b) +/*[clinic end generated code: output=7b2e0c151bd7a5d8 input=c2691e57fb2a98fa]*/ { - PyObject *return_value = NULL; - double x; - - if (PyFloat_CheckExact(arg)) { - x = PyFloat_AS_DOUBLE(arg); + PyObject *g; + + a = PyNumber_Index(a); + if (a == NULL) + return NULL; + b = PyNumber_Index(b); + if (b == NULL) { + Py_DECREF(a); + return NULL; } - else - { - x = PyFloat_AsDouble(arg); - if (x == -1.0 && PyErr_Occurred()) { - goto exit; - } - } - return_value = math_isinf_impl(module, x); - -exit: - return return_value; -} - -PyDoc_STRVAR(math_isclose__doc__, -"isclose($module, /, a, b, *, rel_tol=1e-09, abs_tol=0.0)\n" -"--\n" -"\n" -"Determine whether two floating point numbers are close in value.\n" -"\n" -" rel_tol\n" -" maximum difference for being considered \"close\", relative to the\n" -" magnitude of the input values\n" -" abs_tol\n" -" maximum difference for being considered \"close\", regardless of the\n" -" magnitude of the input values\n" -"\n" -"Return True if a is close in value to b, and False otherwise.\n" -"\n" -"For the values to be considered close, the difference between them\n" -"must be smaller than at least one of the tolerances.\n" -"\n" -"-inf, inf and NaN behave similarly to the IEEE 754 Standard. That\n" -"is, NaN is not close to anything, even itself. inf and -inf are\n" -"only close to themselves."); - -#define MATH_ISCLOSE_METHODDEF \ - {"isclose", (PyCFunction)(void(*)(void))math_isclose, METH_FASTCALL|METH_KEYWORDS, math_isclose__doc__}, + g = _PyLong_GCD(a, b); + Py_DECREF(a); + Py_DECREF(b); + return g; +} -static int -math_isclose_impl(PyObject *module, double a, double b, double rel_tol, - double abs_tol); -static PyObject * -math_isclose(PyObject *module, PyObject *const *args, Py_ssize_t nargs, PyObject *kwnames) +/* Call is_error when errno != 0, and where x is the result libm + * returned. is_error will usually set up an exception and return + * true (1), but may return false (0) without setting up an exception. + */ +static int +is_error(double x) { - PyObject *return_value = NULL; - static const char * const _keywords[] = {"a", "b", "rel_tol", "abs_tol", NULL}; - static _PyArg_Parser _parser = {NULL, _keywords, "isclose", 0}; - PyObject *argsbuf[4]; - Py_ssize_t noptargs = nargs + (kwnames ? PyTuple_GET_SIZE(kwnames) : 0) - 2; - double a; - double b; - double rel_tol = 1e-09; - double abs_tol = 0.0; - int _return_value; - - args = _PyArg_UnpackKeywords(args, nargs, NULL, kwnames, &_parser, 2, 2, 0, argsbuf); - if (!args) { - goto exit; - } - if (PyFloat_CheckExact(args[0])) { - a = PyFloat_AS_DOUBLE(args[0]); - } - else - { - a = PyFloat_AsDouble(args[0]); - if (a == -1.0 && PyErr_Occurred()) { - goto exit; - } - } - if (PyFloat_CheckExact(args[1])) { - b = PyFloat_AS_DOUBLE(args[1]); - } - else - { - b = PyFloat_AsDouble(args[1]); - if (b == -1.0 && PyErr_Occurred()) { - goto exit; - } - } - if (!noptargs) { - goto skip_optional_kwonly; - } - if (args[2]) { - if (PyFloat_CheckExact(args[2])) { - rel_tol = PyFloat_AS_DOUBLE(args[2]); - } + int result = 1; /* presumption of guilt */ + assert(errno); /* non-zero errno is a precondition for calling */ + if (errno == EDOM) + PyErr_SetString(PyExc_ValueError, "math domain error"); + + else if (errno == ERANGE) { + /* ANSI C generally requires libm functions to set ERANGE + * on overflow, but also generally *allows* them to set + * ERANGE on underflow too. There's no consistency about + * the latter across platforms. + * Alas, C99 never requires that errno be set. + * Here we suppress the underflow errors (libm functions + * should return a zero on underflow, and +- HUGE_VAL on + * overflow, so testing the result for zero suffices to + * distinguish the cases). + * + * On some platforms (Ubuntu/ia64) it seems that errno can be + * set to ERANGE for subnormal results that do *not* underflow + * to zero. So to be safe, we'll ignore ERANGE whenever the + * function result is less than one in absolute value. + */ + if (fabs(x) < 1.0) + result = 0; else - { - rel_tol = PyFloat_AsDouble(args[2]); - if (rel_tol == -1.0 && PyErr_Occurred()) { - goto exit; - } - } - if (!--noptargs) { - goto skip_optional_kwonly; - } - } - if (PyFloat_CheckExact(args[3])) { - abs_tol = PyFloat_AS_DOUBLE(args[3]); + PyErr_SetString(PyExc_OverflowError, + "math range error"); } else - { - abs_tol = PyFloat_AsDouble(args[3]); - if (abs_tol == -1.0 && PyErr_Occurred()) { - goto exit; - } - } -skip_optional_kwonly: - _return_value = math_isclose_impl(module, a, b, rel_tol, abs_tol); - if ((_return_value == -1) && PyErr_Occurred()) { - goto exit; - } - return_value = PyBool_FromLong((long)_return_value); - -exit: - return return_value; + /* Unexpected math error */ + PyErr_SetFromErrno(PyExc_ValueError); + return result; } -PyDoc_STRVAR(math_prod__doc__, -"prod($module, iterable, /, *, start=1)\n" -"--\n" -"\n" -"Calculate the product of all the elements in the input iterable.\n" -"\n" -"The default start value for the product is 1.\n" -"\n" -"When the iterable is empty, return the start value. This function is\n" -"intended specifically for use with numeric values and may reject\n" -"non-numeric types."); - -#define MATH_PROD_METHODDEF \ - {"prod", (PyCFunction)(void(*)(void))math_prod, METH_FASTCALL|METH_KEYWORDS, math_prod__doc__}, - -static PyObject * -math_prod_impl(PyObject *module, PyObject *iterable, PyObject *start); +/* + math_1 is used to wrap a libm function f that takes a double + argument and returns a double. + + The error reporting follows these rules, which are designed to do + the right thing on C89/C99 platforms and IEEE 754/non IEEE 754 + platforms. + + - a NaN result from non-NaN inputs causes ValueError to be raised + - an infinite result from finite inputs causes OverflowError to be + raised if can_overflow is 1, or raises ValueError if can_overflow + is 0. + - if the result is finite and errno == EDOM then ValueError is + raised + - if the result is finite and nonzero and errno == ERANGE then + OverflowError is raised + + The last rule is used to catch overflow on platforms which follow + C89 but for which HUGE_VAL is not an infinity. + + For the majority of one-argument functions these rules are enough + to ensure that Python's functions behave as specified in 'Annex F' + of the C99 standard, with the 'invalid' and 'divide-by-zero' + floating-point exceptions mapping to Python's ValueError and the + 'overflow' floating-point exception mapping to OverflowError. + math_1 only works for functions that don't have singularities *and* + the possibility of overflow; fortunately, that covers everything we + care about right now. +*/ static PyObject * -math_prod(PyObject *module, PyObject *const *args, Py_ssize_t nargs, PyObject *kwnames) +math_1_to_whatever(PyObject *arg, double (*func) (double), + PyObject *(*from_double_func) (double), + int can_overflow) { - PyObject *return_value = NULL; - static const char * const _keywords[] = {"", "start", NULL}; - static _PyArg_Parser _parser = {NULL, _keywords, "prod", 0}; - PyObject *argsbuf[2]; - Py_ssize_t noptargs = nargs + (kwnames ? PyTuple_GET_SIZE(kwnames) : 0) - 1; - PyObject *iterable; - PyObject *start = NULL; - - args = _PyArg_UnpackKeywords(args, nargs, NULL, kwnames, &_parser, 1, 1, 0, argsbuf); - if (!args) { - goto exit; + double x, r; + x = PyFloat_AsDouble(arg); + if (x == -1.0 && PyErr_Occurred()) + return NULL; + errno = 0; + r = (*func)(x); + if (Py_IS_NAN(r) && !Py_IS_NAN(x)) { + PyErr_SetString(PyExc_ValueError, + "math domain error"); /* invalid arg */ + return NULL; } - iterable = args[0]; - if (!noptargs) { - goto skip_optional_kwonly; + if (Py_IS_INFINITY(r) && Py_IS_FINITE(x)) { + if (can_overflow) + PyErr_SetString(PyExc_OverflowError, + "math range error"); /* overflow */ + else + PyErr_SetString(PyExc_ValueError, + "math domain error"); /* singularity */ + return NULL; } - start = args[1]; -skip_optional_kwonly: - return_value = math_prod_impl(module, iterable, start); + if (Py_IS_FINITE(r) && errno && is_error(r)) + /* this branch unnecessary on most platforms */ + return NULL; -exit: - return return_value; + return (*from_double_func)(r); } -PyDoc_STRVAR(math_perm__doc__, -"perm($module, n, k=None, /)\n" -"--\n" -"\n" -"Number of ways to choose k items from n items without repetition and with order.\n" -"\n" -"Evaluates to n! / (n - k)! when k <= n and evaluates\n" -"to zero when k > n.\n" -"\n" -"If k is not specified or is None, then k defaults to n\n" -"and the function returns n!.\n" -"\n" -"Raises TypeError if either of the arguments are not integers.\n" -"Raises ValueError if either of the arguments are negative."); - -#define MATH_PERM_METHODDEF \ - {"perm", (PyCFunction)(void(*)(void))math_perm, METH_FASTCALL, math_perm__doc__}, +/* variant of math_1, to be used when the function being wrapped is known to + set errno properly (that is, errno = EDOM for invalid or divide-by-zero, + errno = ERANGE for overflow). */ static PyObject * -math_perm_impl(PyObject *module, PyObject *n, PyObject *k); +math_1a(PyObject *arg, double (*func) (double)) +{ + double x, r; + x = PyFloat_AsDouble(arg); + if (x == -1.0 && PyErr_Occurred()) + return NULL; + errno = 0; + r = (*func)(x); + if (errno && is_error(r)) + return NULL; + return PyFloat_FromDouble(r); +} + +/* + math_2 is used to wrap a libm function f that takes two double + arguments and returns a double. + + The error reporting follows these rules, which are designed to do + the right thing on C89/C99 platforms and IEEE 754/non IEEE 754 + platforms. + + - a NaN result from non-NaN inputs causes ValueError to be raised + - an infinite result from finite inputs causes OverflowError to be + raised. + - if the result is finite and errno == EDOM then ValueError is + raised + - if the result is finite and nonzero and errno == ERANGE then + OverflowError is raised + + The last rule is used to catch overflow on platforms which follow + C89 but for which HUGE_VAL is not an infinity. + + For most two-argument functions (copysign, fmod, hypot, atan2) + these rules are enough to ensure that Python's functions behave as + specified in 'Annex F' of the C99 standard, with the 'invalid' and + 'divide-by-zero' floating-point exceptions mapping to Python's + ValueError and the 'overflow' floating-point exception mapping to + OverflowError. +*/ static PyObject * -math_perm(PyObject *module, PyObject *const *args, Py_ssize_t nargs) +math_1(PyObject *arg, double (*func) (double), int can_overflow) { - PyObject *return_value = NULL; - PyObject *n; - PyObject *k = Py_None; + return math_1_to_whatever(arg, func, PyFloat_FromDouble, can_overflow); +} - if (!_PyArg_CheckPositional("perm", nargs, 1, 2)) { - goto exit; +static PyObject * +math_2(PyObject *const *args, Py_ssize_t nargs, + double (*func) (double, double), const char *funcname) +{ + double x, y, r; + if (!_PyArg_CheckPositional(funcname, nargs, 2, 2)) + return NULL; + x = PyFloat_AsDouble(args[0]); + y = PyFloat_AsDouble(args[1]); + if ((x == -1.0 || y == -1.0) && PyErr_Occurred()) + return NULL; + errno = 0; + r = (*func)(x, y); + if (Py_IS_NAN(r)) { + if (!Py_IS_NAN(x) && !Py_IS_NAN(y)) + errno = EDOM; + else + errno = 0; } - n = args[0]; - if (nargs < 2) { - goto skip_optional; + else if (Py_IS_INFINITY(r)) { + if (Py_IS_FINITE(x) && Py_IS_FINITE(y)) + errno = ERANGE; + else + errno = 0; } - k = args[1]; -skip_optional: - return_value = math_perm_impl(module, n, k); - -exit: - return return_value; + if (errno && is_error(r)) + return NULL; + else + return PyFloat_FromDouble(r); } -PyDoc_STRVAR(math_comb__doc__, -"comb($module, n, k, /)\n" -"--\n" -"\n" -"Number of ways to choose k items from n items without repetition and without order.\n" -"\n" -"Evaluates to n! / (k! * (n - k)!) when k <= n and evaluates\n" -"to zero when k > n.\n" -"\n" -"Also called the binomial coefficient because it is equivalent\n" -"to the coefficient of k-th term in polynomial expansion of the\n" -"expression (1 + x)**n.\n" -"\n" -"Raises TypeError if either of the arguments are not integers.\n" -"Raises ValueError if either of the arguments are negative."); +#define FUNC1(funcname, func, can_overflow, docstring) \ + static PyObject * math_##funcname(PyObject *self, PyObject *args) { \ + return math_1(args, func, can_overflow); \ + }\ + PyDoc_STRVAR(math_##funcname##_doc, docstring); + +#define FUNC1A(funcname, func, docstring) \ + static PyObject * math_##funcname(PyObject *self, PyObject *args) { \ + return math_1a(args, func); \ + }\ + PyDoc_STRVAR(math_##funcname##_doc, docstring); + +#define FUNC2(funcname, func, docstring) \ + static PyObject * math_##funcname(PyObject *self, PyObject *const *args, Py_ssize_t nargs) { \ + return math_2(args, nargs, func, #funcname); \ + }\ + PyDoc_STRVAR(math_##funcname##_doc, docstring); + +FUNC1(acos, acos, 0, + "acos($module, x, /)\n--\n\n" + "Return the arc cosine (measured in radians) of x.\n\n" + "The result is between 0 and pi.") +FUNC1(acosh, m_acosh, 0, + "acosh($module, x, /)\n--\n\n" + "Return the inverse hyperbolic cosine of x.") +FUNC1(asin, asin, 0, + "asin($module, x, /)\n--\n\n" + "Return the arc sine (measured in radians) of x.\n\n" + "The result is between -pi/2 and pi/2.") +FUNC1(asinh, m_asinh, 0, + "asinh($module, x, /)\n--\n\n" + "Return the inverse hyperbolic sine of x.") +FUNC1(atan, atan, 0, + "atan($module, x, /)\n--\n\n" + "Return the arc tangent (measured in radians) of x.\n\n" + "The result is between -pi/2 and pi/2.") +FUNC2(atan2, m_atan2, + "atan2($module, y, x, /)\n--\n\n" + "Return the arc tangent (measured in radians) of y/x.\n\n" + "Unlike atan(y/x), the signs of both x and y are considered.") +FUNC1(atanh, m_atanh, 0, + "atanh($module, x, /)\n--\n\n" + "Return the inverse hyperbolic tangent of x.") -#define MATH_COMB_METHODDEF \ - {"comb", (PyCFunction)(void(*)(void))math_comb, METH_FASTCALL, math_comb__doc__}, +/*[clinic input] +math.ceil -static PyObject * -math_comb_impl(PyObject *module, PyObject *n, PyObject *k); + x as number: object + / + +Return the ceiling of x as an Integral. + +This is the smallest integer >= x. +[clinic start generated code]*/ static PyObject * -math_comb(PyObject *module, PyObject *const *args, Py_ssize_t nargs) +math_ceil(PyObject *module, PyObject *number) +/*[clinic end generated code: output=6c3b8a78bc201c67 input=2725352806399cab]*/ { - PyObject *return_value = NULL; - PyObject *n; - PyObject *k; - - if (!_PyArg_CheckPositional("comb", nargs, 2, 2)) { - goto exit; + _Py_IDENTIFIER(__ceil__); + + if (!PyFloat_CheckExact(number)) { + PyObject *method = _PyObject_LookupSpecial(number, &PyId___ceil__); + if (method != NULL) { + PyObject *result = _PyObject_CallNoArg(method); + Py_DECREF(method); + return result; + } + if (PyErr_Occurred()) + return NULL; } - n = args[0]; - k = args[1]; - return_value = math_comb_impl(module, n, k); + double x = PyFloat_AsDouble(number); + if (x == -1.0 && PyErr_Occurred()) + return NULL; -exit: - return return_value; + return PyLong_FromDouble(ceil(x)); } -PyDoc_STRVAR(math_nextafter__doc__, -"nextafter($module, x, y, /)\n" -"--\n" -"\n" -"Return the next floating-point value after x towards y."); +FUNC2(copysign, copysign, + "copysign($module, x, y, /)\n--\n\n" + "Return a float with the magnitude (absolute value) of x but the sign of y.\n\n" + "On platforms that support signed zeros, copysign(1.0, -0.0)\n" + "returns -1.0.\n") +FUNC1(cos, cos, 0, + "cos($module, x, /)\n--\n\n" + "Return the cosine of x (measured in radians).") +FUNC1(cosh, cosh, 1, + "cosh($module, x, /)\n--\n\n" + "Return the hyperbolic cosine of x.") +FUNC1A(erf, m_erf, + "erf($module, x, /)\n--\n\n" + "Error function at x.") +FUNC1A(erfc, m_erfc, + "erfc($module, x, /)\n--\n\n" + "Complementary error function at x.") +FUNC1(exp, exp, 1, + "exp($module, x, /)\n--\n\n" + "Return e raised to the power of x.") +FUNC1(expm1, m_expm1, 1, + "expm1($module, x, /)\n--\n\n" + "Return exp(x)-1.\n\n" + "This function avoids the loss of precision involved in the direct " + "evaluation of exp(x)-1 for small x.") +FUNC1(fabs, fabs, 0, + "fabs($module, x, /)\n--\n\n" + "Return the absolute value of the float x.") + +/*[clinic input] +math.floor -#define MATH_NEXTAFTER_METHODDEF \ - {"nextafter", (PyCFunction)(void(*)(void))math_nextafter, METH_FASTCALL, math_nextafter__doc__}, + x as number: object + / -static PyObject * -math_nextafter_impl(PyObject *module, double x, double y); +Return the floor of x as an Integral. + +This is the largest integer <= x. +[clinic start generated code]*/ static PyObject * -math_nextafter(PyObject *module, PyObject *const *args, Py_ssize_t nargs) +math_floor(PyObject *module, PyObject *number) +/*[clinic end generated code: output=c6a65c4884884b8a input=63af6b5d7ebcc3d6]*/ { - PyObject *return_value = NULL; - double x; - double y; - - if (!_PyArg_CheckPositional("nextafter", nargs, 2, 2)) { - goto exit; - } - if (PyFloat_CheckExact(args[0])) { - x = PyFloat_AS_DOUBLE(args[0]); - } - else - { - x = PyFloat_AsDouble(args[0]); - if (x == -1.0 && PyErr_Occurred()) { - goto exit; + _Py_IDENTIFIER(__floor__); + + if (!PyFloat_CheckExact(number)) { + PyObject *method = _PyObject_LookupSpecial(number, &PyId___floor__); + if (method != NULL) { + PyObject *result = _PyObject_CallNoArg(method); + Py_DECREF(method); + return result; } + if (PyErr_Occurred()) + return NULL; } - if (PyFloat_CheckExact(args[1])) { - y = PyFloat_AS_DOUBLE(args[1]); - } - else - { - y = PyFloat_AsDouble(args[1]); - if (y == -1.0 && PyErr_Occurred()) { - goto exit; + double x = PyFloat_AsDouble(number); + if (x == -1.0 && PyErr_Occurred()) + return NULL; + + return PyLong_FromDouble(floor(x)); +} + +FUNC1A(gamma, m_tgamma, + "gamma($module, x, /)\n--\n\n" + "Gamma function at x.") +FUNC1A(lgamma, m_lgamma, + "lgamma($module, x, /)\n--\n\n" + "Natural logarithm of absolute value of Gamma function at x.") +FUNC1(log1p, m_log1p, 0, + "log1p($module, x, /)\n--\n\n" + "Return the natural logarithm of 1+x (base e).\n\n" + "The result is computed in a way which is accurate for x near zero.") +FUNC2(remainder, m_remainder, + "remainder($module, x, y, /)\n--\n\n" + "Difference between x and the closest integer multiple of y.\n\n" + "Return x - n*y where n*y is the closest integer multiple of y.\n" + "In the case where x is exactly halfway between two multiples of\n" + "y, the nearest even value of n is used. The result is always exact.") +FUNC1(sin, sin, 0, + "sin($module, x, /)\n--\n\n" + "Return the sine of x (measured in radians).") +FUNC1(sinh, sinh, 1, + "sinh($module, x, /)\n--\n\n" + "Return the hyperbolic sine of x.") +FUNC1(sqrt, sqrt, 0, + "sqrt($module, x, /)\n--\n\n" + "Return the square root of x.") +FUNC1(tan, tan, 0, + "tan($module, x, /)\n--\n\n" + "Return the tangent of x (measured in radians).") +FUNC1(tanh, tanh, 0, + "tanh($module, x, /)\n--\n\n" + "Return the hyperbolic tangent of x.") + +/* Precision summation function as msum() by Raymond Hettinger in + , + enhanced with the exact partials sum and roundoff from Mark + Dickinson's post at . + See those links for more details, proofs and other references. + + Note 1: IEEE 754R floating point semantics are assumed, + but the current implementation does not re-establish special + value semantics across iterations (i.e. handling -Inf + Inf). + + Note 2: No provision is made for intermediate overflow handling; + therefore, sum([1e+308, 1e-308, 1e+308]) returns 1e+308 while + sum([1e+308, 1e+308, 1e-308]) raises an OverflowError due to the + overflow of the first partial sum. + + Note 3: The intermediate values lo, yr, and hi are declared volatile so + aggressive compilers won't algebraically reduce lo to always be exactly 0.0. + Also, the volatile declaration forces the values to be stored in memory as + regular doubles instead of extended long precision (80-bit) values. This + prevents double rounding because any addition or subtraction of two doubles + can be resolved exactly into double-sized hi and lo values. As long as the + hi value gets forced into a double before yr and lo are computed, the extra + bits in downstream extended precision operations (x87 for example) will be + exactly zero and therefore can be losslessly stored back into a double, + thereby preventing double rounding. + + Note 4: A similar implementation is in Modules/cmathmodule.c. + Be sure to update both when making changes. + + Note 5: The signature of math.fsum() differs from builtins.sum() + because the start argument doesn't make sense in the context of + accurate summation. Since the partials table is collapsed before + returning a result, sum(seq2, start=sum(seq1)) may not equal the + accurate result returned by sum(itertools.chain(seq1, seq2)). +*/ + +#define NUM_PARTIALS 32 /* initial partials array size, on stack */ + +/* Extend the partials array p[] by doubling its size. */ +static int /* non-zero on error */ +_fsum_realloc(double **p_ptr, Py_ssize_t n, + double *ps, Py_ssize_t *m_ptr) +{ + void *v = NULL; + Py_ssize_t m = *m_ptr; + + m += m; /* double */ + if (n < m && (size_t)m < ((size_t)PY_SSIZE_T_MAX / sizeof(double))) { + double *p = *p_ptr; + if (p == ps) { + v = PyMem_Malloc(sizeof(double) * m); + if (v != NULL) + memcpy(v, ps, sizeof(double) * n); } + else + v = PyMem_Realloc(p, sizeof(double) * m); } - return_value = math_nextafter_impl(module, x, y); - -exit: - return return_value; + if (v == NULL) { /* size overflow or no memory */ + PyErr_SetString(PyExc_MemoryError, "math.fsum partials"); + return 1; + } + *p_ptr = (double*) v; + *m_ptr = m; + return 0; } -PyDoc_STRVAR(math_ulp__doc__, -"ulp($module, x, /)\n" -"--\n" -"\n" -"Return the value of the least significant bit of the float x."); +/* Full precision summation of a sequence of floats. + + def msum(iterable): + partials = [] # sorted, non-overlapping partial sums + for x in iterable: + i = 0 + for y in partials: + if abs(x) < abs(y): + x, y = y, x + hi = x + y + lo = y - (hi - x) + if lo: + partials[i] = lo + i += 1 + x = hi + partials[i:] = [x] + return sum_exact(partials) + + Rounded x+y stored in hi with the roundoff stored in lo. Together hi+lo + are exactly equal to x+y. The inner loop applies hi/lo summation to each + partial so that the list of partial sums remains exact. + + Sum_exact() adds the partial sums exactly and correctly rounds the final + result (using the round-half-to-even rule). The items in partials remain + non-zero, non-special, non-overlapping and strictly increasing in + magnitude, but possibly not all having the same sign. + + Depends on IEEE 754 arithmetic guarantees and half-even rounding. +*/ + +/*[clinic input] +math.fsum -#define MATH_ULP_METHODDEF \ - {"ulp", (PyCFunction)math_ulp, METH_O, math_ulp__doc__}, + seq: object + / -static double -math_ulp_impl(PyObject *module, double x); +Return an accurate floating point sum of values in the iterable seq. + +Assumes IEEE-754 floating point arithmetic. +[clinic start generated code]*/ static PyObject * -math_ulp(PyObject *module, PyObject *arg) +math_fsum(PyObject *module, PyObject *seq) +/*[clinic end generated code: output=ba5c672b87fe34fc input=c51b7d8caf6f6e82]*/ { - PyObject *return_value = NULL; - double x; - double _return_value; - - if (PyFloat_CheckExact(arg)) { - x = PyFloat_AS_DOUBLE(arg); - } + PyObject *item, *iter, *sum = NULL; + Py_ssize_t i, j, n = 0, m = NUM_PARTIALS; + double x, y, t, ps[NUM_PARTIALS], *p = ps; + double xsave, special_sum = 0.0, inf_sum = 0.0; + volatile double hi, yr, lo; + + iter = PyObject_GetIter(seq); + if (iter == NULL) + return NULL; + + for(;;) { /* for x in iterable */ + assert(0 <= n && n <= m); + assert((m == NUM_PARTIALS && p == ps) || + (m > NUM_PARTIALS && p != NULL)); + + item = PyIter_Next(iter); + if (item == NULL) { + if (PyErr_Occurred()) + goto _fsum_error; + break; + } + ASSIGN_DOUBLE(x, item, error_with_item); + Py_DECREF(item); + + xsave = x; + for (i = j = 0; j < n; j++) { /* for y in partials */ + y = p[j]; + if (fabs(x) < fabs(y)) { + t = x; x = y; y = t; + } + hi = x + y; + yr = hi - x; + lo = y - yr; + if (lo != 0.0) + p[i++] = lo; + x = hi; + } + + n = i; /* ps[i:] = [x] */ + if (x != 0.0) { + if (! Py_IS_FINITE(x)) { + /* a nonfinite x could arise either as + a result of intermediate overflow, or + as a result of a nan or inf in the + summands */ + if (Py_IS_FINITE(xsave)) { + PyErr_SetString(PyExc_OverflowError, + "intermediate overflow in fsum"); + goto _fsum_error; + } + if (Py_IS_INFINITY(xsave)) + inf_sum += xsave; + special_sum += xsave; + /* reset partials */ + n = 0; + } + else if (n >= m && _fsum_realloc(&p, n, ps, &m)) + goto _fsum_error; + else + p[n++] = x; + } + } + + if (special_sum != 0.0) { + if (Py_IS_NAN(inf_sum)) + PyErr_SetString(PyExc_ValueError, + "-inf + inf in fsum"); + else + sum = PyFloat_FromDouble(special_sum); + goto _fsum_error; + } + + hi = 0.0; + if (n > 0) { + hi = p[--n]; + /* sum_exact(ps, hi) from the top, stop when the sum becomes + inexact. */ + while (n > 0) { + x = hi; + y = p[--n]; + assert(fabs(y) < fabs(x)); + hi = x + y; + yr = hi - x; + lo = y - yr; + if (lo != 0.0) + break; + } + /* Make half-even rounding work across multiple partials. + Needed so that sum([1e-16, 1, 1e16]) will round-up the last + digit to two instead of down to zero (the 1e-16 makes the 1 + slightly closer to two). With a potential 1 ULP rounding + error fixed-up, math.fsum() can guarantee commutativity. */ + if (n > 0 && ((lo < 0.0 && p[n-1] < 0.0) || + (lo > 0.0 && p[n-1] > 0.0))) { + y = lo * 2.0; + x = hi + y; + yr = x - hi; + if (y == yr) + hi = x; + } + } + sum = PyFloat_FromDouble(hi); + + _fsum_error: + Py_DECREF(iter); + if (p != ps) + PyMem_Free(p); + return sum; + + error_with_item: + Py_DECREF(item); + goto _fsum_error; +} + +#undef NUM_PARTIALS + + +static unsigned long +count_set_bits(unsigned long n) +{ + unsigned long count = 0; + while (n != 0) { + ++count; + n &= n - 1; /* clear least significant bit */ + } + return count; +} + +/* Integer square root + +Given a nonnegative integer `n`, we want to compute the largest integer +`a` for which `a * a <= n`, or equivalently the integer part of the exact +square root of `n`. + +We use an adaptive-precision pure-integer version of Newton's iteration. Given +a positive integer `n`, the algorithm produces at each iteration an integer +approximation `a` to the square root of `n >> s` for some even integer `s`, +with `s` decreasing as the iterations progress. On the final iteration, `s` is +zero and we have an approximation to the square root of `n` itself. + +At every step, the approximation `a` is strictly within 1.0 of the true square +root, so we have + + (a - 1)**2 < (n >> s) < (a + 1)**2 + +After the final iteration, a check-and-correct step is needed to determine +whether `a` or `a - 1` gives the desired integer square root of `n`. + +The algorithm is remarkable in its simplicity. There's no need for a +per-iteration check-and-correct step, and termination is straightforward: the +number of iterations is known in advance (it's exactly `floor(log2(log2(n)))` +for `n > 1`). The only tricky part of the correctness proof is in establishing +that the bound `(a - 1)**2 < (n >> s) < (a + 1)**2` is maintained from one +iteration to the next. A sketch of the proof of this is given below. + +In addition to the proof sketch, a formal, computer-verified proof +of correctness (using Lean) of an equivalent recursive algorithm can be found +here: + + https://github.com/mdickinson/snippets/blob/master/proofs/isqrt/src/isqrt.lean + + +Here's Python code equivalent to the C implementation below: + + def isqrt(n): + """ + Return the integer part of the square root of the input. + """ + n = operator.index(n) + + if n < 0: + raise ValueError("isqrt() argument must be nonnegative") + if n == 0: + return 0 + + c = (n.bit_length() - 1) // 2 + a = 1 + d = 0 + for s in reversed(range(c.bit_length())): + # Loop invariant: (a-1)**2 < (n >> 2*(c - d)) < (a+1)**2 + e = d + d = c >> s + a = (a << d - e - 1) + (n >> 2*c - e - d + 1) // a + + return a - (a*a > n) + + +Sketch of proof of correctness +------------------------------ + +The delicate part of the correctness proof is showing that the loop invariant +is preserved from one iteration to the next. That is, just before the line + + a = (a << d - e - 1) + (n >> 2*c - e - d + 1) // a + +is executed in the above code, we know that + + (1) (a - 1)**2 < (n >> 2*(c - e)) < (a + 1)**2. + +(since `e` is always the value of `d` from the previous iteration). We must +prove that after that line is executed, we have + + (a - 1)**2 < (n >> 2*(c - d)) < (a + 1)**2 + +To facilitate the proof, we make some changes of notation. Write `m` for +`n >> 2*(c-d)`, and write `b` for the new value of `a`, so + + b = (a << d - e - 1) + (n >> 2*c - e - d + 1) // a + +or equivalently: + + (2) b = (a << d - e - 1) + (m >> d - e + 1) // a + +Then we can rewrite (1) as: + + (3) (a - 1)**2 < (m >> 2*(d - e)) < (a + 1)**2 + +and we must show that (b - 1)**2 < m < (b + 1)**2. + +From this point on, we switch to mathematical notation, so `/` means exact +division rather than integer division and `^` is used for exponentiation. We +use the `√` symbol for the exact square root. In (3), we can remove the +implicit floor operation to give: + + (4) (a - 1)^2 < m / 4^(d - e) < (a + 1)^2 + +Taking square roots throughout (4), scaling by `2^(d-e)`, and rearranging gives + + (5) 0 <= | 2^(d-e)a - √m | < 2^(d-e) + +Squaring and dividing through by `2^(d-e+1) a` gives + + (6) 0 <= 2^(d-e-1) a + m / (2^(d-e+1) a) - √m < 2^(d-e-1) / a + +We'll show below that `2^(d-e-1) <= a`. Given that, we can replace the +right-hand side of (6) with `1`, and now replacing the central +term `m / (2^(d-e+1) a)` with its floor in (6) gives + + (7) -1 < 2^(d-e-1) a + m // 2^(d-e+1) a - √m < 1 + +Or equivalently, from (2): + + (7) -1 < b - √m < 1 + +and rearranging gives that `(b-1)^2 < m < (b+1)^2`, which is what we needed +to prove. + +We're not quite done: we still have to prove the inequality `2^(d - e - 1) <= +a` that was used to get line (7) above. From the definition of `c`, we have +`4^c <= n`, which implies + + (8) 4^d <= m + +also, since `e == d >> 1`, `d` is at most `2e + 1`, from which it follows +that `2d - 2e - 1 <= d` and hence that + + (9) 4^(2d - 2e - 1) <= m + +Dividing both sides by `4^(d - e)` gives + + (10) 4^(d - e - 1) <= m / 4^(d - e) + +But we know from (4) that `m / 4^(d-e) < (a + 1)^2`, hence + + (11) 4^(d - e - 1) < (a + 1)^2 + +Now taking square roots of both sides and observing that both `2^(d-e-1)` and +`a` are integers gives `2^(d - e - 1) <= a`, which is what we needed. This +completes the proof sketch. + +*/ + + +/* Approximate square root of a large 64-bit integer. + + Given `n` satisfying `2**62 <= n < 2**64`, return `a` + satisfying `(a - 1)**2 < n < (a + 1)**2`. */ + +static uint64_t +_approximate_isqrt(uint64_t n) +{ + uint32_t u = 1U + (n >> 62); + u = (u << 1) + (n >> 59) / u; + u = (u << 3) + (n >> 53) / u; + u = (u << 7) + (n >> 41) / u; + return (u << 15) + (n >> 17) / u; +} + +/*[clinic input] +math.isqrt + + n: object + / + +Return the integer part of the square root of the input. +[clinic start generated code]*/ + +static PyObject * +math_isqrt(PyObject *module, PyObject *n) +/*[clinic end generated code: output=35a6f7f980beab26 input=5b6e7ae4fa6c43d6]*/ +{ + int a_too_large, c_bit_length; + size_t c, d; + uint64_t m, u; + PyObject *a = NULL, *b; + + n = PyNumber_Index(n); + if (n == NULL) { + return NULL; + } + + if (_PyLong_Sign(n) < 0) { + PyErr_SetString( + PyExc_ValueError, + "isqrt() argument must be nonnegative"); + goto error; + } + if (_PyLong_Sign(n) == 0) { + Py_DECREF(n); + return PyLong_FromLong(0); + } + + /* c = (n.bit_length() - 1) // 2 */ + c = _PyLong_NumBits(n); + if (c == (size_t)(-1)) { + goto error; + } + c = (c - 1U) / 2U; + + /* Fast path: if c <= 31 then n < 2**64 and we can compute directly with a + fast, almost branch-free algorithm. In the final correction, we use `u*u + - 1 >= m` instead of the simpler `u*u > m` in order to get the correct + result in the corner case where `u=2**32`. */ + if (c <= 31U) { + m = (uint64_t)PyLong_AsUnsignedLongLong(n); + Py_DECREF(n); + if (m == (uint64_t)(-1) && PyErr_Occurred()) { + return NULL; + } + u = _approximate_isqrt(m << (62U - 2U*c)) >> (31U - c); + u -= u * u - 1U >= m; + return PyLong_FromUnsignedLongLong((unsigned long long)u); + } + + /* Slow path: n >= 2**64. We perform the first five iterations in C integer + arithmetic, then switch to using Python long integers. */ + + /* From n >= 2**64 it follows that c.bit_length() >= 6. */ + c_bit_length = 6; + while ((c >> c_bit_length) > 0U) { + ++c_bit_length; + } + + /* Initialise d and a. */ + d = c >> (c_bit_length - 5); + b = _PyLong_Rshift(n, 2U*c - 62U); + if (b == NULL) { + goto error; + } + m = (uint64_t)PyLong_AsUnsignedLongLong(b); + Py_DECREF(b); + if (m == (uint64_t)(-1) && PyErr_Occurred()) { + goto error; + } + u = _approximate_isqrt(m) >> (31U - d); + a = PyLong_FromUnsignedLongLong((unsigned long long)u); + if (a == NULL) { + goto error; + } + + for (int s = c_bit_length - 6; s >= 0; --s) { + PyObject *q; + size_t e = d; + + d = c >> s; + + /* q = (n >> 2*c - e - d + 1) // a */ + q = _PyLong_Rshift(n, 2U*c - d - e + 1U); + if (q == NULL) { + goto error; + } + Py_SETREF(q, PyNumber_FloorDivide(q, a)); + if (q == NULL) { + goto error; + } + + /* a = (a << d - 1 - e) + q */ + Py_SETREF(a, _PyLong_Lshift(a, d - 1U - e)); + if (a == NULL) { + Py_DECREF(q); + goto error; + } + Py_SETREF(a, PyNumber_Add(a, q)); + Py_DECREF(q); + if (a == NULL) { + goto error; + } + } + + /* The correct result is either a or a - 1. Figure out which, and + decrement a if necessary. */ + + /* a_too_large = n < a * a */ + b = PyNumber_Multiply(a, a); + if (b == NULL) { + goto error; + } + a_too_large = PyObject_RichCompareBool(n, b, Py_LT); + Py_DECREF(b); + if (a_too_large == -1) { + goto error; + } + + if (a_too_large) { + Py_SETREF(a, PyNumber_Subtract(a, _PyLong_One)); + } + Py_DECREF(n); + return a; + + error: + Py_XDECREF(a); + Py_DECREF(n); + return NULL; +} + +/* Divide-and-conquer factorial algorithm + * + * Based on the formula and pseudo-code provided at: + * http://www.luschny.de/math/factorial/binarysplitfact.html + * + * Faster algorithms exist, but they're more complicated and depend on + * a fast prime factorization algorithm. + * + * Notes on the algorithm + * ---------------------- + * + * factorial(n) is written in the form 2**k * m, with m odd. k and m are + * computed separately, and then combined using a left shift. + * + * The function factorial_odd_part computes the odd part m (i.e., the greatest + * odd divisor) of factorial(n), using the formula: + * + * factorial_odd_part(n) = + * + * product_{i >= 0} product_{0 < j <= n / 2**i, j odd} j + * + * Example: factorial_odd_part(20) = + * + * (1) * + * (1) * + * (1 * 3 * 5) * + * (1 * 3 * 5 * 7 * 9) + * (1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19) + * + * Here i goes from large to small: the first term corresponds to i=4 (any + * larger i gives an empty product), and the last term corresponds to i=0. + * Each term can be computed from the last by multiplying by the extra odd + * numbers required: e.g., to get from the penultimate term to the last one, + * we multiply by (11 * 13 * 15 * 17 * 19). + * + * To see a hint of why this formula works, here are the same numbers as above + * but with the even parts (i.e., the appropriate powers of 2) included. For + * each subterm in the product for i, we multiply that subterm by 2**i: + * + * factorial(20) = + * + * (16) * + * (8) * + * (4 * 12 * 20) * + * (2 * 6 * 10 * 14 * 18) * + * (1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19) + * + * The factorial_partial_product function computes the product of all odd j in + * range(start, stop) for given start and stop. It's used to compute the + * partial products like (11 * 13 * 15 * 17 * 19) in the example above. It + * operates recursively, repeatedly splitting the range into two roughly equal + * pieces until the subranges are small enough to be computed using only C + * integer arithmetic. + * + * The two-valuation k (i.e., the exponent of the largest power of 2 dividing + * the factorial) is computed independently in the main math_factorial + * function. By standard results, its value is: + * + * two_valuation = n//2 + n//4 + n//8 + .... + * + * It can be shown (e.g., by complete induction on n) that two_valuation is + * equal to n - count_set_bits(n), where count_set_bits(n) gives the number of + * '1'-bits in the binary expansion of n. + */ + +/* factorial_partial_product: Compute product(range(start, stop, 2)) using + * divide and conquer. Assumes start and stop are odd and stop > start. + * max_bits must be >= bit_length(stop - 2). */ + +static PyObject * +factorial_partial_product(unsigned long start, unsigned long stop, + unsigned long max_bits) +{ + unsigned long midpoint, num_operands; + PyObject *left = NULL, *right = NULL, *result = NULL; + + /* If the return value will fit an unsigned long, then we can + * multiply in a tight, fast loop where each multiply is O(1). + * Compute an upper bound on the number of bits required to store + * the answer. + * + * Storing some integer z requires floor(lg(z))+1 bits, which is + * conveniently the value returned by bit_length(z). The + * product x*y will require at most + * bit_length(x) + bit_length(y) bits to store, based + * on the idea that lg product = lg x + lg y. + * + * We know that stop - 2 is the largest number to be multiplied. From + * there, we have: bit_length(answer) <= num_operands * + * bit_length(stop - 2) + */ + + num_operands = (stop - start) / 2; + /* The "num_operands <= 8 * SIZEOF_LONG" check guards against the + * unlikely case of an overflow in num_operands * max_bits. */ + if (num_operands <= 8 * SIZEOF_LONG && + num_operands * max_bits <= 8 * SIZEOF_LONG) { + unsigned long j, total; + for (total = start, j = start + 2; j < stop; j += 2) + total *= j; + return PyLong_FromUnsignedLong(total); + } + + /* find midpoint of range(start, stop), rounded up to next odd number. */ + midpoint = (start + num_operands) | 1; + left = factorial_partial_product(start, midpoint, + _Py_bit_length(midpoint - 2)); + if (left == NULL) + goto error; + right = factorial_partial_product(midpoint, stop, max_bits); + if (right == NULL) + goto error; + result = PyNumber_Multiply(left, right); + + error: + Py_XDECREF(left); + Py_XDECREF(right); + return result; +} + +/* factorial_odd_part: compute the odd part of factorial(n). */ + +static PyObject * +factorial_odd_part(unsigned long n) +{ + long i; + unsigned long v, lower, upper; + PyObject *partial, *tmp, *inner, *outer; + + inner = PyLong_FromLong(1); + if (inner == NULL) + return NULL; + outer = inner; + Py_INCREF(outer); + + upper = 3; + for (i = _Py_bit_length(n) - 2; i >= 0; i--) { + v = n >> i; + if (v <= 2) + continue; + lower = upper; + /* (v + 1) | 1 = least odd integer strictly larger than n / 2**i */ + upper = (v + 1) | 1; + /* Here inner is the product of all odd integers j in the range (0, + n/2**(i+1)]. The factorial_partial_product call below gives the + product of all odd integers j in the range (n/2**(i+1), n/2**i]. */ + partial = factorial_partial_product(lower, upper, _Py_bit_length(upper-2)); + /* inner *= partial */ + if (partial == NULL) + goto error; + tmp = PyNumber_Multiply(inner, partial); + Py_DECREF(partial); + if (tmp == NULL) + goto error; + Py_DECREF(inner); + inner = tmp; + /* Now inner is the product of all odd integers j in the range (0, + n/2**i], giving the inner product in the formula above. */ + + /* outer *= inner; */ + tmp = PyNumber_Multiply(outer, inner); + if (tmp == NULL) + goto error; + Py_DECREF(outer); + outer = tmp; + } + Py_DECREF(inner); + return outer; + + error: + Py_DECREF(outer); + Py_DECREF(inner); + return NULL; +} + + +/* Lookup table for small factorial values */ + +static const unsigned long SmallFactorials[] = { + 1, 1, 2, 6, 24, 120, 720, 5040, 40320, + 362880, 3628800, 39916800, 479001600, +#if SIZEOF_LONG >= 8 + 6227020800, 87178291200, 1307674368000, + 20922789888000, 355687428096000, 6402373705728000, + 121645100408832000, 2432902008176640000 +#endif +}; + +/*[clinic input] +math.factorial + + x as arg: object + / + +Find x!. + +Raise a ValueError if x is negative or non-integral. +[clinic start generated code]*/ + +static PyObject * +math_factorial(PyObject *module, PyObject *arg) +/*[clinic end generated code: output=6686f26fae00e9ca input=6d1c8105c0d91fb4]*/ +{ + long x, two_valuation; + int overflow; + PyObject *result, *odd_part, *pyint_form; + + if (PyFloat_Check(arg)) { + if (PyErr_WarnEx(PyExc_DeprecationWarning, + "Using factorial() with floats is deprecated", + 1) < 0) + { + return NULL; + } + PyObject *lx; + double dx = PyFloat_AS_DOUBLE((PyFloatObject *)arg); + if (!(Py_IS_FINITE(dx) && dx == floor(dx))) { + PyErr_SetString(PyExc_ValueError, + "factorial() only accepts integral values"); + return NULL; + } + lx = PyLong_FromDouble(dx); + if (lx == NULL) + return NULL; + x = PyLong_AsLongAndOverflow(lx, &overflow); + Py_DECREF(lx); + } + else { + pyint_form = PyNumber_Index(arg); + if (pyint_form == NULL) { + return NULL; + } + x = PyLong_AsLongAndOverflow(pyint_form, &overflow); + Py_DECREF(pyint_form); + } + + if (x == -1 && PyErr_Occurred()) { + return NULL; + } + else if (overflow == 1) { + PyErr_Format(PyExc_OverflowError, + "factorial() argument should not exceed %ld", + LONG_MAX); + return NULL; + } + else if (overflow == -1 || x < 0) { + PyErr_SetString(PyExc_ValueError, + "factorial() not defined for negative values"); + return NULL; + } + + /* use lookup table if x is small */ + if (x < (long)Py_ARRAY_LENGTH(SmallFactorials)) + return PyLong_FromUnsignedLong(SmallFactorials[x]); + + /* else express in the form odd_part * 2**two_valuation, and compute as + odd_part << two_valuation. */ + odd_part = factorial_odd_part(x); + if (odd_part == NULL) + return NULL; + two_valuation = x - count_set_bits(x); + result = _PyLong_Lshift(odd_part, two_valuation); + Py_DECREF(odd_part); + return result; +} + + +/*[clinic input] +math.lcm + x as a: object + y as b: object + / +least common multiple of x and y +[clinic start generated code]*/ + +static PyObject * +math_lcm_impl(PyObject *module, PyObject *a, PyObject *b) +/*[clinic end generated code: output=6f83fb6d671074ba input=efb3d7b7334b7118]*/ + + PyObject *g, *m, *f, *ab; + + a = PyNumber_Index(a); + if (a == NULL) { + return NULL; + } + b = PyNumber_Index(b); + if (b == NULL) { + Py_DECREF(a); + return NULL; + } + if (_PyLong_Sign(a) == 0 || _PyLong_Sign(b) == 0) { + Py_DECREF(a); + Py_DECREF(b); + return PyLong_FromLong(0); + } + g = _PyLong_GCD(a, b); + if (g == NULL) { + Py_DECREF(a); + Py_DECREF(b); + return NULL; + } + f = PyNumber_FloorDivide(a, g); + Py_DECREF(g); + Py_DECREF(a); + if (f == NULL) { + Py_DECREF(b); + return NULL; + } + m = PyNumber_Multiply(f, b); + Py_DECREF(f); + Py_DECREF(b); + if (m == NULL) { + return NULL; + } + ab = PyNumber_Absolute(m); + Py_DECREF(m); + return ab; +} + + +/*[clinic input] +math.trunc + + x: object + / + +Truncates the Real x to the nearest Integral toward 0. + +Uses the __trunc__ magic method. +[clinic start generated code]*/ + +static PyObject * +math_trunc(PyObject *module, PyObject *x) +/*[clinic end generated code: output=34b9697b707e1031 input=2168b34e0a09134d]*/ +{ + _Py_IDENTIFIER(__trunc__); + PyObject *trunc, *result; + + if (PyFloat_CheckExact(x)) { + return PyFloat_Type.tp_as_number->nb_int(x); + } + + if (Py_TYPE(x)->tp_dict == NULL) { + if (PyType_Ready(Py_TYPE(x)) < 0) + return NULL; + } + + trunc = _PyObject_LookupSpecial(x, &PyId___trunc__); + if (trunc == NULL) { + if (!PyErr_Occurred()) + PyErr_Format(PyExc_TypeError, + "type %.100s doesn't define __trunc__ method", + Py_TYPE(x)->tp_name); + return NULL; + } + result = _PyObject_CallNoArg(trunc); + Py_DECREF(trunc); + return result; +} + + +/*[clinic input] +math.frexp + + x: double + / + +Return the mantissa and exponent of x, as pair (m, e). + +m is a float and e is an int, such that x = m * 2.**e. +If x is 0, m and e are both 0. Else 0.5 <= abs(m) < 1.0. +[clinic start generated code]*/ + +static PyObject * +math_frexp_impl(PyObject *module, double x) +/*[clinic end generated code: output=03e30d252a15ad4a input=96251c9e208bc6e9]*/ +{ + int i; + /* deal with special cases directly, to sidestep platform + differences */ + if (Py_IS_NAN(x) || Py_IS_INFINITY(x) || !x) { + i = 0; + } + else { + x = frexp(x, &i); + } + return Py_BuildValue("(di)", x, i); +} + + +/*[clinic input] +math.ldexp + + x: double + i: object + / + +Return x * (2**i). + +This is essentially the inverse of frexp(). +[clinic start generated code]*/ + +static PyObject * +math_ldexp_impl(PyObject *module, double x, PyObject *i) +/*[clinic end generated code: output=b6892f3c2df9cc6a input=17d5970c1a40a8c1]*/ +{ + double r; + long exp; + int overflow; + + if (PyLong_Check(i)) { + /* on overflow, replace exponent with either LONG_MAX + or LONG_MIN, depending on the sign. */ + exp = PyLong_AsLongAndOverflow(i, &overflow); + if (exp == -1 && PyErr_Occurred()) + return NULL; + if (overflow) + exp = overflow < 0 ? LONG_MIN : LONG_MAX; + } + else { + PyErr_SetString(PyExc_TypeError, + "Expected an int as second argument to ldexp."); + return NULL; + } + + if (x == 0. || !Py_IS_FINITE(x)) { + /* NaNs, zeros and infinities are returned unchanged */ + r = x; + errno = 0; + } else if (exp > INT_MAX) { + /* overflow */ + r = copysign(Py_HUGE_VAL, x); + errno = ERANGE; + } else if (exp < INT_MIN) { + /* underflow to +-0 */ + r = copysign(0., x); + errno = 0; + } else { + errno = 0; + r = ldexp(x, (int)exp); + if (Py_IS_INFINITY(r)) + errno = ERANGE; + } + + if (errno && is_error(r)) + return NULL; + return PyFloat_FromDouble(r); +} + + +/*[clinic input] +math.modf + + x: double + / + +Return the fractional and integer parts of x. + +Both results carry the sign of x and are floats. +[clinic start generated code]*/ + +static PyObject * +math_modf_impl(PyObject *module, double x) +/*[clinic end generated code: output=90cee0260014c3c0 input=b4cfb6786afd9035]*/ +{ + double y; + /* some platforms don't do the right thing for NaNs and + infinities, so we take care of special cases directly. */ + if (!Py_IS_FINITE(x)) { + if (Py_IS_INFINITY(x)) + return Py_BuildValue("(dd)", copysign(0., x), x); + else if (Py_IS_NAN(x)) + return Py_BuildValue("(dd)", x, x); + } + + errno = 0; + x = modf(x, &y); + return Py_BuildValue("(dd)", x, y); +} + + +/* A decent logarithm is easy to compute even for huge ints, but libm can't + do that by itself -- loghelper can. func is log or log10, and name is + "log" or "log10". Note that overflow of the result isn't possible: an int + can contain no more than INT_MAX * SHIFT bits, so has value certainly less + than 2**(2**64 * 2**16) == 2**2**80, and log2 of that is 2**80, which is + small enough to fit in an IEEE single. log and log10 are even smaller. + However, intermediate overflow is possible for an int if the number of bits + in that int is larger than PY_SSIZE_T_MAX. */ + +static PyObject* +loghelper(PyObject* arg, double (*func)(double), const char *funcname) +{ + /* If it is int, do it ourselves. */ + if (PyLong_Check(arg)) { + double x, result; + Py_ssize_t e; + + /* Negative or zero inputs give a ValueError. */ + if (Py_SIZE(arg) <= 0) { + PyErr_SetString(PyExc_ValueError, + "math domain error"); + return NULL; + } + + x = PyLong_AsDouble(arg); + if (x == -1.0 && PyErr_Occurred()) { + if (!PyErr_ExceptionMatches(PyExc_OverflowError)) + return NULL; + /* Here the conversion to double overflowed, but it's possible + to compute the log anyway. Clear the exception and continue. */ + PyErr_Clear(); + x = _PyLong_Frexp((PyLongObject *)arg, &e); + if (x == -1.0 && PyErr_Occurred()) + return NULL; + /* Value is ~= x * 2**e, so the log ~= log(x) + log(2) * e. */ + result = func(x) + func(2.0) * e; + } + else + /* Successfully converted x to a double. */ + result = func(x); + return PyFloat_FromDouble(result); + } + + /* Else let libm handle it by itself. */ + return math_1(arg, func, 0); +} + + +/*[clinic input] +math.log + + x: object + [ + base: object(c_default="NULL") = math.e + ] + / + +Return the logarithm of x to the given base. + +If the base not specified, returns the natural logarithm (base e) of x. +[clinic start generated code]*/ + +static PyObject * +math_log_impl(PyObject *module, PyObject *x, int group_right_1, + PyObject *base) +/*[clinic end generated code: output=7b5a39e526b73fc9 input=0f62d5726cbfebbd]*/ +{ + PyObject *num, *den; + PyObject *ans; + + num = loghelper(x, m_log, "log"); + if (num == NULL || base == NULL) + return num; + + den = loghelper(base, m_log, "log"); + if (den == NULL) { + Py_DECREF(num); + return NULL; + } + + ans = PyNumber_TrueDivide(num, den); + Py_DECREF(num); + Py_DECREF(den); + return ans; +} + + +/*[clinic input] +math.log2 + + x: object + / + +Return the base 2 logarithm of x. +[clinic start generated code]*/ + +static PyObject * +math_log2(PyObject *module, PyObject *x) +/*[clinic end generated code: output=5425899a4d5d6acb input=08321262bae4f39b]*/ +{ + return loghelper(x, m_log2, "log2"); +} + + +/*[clinic input] +math.log10 + + x: object + / + +Return the base 10 logarithm of x. +[clinic start generated code]*/ + +static PyObject * +math_log10(PyObject *module, PyObject *x) +/*[clinic end generated code: output=be72a64617df9c6f input=b2469d02c6469e53]*/ +{ + return loghelper(x, m_log10, "log10"); +} + + +/*[clinic input] +math.fmod + + x: double + y: double + / + +Return fmod(x, y), according to platform C. + +x % y may differ. +[clinic start generated code]*/ + +static PyObject * +math_fmod_impl(PyObject *module, double x, double y) +/*[clinic end generated code: output=7559d794343a27b5 input=4f84caa8cfc26a03]*/ +{ + double r; + /* fmod(x, +/-Inf) returns x for finite x. */ + if (Py_IS_INFINITY(y) && Py_IS_FINITE(x)) + return PyFloat_FromDouble(x); + errno = 0; + r = fmod(x, y); + if (Py_IS_NAN(r)) { + if (!Py_IS_NAN(x) && !Py_IS_NAN(y)) + errno = EDOM; + else + errno = 0; + } + if (errno && is_error(r)) + return NULL; + else + return PyFloat_FromDouble(r); +} + +/* +Given an *n* length *vec* of values and a value *max*, compute: + + max * sqrt(sum((x / max) ** 2 for x in vec)) + +The value of the *max* variable must be non-negative and +equal to the absolute value of the largest magnitude +entry in the vector. If n==0, then *max* should be 0.0. +If an infinity is present in the vec, *max* should be INF. + +The *found_nan* variable indicates whether some member of +the *vec* is a NaN. + +To improve accuracy and to increase the number of cases where +vector_norm() is commutative, we use a variant of Neumaier +summation specialized to exploit that we always know that +|csum| >= |x|. + +The *csum* variable tracks the cumulative sum and *frac* tracks +the cumulative fractional errors at each step. Since this +variant assumes that |csum| >= |x| at each step, we establish +the precondition by starting the accumulation from 1.0 which +represents the largest possible value of (x/max)**2. + +After the loop is finished, the initial 1.0 is subtracted out +for a net zero effect on the final sum. Since *csum* will be +greater than 1.0, the subtraction of 1.0 will not cause +fractional digits to be dropped from *csum*. + +*/ + +static inline double +vector_norm(Py_ssize_t n, double *vec, double max, int found_nan) +{ + double x, csum = 1.0, oldcsum, frac = 0.0; + Py_ssize_t i; + + if (Py_IS_INFINITY(max)) { + return max; + } + if (found_nan) { + return Py_NAN; + } + if (max == 0.0 || n <= 1) { + return max; + } + for (i=0 ; i < n ; i++) { + x = vec[i]; + assert(Py_IS_FINITE(x) && fabs(x) <= max); + x /= max; + x = x*x; + oldcsum = csum; + csum += x; + assert(csum >= x); + frac += (oldcsum - csum) + x; + } + return max * sqrt(csum - 1.0 + frac); +} + +#define NUM_STACK_ELEMS 16 + +/*[clinic input] +math.dist + + p: object + q: object + / + +Return the Euclidean distance between two points p and q. + +The points should be specified as sequences (or iterables) of +coordinates. Both inputs must have the same dimension. + +Roughly equivalent to: + sqrt(sum((px - qx) ** 2.0 for px, qx in zip(p, q))) +[clinic start generated code]*/ + +static PyObject * +math_dist_impl(PyObject *module, PyObject *p, PyObject *q) +/*[clinic end generated code: output=56bd9538d06bbcfe input=74e85e1b6092e68e]*/ +{ + PyObject *item; + double max = 0.0; + double x, px, qx, result; + Py_ssize_t i, m, n; + int found_nan = 0, p_allocated = 0, q_allocated = 0; + double diffs_on_stack[NUM_STACK_ELEMS]; + double *diffs = diffs_on_stack; + + if (!PyTuple_Check(p)) { + p = PySequence_Tuple(p); + if (p == NULL) { + return NULL; + } + p_allocated = 1; + } + if (!PyTuple_Check(q)) { + q = PySequence_Tuple(q); + if (q == NULL) { + if (p_allocated) { + Py_DECREF(p); + } + return NULL; + } + q_allocated = 1; + } + + m = PyTuple_GET_SIZE(p); + n = PyTuple_GET_SIZE(q); + if (m != n) { + PyErr_SetString(PyExc_ValueError, + "both points must have the same number of dimensions"); + return NULL; + + } + if (n > NUM_STACK_ELEMS) { + diffs = (double *) PyObject_Malloc(n * sizeof(double)); + if (diffs == NULL) { + return PyErr_NoMemory(); + } + } + for (i=0 ; i max) { + max = x; + } + } + result = vector_norm(n, diffs, max, found_nan); + if (diffs != diffs_on_stack) { + PyObject_Free(diffs); + } + if (p_allocated) { + Py_DECREF(p); + } + if (q_allocated) { + Py_DECREF(q); + } + return PyFloat_FromDouble(result); + + error_exit: + if (diffs != diffs_on_stack) { + PyObject_Free(diffs); + } + if (p_allocated) { + Py_DECREF(p); + } + if (q_allocated) { + Py_DECREF(q); + } + return NULL; +} + +/* AC: cannot convert yet, waiting for *args support */ +static PyObject * +math_hypot(PyObject *self, PyObject *const *args, Py_ssize_t nargs) +{ + Py_ssize_t i; + PyObject *item; + double max = 0.0; + double x, result; + int found_nan = 0; + double coord_on_stack[NUM_STACK_ELEMS]; + double *coordinates = coord_on_stack; + + if (nargs > NUM_STACK_ELEMS) { + coordinates = (double *) PyObject_Malloc(nargs * sizeof(double)); + if (coordinates == NULL) { + return PyErr_NoMemory(); + } + } + for (i = 0; i < nargs; i++) { + item = args[i]; + ASSIGN_DOUBLE(x, item, error_exit); + x = fabs(x); + coordinates[i] = x; + found_nan |= Py_IS_NAN(x); + if (x > max) { + max = x; + } + } + result = vector_norm(nargs, coordinates, max, found_nan); + if (coordinates != coord_on_stack) { + PyObject_Free(coordinates); + } + return PyFloat_FromDouble(result); + + error_exit: + if (coordinates != coord_on_stack) { + PyObject_Free(coordinates); + } + return NULL; +} + +#undef NUM_STACK_ELEMS + +PyDoc_STRVAR(math_hypot_doc, + "hypot(*coordinates) -> value\n\n\ +Multidimensional Euclidean distance from the origin to a point.\n\ +\n\ +Roughly equivalent to:\n\ + sqrt(sum(x**2 for x in coordinates))\n\ +\n\ +For a two dimensional point (x, y), gives the hypotenuse\n\ +using the Pythagorean theorem: sqrt(x*x + y*y).\n\ +\n\ +For example, the hypotenuse of a 3/4/5 right triangle is:\n\ +\n\ + >>> hypot(3.0, 4.0)\n\ + 5.0\n\ +"); + +/* pow can't use math_2, but needs its own wrapper: the problem is + that an infinite result can arise either as a result of overflow + (in which case OverflowError should be raised) or as a result of + e.g. 0.**-5. (for which ValueError needs to be raised.) +*/ + +/*[clinic input] +math.pow + + x: double + y: double + / + +Return x**y (x to the power of y). +[clinic start generated code]*/ + +static PyObject * +math_pow_impl(PyObject *module, double x, double y) +/*[clinic end generated code: output=fff93e65abccd6b0 input=c26f1f6075088bfd]*/ +{ + double r; + int odd_y; + + /* deal directly with IEEE specials, to cope with problems on various + platforms whose semantics don't exactly match C99 */ + r = 0.; /* silence compiler warning */ + if (!Py_IS_FINITE(x) || !Py_IS_FINITE(y)) { + errno = 0; + if (Py_IS_NAN(x)) + r = y == 0. ? 1. : x; /* NaN**0 = 1 */ + else if (Py_IS_NAN(y)) + r = x == 1. ? 1. : y; /* 1**NaN = 1 */ + else if (Py_IS_INFINITY(x)) { + odd_y = Py_IS_FINITE(y) && fmod(fabs(y), 2.0) == 1.0; + if (y > 0.) + r = odd_y ? x : fabs(x); + else if (y == 0.) + r = 1.; + else /* y < 0. */ + r = odd_y ? copysign(0., x) : 0.; + } + else if (Py_IS_INFINITY(y)) { + if (fabs(x) == 1.0) + r = 1.; + else if (y > 0. && fabs(x) > 1.0) + r = y; + else if (y < 0. && fabs(x) < 1.0) { + r = -y; /* result is +inf */ + if (x == 0.) /* 0**-inf: divide-by-zero */ + errno = EDOM; + } + else + r = 0.; + } + } + else { + /* let libm handle finite**finite */ + errno = 0; + r = pow(x, y); + /* a NaN result should arise only from (-ve)**(finite + non-integer); in this case we want to raise ValueError. */ + if (!Py_IS_FINITE(r)) { + if (Py_IS_NAN(r)) { + errno = EDOM; + } + /* + an infinite result here arises either from: + (A) (+/-0.)**negative (-> divide-by-zero) + (B) overflow of x**y with x and y finite + */ + else if (Py_IS_INFINITY(r)) { + if (x == 0.) + errno = EDOM; + else + errno = ERANGE; + } + } + } + + if (errno && is_error(r)) + return NULL; else - { - x = PyFloat_AsDouble(arg); - if (x == -1.0 && PyErr_Occurred()) { - goto exit; + return PyFloat_FromDouble(r); +} + + +static const double degToRad = Py_MATH_PI / 180.0; +static const double radToDeg = 180.0 / Py_MATH_PI; + +/*[clinic input] +math.degrees + + x: double + / + +Convert angle x from radians to degrees. +[clinic start generated code]*/ + +static PyObject * +math_degrees_impl(PyObject *module, double x) +/*[clinic end generated code: output=7fea78b294acd12f input=81e016555d6e3660]*/ +{ + return PyFloat_FromDouble(x * radToDeg); +} + + +/*[clinic input] +math.radians + + x: double + / + +Convert angle x from degrees to radians. +[clinic start generated code]*/ + +static PyObject * +math_radians_impl(PyObject *module, double x) +/*[clinic end generated code: output=34daa47caf9b1590 input=91626fc489fe3d63]*/ +{ + return PyFloat_FromDouble(x * degToRad); +} + + +/*[clinic input] +math.isfinite + + x: double + / + +Return True if x is neither an infinity nor a NaN, and False otherwise. +[clinic start generated code]*/ + +static PyObject * +math_isfinite_impl(PyObject *module, double x) +/*[clinic end generated code: output=8ba1f396440c9901 input=46967d254812e54a]*/ +{ + return PyBool_FromLong((long)Py_IS_FINITE(x)); +} + + +/*[clinic input] +math.isnan + + x: double + / + +Return True if x is a NaN (not a number), and False otherwise. +[clinic start generated code]*/ + +static PyObject * +math_isnan_impl(PyObject *module, double x) +/*[clinic end generated code: output=f537b4d6df878c3e input=935891e66083f46a]*/ +{ + return PyBool_FromLong((long)Py_IS_NAN(x)); +} + + +/*[clinic input] +math.isinf + + x: double + / + +Return True if x is a positive or negative infinity, and False otherwise. +[clinic start generated code]*/ + +static PyObject * +math_isinf_impl(PyObject *module, double x) +/*[clinic end generated code: output=9f00cbec4de7b06b input=32630e4212cf961f]*/ +{ + return PyBool_FromLong((long)Py_IS_INFINITY(x)); +} + + +/*[clinic input] +math.isclose -> bool + + a: double + b: double + * + rel_tol: double = 1e-09 + maximum difference for being considered "close", relative to the + magnitude of the input values + abs_tol: double = 0.0 + maximum difference for being considered "close", regardless of the + magnitude of the input values + +Determine whether two floating point numbers are close in value. + +Return True if a is close in value to b, and False otherwise. + +For the values to be considered close, the difference between them +must be smaller than at least one of the tolerances. + +-inf, inf and NaN behave similarly to the IEEE 754 Standard. That +is, NaN is not close to anything, even itself. inf and -inf are +only close to themselves. +[clinic start generated code]*/ + +static int +math_isclose_impl(PyObject *module, double a, double b, double rel_tol, + double abs_tol) +/*[clinic end generated code: output=b73070207511952d input=f28671871ea5bfba]*/ +{ + double diff = 0.0; + + /* sanity check on the inputs */ + if (rel_tol < 0.0 || abs_tol < 0.0 ) { + PyErr_SetString(PyExc_ValueError, + "tolerances must be non-negative"); + return -1; + } + + if ( a == b ) { + /* short circuit exact equality -- needed to catch two infinities of + the same sign. And perhaps speeds things up a bit sometimes. + */ + return 1; + } + + /* This catches the case of two infinities of opposite sign, or + one infinity and one finite number. Two infinities of opposite + sign would otherwise have an infinite relative tolerance. + Two infinities of the same sign are caught by the equality check + above. + */ + + if (Py_IS_INFINITY(a) || Py_IS_INFINITY(b)) { + return 0; + } + + /* now do the regular computation + this is essentially the "weak" test from the Boost library + */ + + diff = fabs(b - a); + + return (((diff <= fabs(rel_tol * b)) || + (diff <= fabs(rel_tol * a))) || + (diff <= abs_tol)); +} + +static inline int +_check_long_mult_overflow(long a, long b) { + + /* From Python2's int_mul code: + + Integer overflow checking for * is painful: Python tried a couple ways, but + they didn't work on all platforms, or failed in endcases (a product of + -sys.maxint-1 has been a particular pain). + + Here's another way: + + The native long product x*y is either exactly right or *way* off, being + just the last n bits of the true product, where n is the number of bits + in a long (the delivered product is the true product plus i*2**n for + some integer i). + + The native double product (double)x * (double)y is subject to three + rounding errors: on a sizeof(long)==8 box, each cast to double can lose + info, and even on a sizeof(long)==4 box, the multiplication can lose info. + But, unlike the native long product, it's not in *range* trouble: even + if sizeof(long)==32 (256-bit longs), the product easily fits in the + dynamic range of a double. So the leading 50 (or so) bits of the double + product are correct. + + We check these two ways against each other, and declare victory if they're + approximately the same. Else, because the native long product is the only + one that can lose catastrophic amounts of information, it's the native long + product that must have overflowed. + + */ + + long longprod = (long)((unsigned long)a * b); + double doubleprod = (double)a * (double)b; + double doubled_longprod = (double)longprod; + + if (doubled_longprod == doubleprod) { + return 0; + } + + const double diff = doubled_longprod - doubleprod; + const double absdiff = diff >= 0.0 ? diff : -diff; + const double absprod = doubleprod >= 0.0 ? doubleprod : -doubleprod; + + if (32.0 * absdiff <= absprod) { + return 0; + } + + return 1; +} + +/*[clinic input] +math.prod + + iterable: object + / + * + start: object(c_default="NULL") = 1 + +Calculate the product of all the elements in the input iterable. + +The default start value for the product is 1. + +When the iterable is empty, return the start value. This function is +intended specifically for use with numeric values and may reject +non-numeric types. +[clinic start generated code]*/ + +static PyObject * +math_prod_impl(PyObject *module, PyObject *iterable, PyObject *start) +/*[clinic end generated code: output=36153bedac74a198 input=4c5ab0682782ed54]*/ +{ + PyObject *result = start; + PyObject *temp, *item, *iter; + + iter = PyObject_GetIter(iterable); + if (iter == NULL) { + return NULL; + } + + if (result == NULL) { + result = PyLong_FromLong(1); + if (result == NULL) { + Py_DECREF(iter); + return NULL; } + } else { + Py_INCREF(result); } - _return_value = math_ulp_impl(module, x); - if ((_return_value == -1.0) && PyErr_Occurred()) { - goto exit; +#ifndef SLOW_PROD + /* Fast paths for integers keeping temporary products in C. + * Assumes all inputs are the same type. + * If the assumption fails, default to use PyObjects instead. + */ + if (PyLong_CheckExact(result)) { + int overflow; + long i_result = PyLong_AsLongAndOverflow(result, &overflow); + /* If this already overflowed, don't even enter the loop. */ + if (overflow == 0) { + Py_DECREF(result); + result = NULL; + } + /* Loop over all the items in the iterable until we finish, we overflow + * or we found a non integer element */ + while(result == NULL) { + item = PyIter_Next(iter); + if (item == NULL) { + Py_DECREF(iter); + if (PyErr_Occurred()) { + return NULL; + } + return PyLong_FromLong(i_result); + } + if (PyLong_CheckExact(item)) { + long b = PyLong_AsLongAndOverflow(item, &overflow); + if (overflow == 0 && !_check_long_mult_overflow(i_result, b)) { + long x = i_result * b; + i_result = x; + Py_DECREF(item); + continue; + } + } + /* Either overflowed or is not an int. + * Restore real objects and process normally */ + result = PyLong_FromLong(i_result); + if (result == NULL) { + Py_DECREF(item); + Py_DECREF(iter); + return NULL; + } + temp = PyNumber_Multiply(result, item); + Py_DECREF(result); + Py_DECREF(item); + result = temp; + if (result == NULL) { + Py_DECREF(iter); + return NULL; + } + } + } + + /* Fast paths for floats keeping temporary products in C. + * Assumes all inputs are the same type. + * If the assumption fails, default to use PyObjects instead. + */ + if (PyFloat_CheckExact(result)) { + double f_result = PyFloat_AS_DOUBLE(result); + Py_DECREF(result); + result = NULL; + while(result == NULL) { + item = PyIter_Next(iter); + if (item == NULL) { + Py_DECREF(iter); + if (PyErr_Occurred()) { + return NULL; + } + return PyFloat_FromDouble(f_result); + } + if (PyFloat_CheckExact(item)) { + f_result *= PyFloat_AS_DOUBLE(item); + Py_DECREF(item); + continue; + } + if (PyLong_CheckExact(item)) { + long value; + int overflow; + value = PyLong_AsLongAndOverflow(item, &overflow); + if (!overflow) { + f_result *= (double)value; + Py_DECREF(item); + continue; + } + } + result = PyFloat_FromDouble(f_result); + if (result == NULL) { + Py_DECREF(item); + Py_DECREF(iter); + return NULL; + } + temp = PyNumber_Multiply(result, item); + Py_DECREF(result); + Py_DECREF(item); + result = temp; + if (result == NULL) { + Py_DECREF(iter); + return NULL; + } + } + } +#endif + /* Consume rest of the iterable (if any) that could not be handled + * by specialized functions above.*/ + for(;;) { + item = PyIter_Next(iter); + if (item == NULL) { + /* error, or end-of-sequence */ + if (PyErr_Occurred()) { + Py_DECREF(result); + result = NULL; + } + break; + } + temp = PyNumber_Multiply(result, item); + Py_DECREF(result); + Py_DECREF(item); + result = temp; + if (result == NULL) + break; + } + Py_DECREF(iter); + return result; +} + + +/*[clinic input] +math.perm + + n: object + k: object = None + / + +Number of ways to choose k items from n items without repetition and with order. + +Evaluates to n! / (n - k)! when k <= n and evaluates +to zero when k > n. + +If k is not specified or is None, then k defaults to n +and the function returns n!. + +Raises TypeError if either of the arguments are not integers. +Raises ValueError if either of the arguments are negative. +[clinic start generated code]*/ + +static PyObject * +math_perm_impl(PyObject *module, PyObject *n, PyObject *k) +/*[clinic end generated code: output=e021a25469653e23 input=5311c5a00f359b53]*/ +{ + PyObject *result = NULL, *factor = NULL; + int overflow, cmp; + long long i, factors; + + if (k == Py_None) { + return math_factorial(module, n); + } + n = PyNumber_Index(n); + if (n == NULL) { + return NULL; + } + if (!PyLong_CheckExact(n)) { + Py_SETREF(n, _PyLong_Copy((PyLongObject *)n)); + if (n == NULL) { + return NULL; + } + } + k = PyNumber_Index(k); + if (k == NULL) { + Py_DECREF(n); + return NULL; + } + if (!PyLong_CheckExact(k)) { + Py_SETREF(k, _PyLong_Copy((PyLongObject *)k)); + if (k == NULL) { + Py_DECREF(n); + return NULL; + } + } + + if (Py_SIZE(n) < 0) { + PyErr_SetString(PyExc_ValueError, + "n must be a non-negative integer"); + goto error; + } + if (Py_SIZE(k) < 0) { + PyErr_SetString(PyExc_ValueError, + "k must be a non-negative integer"); + goto error; + } + + cmp = PyObject_RichCompareBool(n, k, Py_LT); + if (cmp != 0) { + if (cmp > 0) { + result = PyLong_FromLong(0); + goto done; + } + goto error; + } + + factors = PyLong_AsLongLongAndOverflow(k, &overflow); + if (overflow > 0) { + PyErr_Format(PyExc_OverflowError, + "k must not exceed %lld", + LLONG_MAX); + goto error; + } + else if (factors == -1) { + /* k is nonnegative, so a return value of -1 can only indicate error */ + goto error; + } + + if (factors == 0) { + result = PyLong_FromLong(1); + goto done; + } + + result = n; + Py_INCREF(result); + if (factors == 1) { + goto done; + } + + factor = n; + Py_INCREF(factor); + for (i = 1; i < factors; ++i) { + Py_SETREF(factor, PyNumber_Subtract(factor, _PyLong_One)); + if (factor == NULL) { + goto error; + } + Py_SETREF(result, PyNumber_Multiply(result, factor)); + if (result == NULL) { + goto error; + } + } + Py_DECREF(factor); + +done: + Py_DECREF(n); + Py_DECREF(k); + return result; + +error: + Py_XDECREF(factor); + Py_XDECREF(result); + Py_DECREF(n); + Py_DECREF(k); + return NULL; +} + + +/*[clinic input] +math.comb + + n: object + k: object + / + +Number of ways to choose k items from n items without repetition and without order. + +Evaluates to n! / (k! * (n - k)!) when k <= n and evaluates +to zero when k > n. + +Also called the binomial coefficient because it is equivalent +to the coefficient of k-th term in polynomial expansion of the +expression (1 + x)**n. + +Raises TypeError if either of the arguments are not integers. +Raises ValueError if either of the arguments are negative. + +[clinic start generated code]*/ + +static PyObject * +math_comb_impl(PyObject *module, PyObject *n, PyObject *k) +/*[clinic end generated code: output=bd2cec8d854f3493 input=9a05315af2518709]*/ +{ + PyObject *result = NULL, *factor = NULL, *temp; + int overflow, cmp; + long long i, factors; + + n = PyNumber_Index(n); + if (n == NULL) { + return NULL; + } + if (!PyLong_CheckExact(n)) { + Py_SETREF(n, _PyLong_Copy((PyLongObject *)n)); + if (n == NULL) { + return NULL; + } + } + k = PyNumber_Index(k); + if (k == NULL) { + Py_DECREF(n); + return NULL; + } + if (!PyLong_CheckExact(k)) { + Py_SETREF(k, _PyLong_Copy((PyLongObject *)k)); + if (k == NULL) { + Py_DECREF(n); + return NULL; + } + } + + if (Py_SIZE(n) < 0) { + PyErr_SetString(PyExc_ValueError, + "n must be a non-negative integer"); + goto error; + } + if (Py_SIZE(k) < 0) { + PyErr_SetString(PyExc_ValueError, + "k must be a non-negative integer"); + goto error; + } + + /* k = min(k, n - k) */ + temp = PyNumber_Subtract(n, k); + if (temp == NULL) { + goto error; + } + if (Py_SIZE(temp) < 0) { + Py_DECREF(temp); + result = PyLong_FromLong(0); + goto done; + } + cmp = PyObject_RichCompareBool(temp, k, Py_LT); + if (cmp > 0) { + Py_SETREF(k, temp); + } + else { + Py_DECREF(temp); + if (cmp < 0) { + goto error; + } + } + + factors = PyLong_AsLongLongAndOverflow(k, &overflow); + if (overflow > 0) { + PyErr_Format(PyExc_OverflowError, + "min(n - k, k) must not exceed %lld", + LLONG_MAX); + goto error; + } + if (factors == -1) { + /* k is nonnegative, so a return value of -1 can only indicate error */ + goto error; + } + + if (factors == 0) { + result = PyLong_FromLong(1); + goto done; + } + + result = n; + Py_INCREF(result); + if (factors == 1) { + goto done; + } + + factor = n; + Py_INCREF(factor); + for (i = 1; i < factors; ++i) { + Py_SETREF(factor, PyNumber_Subtract(factor, _PyLong_One)); + if (factor == NULL) { + goto error; + } + Py_SETREF(result, PyNumber_Multiply(result, factor)); + if (result == NULL) { + goto error; + } + + temp = PyLong_FromUnsignedLongLong((unsigned long long)i + 1); + if (temp == NULL) { + goto error; + } + Py_SETREF(result, PyNumber_FloorDivide(result, temp)); + Py_DECREF(temp); + if (result == NULL) { + goto error; + } + } + Py_DECREF(factor); + +done: + Py_DECREF(n); + Py_DECREF(k); + return result; + +error: + Py_XDECREF(factor); + Py_XDECREF(result); + Py_DECREF(n); + Py_DECREF(k); + return NULL; +} + + +/*[clinic input] +math.nextafter + + x: double + y: double + / + +Return the next floating-point value after x towards y. +[clinic start generated code]*/ + +static PyObject * +math_nextafter_impl(PyObject *module, double x, double y) +/*[clinic end generated code: output=750c8266c1c540ce input=02b2d50cd1d9f9b6]*/ +{ +#if defined(_AIX) + if (x == y) { + /* On AIX 7.1, libm nextafter(-0.0, +0.0) returns -0.0. + Bug fixed in bos.adt.libm 7.2.2.0 by APAR IV95512. */ + return PyFloat_FromDouble(y); } - return_value = PyFloat_FromDouble(_return_value); +#endif + return PyFloat_FromDouble(nextafter(x, y)); +} + + +/*[clinic input] +math.ulp -> double + + x: double + / + +Return the value of the least significant bit of the float x. +[clinic start generated code]*/ + +static double +math_ulp_impl(PyObject *module, double x) +/*[clinic end generated code: output=f5207867a9384dd4 input=31f9bfbbe373fcaa]*/ +{ + if (Py_IS_NAN(x)) { + return x; + } + x = fabs(x); + if (Py_IS_INFINITY(x)) { + return x; + } + double inf = m_inf(); + double x2 = nextafter(x, inf); + if (Py_IS_INFINITY(x2)) { + /* special case: x is the largest positive representable float */ + x2 = nextafter(x, -inf); + return x - x2; + } + return x2 - x; +} -exit: - return return_value; + +static PyMethodDef math_methods[] = { + {"acos", math_acos, METH_O, math_acos_doc}, + {"acosh", math_acosh, METH_O, math_acosh_doc}, + {"asin", math_asin, METH_O, math_asin_doc}, + {"asinh", math_asinh, METH_O, math_asinh_doc}, + {"atan", math_atan, METH_O, math_atan_doc}, + {"atan2", (PyCFunction)(void(*)(void))math_atan2, METH_FASTCALL, math_atan2_doc}, + {"atanh", math_atanh, METH_O, math_atanh_doc}, + MATH_CEIL_METHODDEF + {"copysign", (PyCFunction)(void(*)(void))math_copysign, METH_FASTCALL, math_copysign_doc}, + {"cos", math_cos, METH_O, math_cos_doc}, + {"cosh", math_cosh, METH_O, math_cosh_doc}, + MATH_DEGREES_METHODDEF + MATH_DIST_METHODDEF + {"erf", math_erf, METH_O, math_erf_doc}, + {"erfc", math_erfc, METH_O, math_erfc_doc}, + {"exp", math_exp, METH_O, math_exp_doc}, + {"expm1", math_expm1, METH_O, math_expm1_doc}, + {"fabs", math_fabs, METH_O, math_fabs_doc}, + MATH_FACTORIAL_METHODDEF + MATH_FLOOR_METHODDEF + MATH_FMOD_METHODDEF + MATH_FREXP_METHODDEF + MATH_FSUM_METHODDEF + {"gamma", math_gamma, METH_O, math_gamma_doc}, + MATH_GCD_METHODDEF + {"hypot", (PyCFunction)(void(*)(void))math_hypot, METH_FASTCALL, math_hypot_doc}, + MATH_ISCLOSE_METHODDEF + MATH_ISFINITE_METHODDEF + MATH_ISINF_METHODDEF + MATH_ISNAN_METHODDEF + MATH_ISQRT_METHODDEF + MATH_LDEXP_METHODDEF + {"lgamma", math_lgamma, METH_O, math_lgamma_doc}, + MATH_LOG_METHODDEF + {"log1p", math_log1p, METH_O, math_log1p_doc}, + MATH_LOG10_METHODDEF + MATH_LOG2_METHODDEF + MATH_MODF_METHODDEF + MATH_POW_METHODDEF + MATH_RADIANS_METHODDEF + {"remainder", (PyCFunction)(void(*)(void))math_remainder, METH_FASTCALL, math_remainder_doc}, + {"sin", math_sin, METH_O, math_sin_doc}, + {"sinh", math_sinh, METH_O, math_sinh_doc}, + {"sqrt", math_sqrt, METH_O, math_sqrt_doc}, + {"tan", math_tan, METH_O, math_tan_doc}, + {"tanh", math_tanh, METH_O, math_tanh_doc}, + MATH_TRUNC_METHODDEF + MATH_PROD_METHODDEF + MATH_PERM_METHODDEF + MATH_COMB_METHODDEF + MATH_NEXTAFTER_METHODDEF + MATH_ULP_METHODDEF + {NULL, NULL} /* sentinel */ +}; + + +PyDoc_STRVAR(module_doc, +"This module provides access to the mathematical functions\n" +"defined by the C standard."); + + +static struct PyModuleDef mathmodule = { + PyModuleDef_HEAD_INIT, + "math", + module_doc, + -1, + math_methods, + NULL, + NULL, + NULL, + NULL +}; + +PyMODINIT_FUNC +PyInit_math(void) +{ + PyObject *m; + + m = PyModule_Create(&mathmodule); + if (m == NULL) + goto finally; + + PyModule_AddObject(m, "pi", PyFloat_FromDouble(Py_MATH_PI)); + PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E)); + PyModule_AddObject(m, "tau", PyFloat_FromDouble(Py_MATH_TAU)); /* 2pi */ + PyModule_AddObject(m, "inf", PyFloat_FromDouble(m_inf())); +#if !defined(PY_NO_SHORT_FLOAT_REPR) || defined(Py_NAN) + PyModule_AddObject(m, "nan", PyFloat_FromDouble(m_nan())); +#endif + + finally: + return m; } -/*[clinic end generated code: output=f8daa185c043a7b7 input=a9049054013a1b77]*/ From 76bee942fb7b3e24774c6fcdf8ac70ce00b41506 Mon Sep 17 00:00:00 2001 From: ananthan-123 Date: Wed, 19 Feb 2020 21:30:30 +0530 Subject: [PATCH 16/18] Update mathmodule.c --- Modules/mathmodule.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c index 379e16643743ce..f74b7e1a34203b 100644 --- a/Modules/mathmodule.c +++ b/Modules/mathmodule.c @@ -2027,7 +2027,7 @@ least common multiple of x and y static PyObject * math_lcm_impl(PyObject *module, PyObject *a, PyObject *b) -/*[clinic end generated code: output=6f83fb6d671074ba input=bd41b785dc2a4ff1]*/ +/*[clinic end generated code: output=6f83fb6d671074ba input=efb3d7b7334b7118]*/ { PyObject *g, *m, *f, *ab; From 8e3325e5db3b22f599840b4348c650ebd44b14df Mon Sep 17 00:00:00 2001 From: ananthan-123 Date: Wed, 19 Feb 2020 22:42:22 +0530 Subject: [PATCH 17/18] Update mathmodule.c.h --- Modules/clinic/mathmodule.c.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Modules/clinic/mathmodule.c.h b/Modules/clinic/mathmodule.c.h index 93158246d389b5..4bd659a1b09ae4 100644 --- a/Modules/clinic/mathmodule.c.h +++ b/Modules/clinic/mathmodule.c.h @@ -54,7 +54,7 @@ raised for division by zero and mod by zero. #include "Python.h" #include "pycore_dtoa.h" -#include "_math.h" +#include "math.h" #include "clinic/mathmodule.c.h" From d3edc08797d4c51d24917bc2a97ca7b9fc78c037 Mon Sep 17 00:00:00 2001 From: ananthan-123 Date: Wed, 19 Feb 2020 22:48:46 +0530 Subject: [PATCH 18/18] Update mathmodule.c.h --- Modules/clinic/mathmodule.c.h | 4047 ++++++--------------------------- 1 file changed, 748 insertions(+), 3299 deletions(-) diff --git a/Modules/clinic/mathmodule.c.h b/Modules/clinic/mathmodule.c.h index 4bd659a1b09ae4..df45a1a0c5e855 100644 --- a/Modules/clinic/mathmodule.c.h +++ b/Modules/clinic/mathmodule.c.h @@ -1,3479 +1,928 @@ -/* Math module -- standard C math library functions, pi and e */ +/*[clinic input] +preserve +[clinic start generated code]*/ -/* Here are some comments from Tim Peters, extracted from the - discussion attached to http://bugs.python.org/issue1640. They - describe the general aims of the math module with respect to - special values, IEEE-754 floating-point exceptions, and Python - exceptions. +PyDoc_STRVAR(math_gcd__doc__, +"gcd($module, x, y, /)\n" +"--\n" +"\n" +"greatest common divisor of x and y"); -These are the "spirit of 754" rules: +#define MATH_GCD_METHODDEF \ + {"gcd", (PyCFunction)(void(*)(void))math_gcd, METH_FASTCALL, math_gcd__doc__}, -1. If the mathematical result is a real number, but of magnitude too -large to approximate by a machine float, overflow is signaled and the -result is an infinity (with the appropriate sign). +static PyObject * +math_gcd_impl(PyObject *module, PyObject *a, PyObject *b); -2. If the mathematical result is a real number, but of magnitude too -small to approximate by a machine float, underflow is signaled and the -result is a zero (with the appropriate sign). +static PyObject * +math_gcd(PyObject *module, PyObject *const *args, Py_ssize_t nargs) +{ + PyObject *return_value = NULL; + PyObject *a; + PyObject *b; -3. At a singularity (a value x such that the limit of f(y) as y -approaches x exists and is an infinity), "divide by zero" is signaled -and the result is an infinity (with the appropriate sign). This is -complicated a little by that the left-side and right-side limits may -not be the same; e.g., 1/x approaches +inf or -inf as x approaches 0 -from the positive or negative directions. In that specific case, the -sign of the zero determines the result of 1/0. + if (!_PyArg_CheckPositional("gcd", nargs, 2, 2)) { + goto exit; + } + a = args[0]; + b = args[1]; + return_value = math_gcd_impl(module, a, b); -4. At a point where a function has no defined result in the extended -reals (i.e., the reals plus an infinity or two), invalid operation is -signaled and a NaN is returned. +exit: + return return_value; +} -And these are what Python has historically /tried/ to do (but not -always successfully, as platform libm behavior varies a lot): +PyDoc_STRVAR(math_ceil__doc__, +"ceil($module, x, /)\n" +"--\n" +"\n" +"Return the ceiling of x as an Integral.\n" +"\n" +"This is the smallest integer >= x."); -For #1, raise OverflowError. +#define MATH_CEIL_METHODDEF \ + {"ceil", (PyCFunction)math_ceil, METH_O, math_ceil__doc__}, -For #2, return a zero (with the appropriate sign if that happens by -accident ;-)). +PyDoc_STRVAR(math_floor__doc__, +"floor($module, x, /)\n" +"--\n" +"\n" +"Return the floor of x as an Integral.\n" +"\n" +"This is the largest integer <= x."); -For #3 and #4, raise ValueError. It may have made sense to raise -Python's ZeroDivisionError in #3, but historically that's only been -raised for division by zero and mod by zero. +#define MATH_FLOOR_METHODDEF \ + {"floor", (PyCFunction)math_floor, METH_O, math_floor__doc__}, -*/ +PyDoc_STRVAR(math_fsum__doc__, +"fsum($module, seq, /)\n" +"--\n" +"\n" +"Return an accurate floating point sum of values in the iterable seq.\n" +"\n" +"Assumes IEEE-754 floating point arithmetic."); -/* - In general, on an IEEE-754 platform the aim is to follow the C99 - standard, including Annex 'F', whenever possible. Where the - standard recommends raising the 'divide-by-zero' or 'invalid' - floating-point exceptions, Python should raise a ValueError. Where - the standard recommends raising 'overflow', Python should raise an - OverflowError. In all other circumstances a value should be - returned. - */ +#define MATH_FSUM_METHODDEF \ + {"fsum", (PyCFunction)math_fsum, METH_O, math_fsum__doc__}, -#include "Python.h" -#include "pycore_dtoa.h" -#include "math.h" +PyDoc_STRVAR(math_isqrt__doc__, +"isqrt($module, n, /)\n" +"--\n" +"\n" +"Return the integer part of the square root of the input."); -#include "clinic/mathmodule.c.h" +#define MATH_ISQRT_METHODDEF \ + {"isqrt", (PyCFunction)math_isqrt, METH_O, math_isqrt__doc__}, -/*[clinic input] -module math -[clinic start generated code]*/ -/*[clinic end generated code: output=da39a3ee5e6b4b0d input=76bc7002685dd942]*/ - - -/* - sin(pi*x), giving accurate results for all finite x (especially x - integral or close to an integer). This is here for use in the - reflection formula for the gamma function. It conforms to IEEE - 754-2008 for finite arguments, but not for infinities or nans. -*/ - -static const double pi = 3.141592653589793238462643383279502884197; -static const double logpi = 1.144729885849400174143427351353058711647; -#if !defined(HAVE_ERF) || !defined(HAVE_ERFC) -static const double sqrtpi = 1.772453850905516027298167483341145182798; -#endif /* !defined(HAVE_ERF) || !defined(HAVE_ERFC) */ - - -/* Version of PyFloat_AsDouble() with in-line fast paths - for exact floats and integers. Gives a substantial - speed improvement for extracting float arguments. -*/ - -#define ASSIGN_DOUBLE(target_var, obj, error_label) \ - if (PyFloat_CheckExact(obj)) { \ - target_var = PyFloat_AS_DOUBLE(obj); \ - } \ - else if (PyLong_CheckExact(obj)) { \ - target_var = PyLong_AsDouble(obj); \ - if (target_var == -1.0 && PyErr_Occurred()) { \ - goto error_label; \ - } \ - } \ - else { \ - target_var = PyFloat_AsDouble(obj); \ - if (target_var == -1.0 && PyErr_Occurred()) { \ - goto error_label; \ - } \ - } +PyDoc_STRVAR(math_factorial__doc__, +"factorial($module, x, /)\n" +"--\n" +"\n" +"Find x!.\n" +"\n" +"Raise a ValueError if x is negative or non-integral."); -static double -m_sinpi(double x) +#define MATH_FACTORIAL_METHODDEF \ + {"factorial", (PyCFunction)math_factorial, METH_O, math_factorial__doc__}, + +PyDoc_STRVAR(math_lcm__doc__, +"lcm($module, x, y, /)\n" +"--\n" +"\n" +"least common multiple of x and y"); + +#define MATH_LCM_METHODDEF \ + {"lcm", (PyCFunction)(void(*)(void))math_lcm, METH_FASTCALL, math_lcm__doc__}, + +static PyObject * +math_lcm_impl(PyObject *module, PyObject *a, PyObject *b); + +static PyObject * +math_lcm(PyObject *module, PyObject *const *args, Py_ssize_t nargs) { - double y, r; - int n; - /* this function should only ever be called for finite arguments */ - assert(Py_IS_FINITE(x)); - y = fmod(fabs(x), 2.0); - n = (int)round(2.0*y); - assert(0 <= n && n <= 4); - switch (n) { - case 0: - r = sin(pi*y); - break; - case 1: - r = cos(pi*(y-0.5)); - break; - case 2: - /* N.B. -sin(pi*(y-1.0)) is *not* equivalent: it would give - -0.0 instead of 0.0 when y == 1.0. */ - r = sin(pi*(1.0-y)); - break; - case 3: - r = -cos(pi*(y-1.5)); - break; - case 4: - r = sin(pi*(y-2.0)); - break; - default: - Py_UNREACHABLE(); + PyObject *return_value = NULL; + PyObject *a; + PyObject *b; + + if (!_PyArg_CheckPositional("lcm", nargs, 2, 2)) { + goto exit; } - return copysign(1.0, x)*r; + a = args[0]; + b = args[1]; + return_value = math_lcm_impl(module, a, b); + +exit: + return return_value; } -/* Implementation of the real gamma function. In extensive but non-exhaustive - random tests, this function proved accurate to within <= 10 ulps across the - entire float domain. Note that accuracy may depend on the quality of the - system math functions, the pow function in particular. Special cases - follow C99 annex F. The parameters and method are tailored to platforms - whose double format is the IEEE 754 binary64 format. - - Method: for x > 0.0 we use the Lanczos approximation with parameters N=13 - and g=6.024680040776729583740234375; these parameters are amongst those - used by the Boost library. Following Boost (again), we re-express the - Lanczos sum as a rational function, and compute it that way. The - coefficients below were computed independently using MPFR, and have been - double-checked against the coefficients in the Boost source code. - - For x < 0.0 we use the reflection formula. - - There's one minor tweak that deserves explanation: Lanczos' formula for - Gamma(x) involves computing pow(x+g-0.5, x-0.5) / exp(x+g-0.5). For many x - values, x+g-0.5 can be represented exactly. However, in cases where it - can't be represented exactly the small error in x+g-0.5 can be magnified - significantly by the pow and exp calls, especially for large x. A cheap - correction is to multiply by (1 + e*g/(x+g-0.5)), where e is the error - involved in the computation of x+g-0.5 (that is, e = computed value of - x+g-0.5 - exact value of x+g-0.5). Here's the proof: - - Correction factor - ----------------- - Write x+g-0.5 = y-e, where y is exactly representable as an IEEE 754 - double, and e is tiny. Then: - - pow(x+g-0.5,x-0.5)/exp(x+g-0.5) = pow(y-e, x-0.5)/exp(y-e) - = pow(y, x-0.5)/exp(y) * C, - - where the correction_factor C is given by - - C = pow(1-e/y, x-0.5) * exp(e) - - Since e is tiny, pow(1-e/y, x-0.5) ~ 1-(x-0.5)*e/y, and exp(x) ~ 1+e, so: - - C ~ (1-(x-0.5)*e/y) * (1+e) ~ 1 + e*(y-(x-0.5))/y - - But y-(x-0.5) = g+e, and g+e ~ g. So we get C ~ 1 + e*g/y, and - - pow(x+g-0.5,x-0.5)/exp(x+g-0.5) ~ pow(y, x-0.5)/exp(y) * (1 + e*g/y), - - Note that for accuracy, when computing r*C it's better to do - - r + e*g/y*r; - - than - - r * (1 + e*g/y); - - since the addition in the latter throws away most of the bits of - information in e*g/y. -*/ - -#define LANCZOS_N 13 -static const double lanczos_g = 6.024680040776729583740234375; -static const double lanczos_g_minus_half = 5.524680040776729583740234375; -static const double lanczos_num_coeffs[LANCZOS_N] = { - 23531376880.410759688572007674451636754734846804940, - 42919803642.649098768957899047001988850926355848959, - 35711959237.355668049440185451547166705960488635843, - 17921034426.037209699919755754458931112671403265390, - 6039542586.3520280050642916443072979210699388420708, - 1439720407.3117216736632230727949123939715485786772, - 248874557.86205415651146038641322942321632125127801, - 31426415.585400194380614231628318205362874684987640, - 2876370.6289353724412254090516208496135991145378768, - 186056.26539522349504029498971604569928220784236328, - 8071.6720023658162106380029022722506138218516325024, - 210.82427775157934587250973392071336271166969580291, - 2.5066282746310002701649081771338373386264310793408 -}; - -/* denominator is x*(x+1)*...*(x+LANCZOS_N-2) */ -static const double lanczos_den_coeffs[LANCZOS_N] = { - 0.0, 39916800.0, 120543840.0, 150917976.0, 105258076.0, 45995730.0, - 13339535.0, 2637558.0, 357423.0, 32670.0, 1925.0, 66.0, 1.0}; - -/* gamma values for small positive integers, 1 though NGAMMA_INTEGRAL */ -#define NGAMMA_INTEGRAL 23 -static const double gamma_integral[NGAMMA_INTEGRAL] = { - 1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0, 362880.0, - 3628800.0, 39916800.0, 479001600.0, 6227020800.0, 87178291200.0, - 1307674368000.0, 20922789888000.0, 355687428096000.0, - 6402373705728000.0, 121645100408832000.0, 2432902008176640000.0, - 51090942171709440000.0, 1124000727777607680000.0, -}; - -/* Lanczos' sum L_g(x), for positive x */ +PyDoc_STRVAR(math_trunc__doc__, +"trunc($module, x, /)\n" +"--\n" +"\n" +"Truncates the Real x to the nearest Integral toward 0.\n" +"\n" +"Uses the __trunc__ magic method."); -static double -lanczos_sum(double x) +#define MATH_TRUNC_METHODDEF \ + {"trunc", (PyCFunction)math_trunc, METH_O, math_trunc__doc__}, + +PyDoc_STRVAR(math_frexp__doc__, +"frexp($module, x, /)\n" +"--\n" +"\n" +"Return the mantissa and exponent of x, as pair (m, e).\n" +"\n" +"m is a float and e is an int, such that x = m * 2.**e.\n" +"If x is 0, m and e are both 0. Else 0.5 <= abs(m) < 1.0."); + +#define MATH_FREXP_METHODDEF \ + {"frexp", (PyCFunction)math_frexp, METH_O, math_frexp__doc__}, + +static PyObject * +math_frexp_impl(PyObject *module, double x); + +static PyObject * +math_frexp(PyObject *module, PyObject *arg) { - double num = 0.0, den = 0.0; - int i; - assert(x > 0.0); - /* evaluate the rational function lanczos_sum(x). For large - x, the obvious algorithm risks overflow, so we instead - rescale the denominator and numerator of the rational - function by x**(1-LANCZOS_N) and treat this as a - rational function in 1/x. This also reduces the error for - larger x values. The choice of cutoff point (5.0 below) is - somewhat arbitrary; in tests, smaller cutoff values than - this resulted in lower accuracy. */ - if (x < 5.0) { - for (i = LANCZOS_N; --i >= 0; ) { - num = num * x + lanczos_num_coeffs[i]; - den = den * x + lanczos_den_coeffs[i]; - } + PyObject *return_value = NULL; + double x; + + if (PyFloat_CheckExact(arg)) { + x = PyFloat_AS_DOUBLE(arg); } - else { - for (i = 0; i < LANCZOS_N; i++) { - num = num / x + lanczos_num_coeffs[i]; - den = den / x + lanczos_den_coeffs[i]; + else + { + x = PyFloat_AsDouble(arg); + if (x == -1.0 && PyErr_Occurred()) { + goto exit; } } - return num/den; -} - -/* Constant for +infinity, generated in the same way as float('inf'). */ + return_value = math_frexp_impl(module, x); -static double -m_inf(void) -{ -#ifndef PY_NO_SHORT_FLOAT_REPR - return _Py_dg_infinity(0); -#else - return Py_HUGE_VAL; -#endif +exit: + return return_value; } -/* Constant nan value, generated in the same way as float('nan'). */ -/* We don't currently assume that Py_NAN is defined everywhere. */ +PyDoc_STRVAR(math_ldexp__doc__, +"ldexp($module, x, i, /)\n" +"--\n" +"\n" +"Return x * (2**i).\n" +"\n" +"This is essentially the inverse of frexp()."); -#if !defined(PY_NO_SHORT_FLOAT_REPR) || defined(Py_NAN) - -static double -m_nan(void) -{ -#ifndef PY_NO_SHORT_FLOAT_REPR - return _Py_dg_stdnan(0); -#else - return Py_NAN; -#endif -} +#define MATH_LDEXP_METHODDEF \ + {"ldexp", (PyCFunction)(void(*)(void))math_ldexp, METH_FASTCALL, math_ldexp__doc__}, -#endif +static PyObject * +math_ldexp_impl(PyObject *module, double x, PyObject *i); -static double -m_tgamma(double x) +static PyObject * +math_ldexp(PyObject *module, PyObject *const *args, Py_ssize_t nargs) { - double absx, r, y, z, sqrtpow; - - /* special cases */ - if (!Py_IS_FINITE(x)) { - if (Py_IS_NAN(x) || x > 0.0) - return x; /* tgamma(nan) = nan, tgamma(inf) = inf */ - else { - errno = EDOM; - return Py_NAN; /* tgamma(-inf) = nan, invalid */ - } - } - if (x == 0.0) { - errno = EDOM; - /* tgamma(+-0.0) = +-inf, divide-by-zero */ - return copysign(Py_HUGE_VAL, x); - } + PyObject *return_value = NULL; + double x; + PyObject *i; - /* integer arguments */ - if (x == floor(x)) { - if (x < 0.0) { - errno = EDOM; /* tgamma(n) = nan, invalid for */ - return Py_NAN; /* negative integers n */ - } - if (x <= NGAMMA_INTEGRAL) - return gamma_integral[(int)x - 1]; + if (!_PyArg_CheckPositional("ldexp", nargs, 2, 2)) { + goto exit; } - absx = fabs(x); - - /* tiny arguments: tgamma(x) ~ 1/x for x near 0 */ - if (absx < 1e-20) { - r = 1.0/x; - if (Py_IS_INFINITY(r)) - errno = ERANGE; - return r; + if (PyFloat_CheckExact(args[0])) { + x = PyFloat_AS_DOUBLE(args[0]); } - - /* large arguments: assuming IEEE 754 doubles, tgamma(x) overflows for - x > 200, and underflows to +-0.0 for x < -200, not a negative - integer. */ - if (absx > 200.0) { - if (x < 0.0) { - return 0.0/m_sinpi(x); - } - else { - errno = ERANGE; - return Py_HUGE_VAL; + else + { + x = PyFloat_AsDouble(args[0]); + if (x == -1.0 && PyErr_Occurred()) { + goto exit; } } + i = args[1]; + return_value = math_ldexp_impl(module, x, i); - y = absx + lanczos_g_minus_half; - /* compute error in sum */ - if (absx > lanczos_g_minus_half) { - /* note: the correction can be foiled by an optimizing - compiler that (incorrectly) thinks that an expression like - a + b - a - b can be optimized to 0.0. This shouldn't - happen in a standards-conforming compiler. */ - double q = y - absx; - z = q - lanczos_g_minus_half; - } - else { - double q = y - lanczos_g_minus_half; - z = q - absx; - } - z = z * lanczos_g / y; - if (x < 0.0) { - r = -pi / m_sinpi(absx) / absx * exp(y) / lanczos_sum(absx); - r -= z * r; - if (absx < 140.0) { - r /= pow(y, absx - 0.5); - } - else { - sqrtpow = pow(y, absx / 2.0 - 0.25); - r /= sqrtpow; - r /= sqrtpow; - } - } - else { - r = lanczos_sum(absx) / exp(y); - r += z * r; - if (absx < 140.0) { - r *= pow(y, absx - 0.5); - } - else { - sqrtpow = pow(y, absx / 2.0 - 0.25); - r *= sqrtpow; - r *= sqrtpow; - } - } - if (Py_IS_INFINITY(r)) - errno = ERANGE; - return r; +exit: + return return_value; } -/* - lgamma: natural log of the absolute value of the Gamma function. - For large arguments, Lanczos' formula works extremely well here. -*/ +PyDoc_STRVAR(math_modf__doc__, +"modf($module, x, /)\n" +"--\n" +"\n" +"Return the fractional and integer parts of x.\n" +"\n" +"Both results carry the sign of x and are floats."); -static double -m_lgamma(double x) +#define MATH_MODF_METHODDEF \ + {"modf", (PyCFunction)math_modf, METH_O, math_modf__doc__}, + +static PyObject * +math_modf_impl(PyObject *module, double x); + +static PyObject * +math_modf(PyObject *module, PyObject *arg) { - double r; - double absx; + PyObject *return_value = NULL; + double x; - /* special cases */ - if (!Py_IS_FINITE(x)) { - if (Py_IS_NAN(x)) - return x; /* lgamma(nan) = nan */ - else - return Py_HUGE_VAL; /* lgamma(+-inf) = +inf */ + if (PyFloat_CheckExact(arg)) { + x = PyFloat_AS_DOUBLE(arg); } - - /* integer arguments */ - if (x == floor(x) && x <= 2.0) { - if (x <= 0.0) { - errno = EDOM; /* lgamma(n) = inf, divide-by-zero for */ - return Py_HUGE_VAL; /* integers n <= 0 */ - } - else { - return 0.0; /* lgamma(1) = lgamma(2) = 0.0 */ + else + { + x = PyFloat_AsDouble(arg); + if (x == -1.0 && PyErr_Occurred()) { + goto exit; } } + return_value = math_modf_impl(module, x); - absx = fabs(x); - /* tiny arguments: lgamma(x) ~ -log(fabs(x)) for small x */ - if (absx < 1e-20) - return -log(absx); - - /* Lanczos' formula. We could save a fraction of a ulp in accuracy by - having a second set of numerator coefficients for lanczos_sum that - absorbed the exp(-lanczos_g) term, and throwing out the lanczos_g - subtraction below; it's probably not worth it. */ - r = log(lanczos_sum(absx)) - lanczos_g; - r += (absx - 0.5) * (log(absx + lanczos_g - 0.5) - 1); - if (x < 0.0) - /* Use reflection formula to get value for negative x. */ - r = logpi - log(fabs(m_sinpi(absx))) - log(absx) - r; - if (Py_IS_INFINITY(r)) - errno = ERANGE; - return r; +exit: + return return_value; } -#if !defined(HAVE_ERF) || !defined(HAVE_ERFC) - -/* - Implementations of the error function erf(x) and the complementary error - function erfc(x). +PyDoc_STRVAR(math_log__doc__, +"log(x, [base=math.e])\n" +"Return the logarithm of x to the given base.\n" +"\n" +"If the base not specified, returns the natural logarithm (base e) of x."); - Method: we use a series approximation for erf for small x, and a continued - fraction approximation for erfc(x) for larger x; - combined with the relations erf(-x) = -erf(x) and erfc(x) = 1.0 - erf(x), - this gives us erf(x) and erfc(x) for all x. +#define MATH_LOG_METHODDEF \ + {"log", (PyCFunction)math_log, METH_VARARGS, math_log__doc__}, - The series expansion used is: - - erf(x) = x*exp(-x*x)/sqrt(pi) * [ - 2/1 + 4/3 x**2 + 8/15 x**4 + 16/105 x**6 + ...] +static PyObject * +math_log_impl(PyObject *module, PyObject *x, int group_right_1, + PyObject *base); - The coefficient of x**(2k-2) here is 4**k*factorial(k)/factorial(2*k). - This series converges well for smallish x, but slowly for larger x. +static PyObject * +math_log(PyObject *module, PyObject *args) +{ + PyObject *return_value = NULL; + PyObject *x; + int group_right_1 = 0; + PyObject *base = NULL; - The continued fraction expansion used is: + switch (PyTuple_GET_SIZE(args)) { + case 1: + if (!PyArg_ParseTuple(args, "O:log", &x)) { + goto exit; + } + break; + case 2: + if (!PyArg_ParseTuple(args, "OO:log", &x, &base)) { + goto exit; + } + group_right_1 = 1; + break; + default: + PyErr_SetString(PyExc_TypeError, "math.log requires 1 to 2 arguments"); + goto exit; + } + return_value = math_log_impl(module, x, group_right_1, base); - erfc(x) = x*exp(-x*x)/sqrt(pi) * [1/(0.5 + x**2 -) 0.5/(2.5 + x**2 - ) - 3.0/(4.5 + x**2 - ) 7.5/(6.5 + x**2 - ) ...] +exit: + return return_value; +} - after the first term, the general term has the form: +PyDoc_STRVAR(math_log2__doc__, +"log2($module, x, /)\n" +"--\n" +"\n" +"Return the base 2 logarithm of x."); - k*(k-0.5)/(2*k+0.5 + x**2 - ...). +#define MATH_LOG2_METHODDEF \ + {"log2", (PyCFunction)math_log2, METH_O, math_log2__doc__}, - This expansion converges fast for larger x, but convergence becomes - infinitely slow as x approaches 0.0. The (somewhat naive) continued - fraction evaluation algorithm used below also risks overflow for large x; - but for large x, erfc(x) == 0.0 to within machine precision. (For - example, erfc(30.0) is approximately 2.56e-393). +PyDoc_STRVAR(math_log10__doc__, +"log10($module, x, /)\n" +"--\n" +"\n" +"Return the base 10 logarithm of x."); - Parameters: use series expansion for abs(x) < ERF_SERIES_CUTOFF and - continued fraction expansion for ERF_SERIES_CUTOFF <= abs(x) < - ERFC_CONTFRAC_CUTOFF. ERFC_SERIES_TERMS and ERFC_CONTFRAC_TERMS are the - numbers of terms to use for the relevant expansions. */ +#define MATH_LOG10_METHODDEF \ + {"log10", (PyCFunction)math_log10, METH_O, math_log10__doc__}, -#define ERF_SERIES_CUTOFF 1.5 -#define ERF_SERIES_TERMS 25 -#define ERFC_CONTFRAC_CUTOFF 30.0 -#define ERFC_CONTFRAC_TERMS 50 +PyDoc_STRVAR(math_fmod__doc__, +"fmod($module, x, y, /)\n" +"--\n" +"\n" +"Return fmod(x, y), according to platform C.\n" +"\n" +"x % y may differ."); -/* - Error function, via power series. +#define MATH_FMOD_METHODDEF \ + {"fmod", (PyCFunction)(void(*)(void))math_fmod, METH_FASTCALL, math_fmod__doc__}, - Given a finite float x, return an approximation to erf(x). - Converges reasonably fast for small x. -*/ +static PyObject * +math_fmod_impl(PyObject *module, double x, double y); -static double -m_erf_series(double x) +static PyObject * +math_fmod(PyObject *module, PyObject *const *args, Py_ssize_t nargs) { - double x2, acc, fk, result; - int i, saved_errno; - - x2 = x * x; - acc = 0.0; - fk = (double)ERF_SERIES_TERMS + 0.5; - for (i = 0; i < ERF_SERIES_TERMS; i++) { - acc = 2.0 + x2 * acc / fk; - fk -= 1.0; - } - /* Make sure the exp call doesn't affect errno; - see m_erfc_contfrac for more. */ - saved_errno = errno; - result = acc * x * exp(-x2) / sqrtpi; - errno = saved_errno; - return result; -} - -/* - Complementary error function, via continued fraction expansion. - - Given a positive float x, return an approximation to erfc(x). Converges - reasonably fast for x large (say, x > 2.0), and should be safe from - overflow if x and nterms are not too large. On an IEEE 754 machine, with x - <= 30.0, we're safe up to nterms = 100. For x >= 30.0, erfc(x) is smaller - than the smallest representable nonzero float. */ + PyObject *return_value = NULL; + double x; + double y; -static double -m_erfc_contfrac(double x) -{ - double x2, a, da, p, p_last, q, q_last, b, result; - int i, saved_errno; - - if (x >= ERFC_CONTFRAC_CUTOFF) - return 0.0; - - x2 = x*x; - a = 0.0; - da = 0.5; - p = 1.0; p_last = 0.0; - q = da + x2; q_last = 1.0; - for (i = 0; i < ERFC_CONTFRAC_TERMS; i++) { - double temp; - a += da; - da += 2.0; - b = da + x2; - temp = p; p = b*p - a*p_last; p_last = temp; - temp = q; q = b*q - a*q_last; q_last = temp; + if (!_PyArg_CheckPositional("fmod", nargs, 2, 2)) { + goto exit; + } + if (PyFloat_CheckExact(args[0])) { + x = PyFloat_AS_DOUBLE(args[0]); + } + else + { + x = PyFloat_AsDouble(args[0]); + if (x == -1.0 && PyErr_Occurred()) { + goto exit; + } + } + if (PyFloat_CheckExact(args[1])) { + y = PyFloat_AS_DOUBLE(args[1]); + } + else + { + y = PyFloat_AsDouble(args[1]); + if (y == -1.0 && PyErr_Occurred()) { + goto exit; + } } - /* Issue #8986: On some platforms, exp sets errno on underflow to zero; - save the current errno value so that we can restore it later. */ - saved_errno = errno; - result = p / q * x * exp(-x2) / sqrtpi; - errno = saved_errno; - return result; + return_value = math_fmod_impl(module, x, y); + +exit: + return return_value; } -#endif /* !defined(HAVE_ERF) || !defined(HAVE_ERFC) */ +PyDoc_STRVAR(math_dist__doc__, +"dist($module, p, q, /)\n" +"--\n" +"\n" +"Return the Euclidean distance between two points p and q.\n" +"\n" +"The points should be specified as sequences (or iterables) of\n" +"coordinates. Both inputs must have the same dimension.\n" +"\n" +"Roughly equivalent to:\n" +" sqrt(sum((px - qx) ** 2.0 for px, qx in zip(p, q)))"); -/* Error function erf(x), for general x */ +#define MATH_DIST_METHODDEF \ + {"dist", (PyCFunction)(void(*)(void))math_dist, METH_FASTCALL, math_dist__doc__}, -static double -m_erf(double x) +static PyObject * +math_dist_impl(PyObject *module, PyObject *p, PyObject *q); + +static PyObject * +math_dist(PyObject *module, PyObject *const *args, Py_ssize_t nargs) { -#ifdef HAVE_ERF - return erf(x); -#else - double absx, cf; - - if (Py_IS_NAN(x)) - return x; - absx = fabs(x); - if (absx < ERF_SERIES_CUTOFF) - return m_erf_series(x); - else { - cf = m_erfc_contfrac(absx); - return x > 0.0 ? 1.0 - cf : cf - 1.0; + PyObject *return_value = NULL; + PyObject *p; + PyObject *q; + + if (!_PyArg_CheckPositional("dist", nargs, 2, 2)) { + goto exit; } -#endif + p = args[0]; + q = args[1]; + return_value = math_dist_impl(module, p, q); + +exit: + return return_value; } -/* Complementary error function erfc(x), for general x. */ +PyDoc_STRVAR(math_pow__doc__, +"pow($module, x, y, /)\n" +"--\n" +"\n" +"Return x**y (x to the power of y)."); -static double -m_erfc(double x) -{ -#ifdef HAVE_ERFC - return erfc(x); -#else - double absx, cf; - - if (Py_IS_NAN(x)) - return x; - absx = fabs(x); - if (absx < ERF_SERIES_CUTOFF) - return 1.0 - m_erf_series(x); - else { - cf = m_erfc_contfrac(absx); - return x > 0.0 ? cf : 2.0 - cf; - } -#endif -} +#define MATH_POW_METHODDEF \ + {"pow", (PyCFunction)(void(*)(void))math_pow, METH_FASTCALL, math_pow__doc__}, -/* - wrapper for atan2 that deals directly with special cases before - delegating to the platform libm for the remaining cases. This - is necessary to get consistent behaviour across platforms. - Windows, FreeBSD and alpha Tru64 are amongst platforms that don't - always follow C99. -*/ +static PyObject * +math_pow_impl(PyObject *module, double x, double y); -static double -m_atan2(double y, double x) +static PyObject * +math_pow(PyObject *module, PyObject *const *args, Py_ssize_t nargs) { - if (Py_IS_NAN(x) || Py_IS_NAN(y)) - return Py_NAN; - if (Py_IS_INFINITY(y)) { - if (Py_IS_INFINITY(x)) { - if (copysign(1., x) == 1.) - /* atan2(+-inf, +inf) == +-pi/4 */ - return copysign(0.25*Py_MATH_PI, y); - else - /* atan2(+-inf, -inf) == +-pi*3/4 */ - return copysign(0.75*Py_MATH_PI, y); + PyObject *return_value = NULL; + double x; + double y; + + if (!_PyArg_CheckPositional("pow", nargs, 2, 2)) { + goto exit; + } + if (PyFloat_CheckExact(args[0])) { + x = PyFloat_AS_DOUBLE(args[0]); + } + else + { + x = PyFloat_AsDouble(args[0]); + if (x == -1.0 && PyErr_Occurred()) { + goto exit; } - /* atan2(+-inf, x) == +-pi/2 for finite x */ - return copysign(0.5*Py_MATH_PI, y); } - if (Py_IS_INFINITY(x) || y == 0.) { - if (copysign(1., x) == 1.) - /* atan2(+-y, +inf) = atan2(+-0, +x) = +-0. */ - return copysign(0., y); - else - /* atan2(+-y, -inf) = atan2(+-0., -x) = +-pi. */ - return copysign(Py_MATH_PI, y); + if (PyFloat_CheckExact(args[1])) { + y = PyFloat_AS_DOUBLE(args[1]); + } + else + { + y = PyFloat_AsDouble(args[1]); + if (y == -1.0 && PyErr_Occurred()) { + goto exit; + } } - return atan2(y, x); + return_value = math_pow_impl(module, x, y); + +exit: + return return_value; } +PyDoc_STRVAR(math_degrees__doc__, +"degrees($module, x, /)\n" +"--\n" +"\n" +"Convert angle x from radians to degrees."); -/* IEEE 754-style remainder operation: x - n*y where n*y is the nearest - multiple of y to x, taking n even in the case of a tie. Assuming an IEEE 754 - binary floating-point format, the result is always exact. */ +#define MATH_DEGREES_METHODDEF \ + {"degrees", (PyCFunction)math_degrees, METH_O, math_degrees__doc__}, -static double -m_remainder(double x, double y) -{ - /* Deal with most common case first. */ - if (Py_IS_FINITE(x) && Py_IS_FINITE(y)) { - double absx, absy, c, m, r; +static PyObject * +math_degrees_impl(PyObject *module, double x); - if (y == 0.0) { - return Py_NAN; - } +static PyObject * +math_degrees(PyObject *module, PyObject *arg) +{ + PyObject *return_value = NULL; + double x; - absx = fabs(x); - absy = fabs(y); - m = fmod(absx, absy); - - /* - Warning: some subtlety here. What we *want* to know at this point is - whether the remainder m is less than, equal to, or greater than half - of absy. However, we can't do that comparison directly because we - can't be sure that 0.5*absy is representable (the multiplication - might incur precision loss due to underflow). So instead we compare - m with the complement c = absy - m: m < 0.5*absy if and only if m < - c, and so on. The catch is that absy - m might also not be - representable, but it turns out that it doesn't matter: - - - if m > 0.5*absy then absy - m is exactly representable, by - Sterbenz's lemma, so m > c - - if m == 0.5*absy then again absy - m is exactly representable - and m == c - - if m < 0.5*absy then either (i) 0.5*absy is exactly representable, - in which case 0.5*absy < absy - m, so 0.5*absy <= c and hence m < - c, or (ii) absy is tiny, either subnormal or in the lowest normal - binade. Then absy - m is exactly representable and again m < c. - */ - - c = absy - m; - if (m < c) { - r = m; - } - else if (m > c) { - r = -c; + if (PyFloat_CheckExact(arg)) { + x = PyFloat_AS_DOUBLE(arg); + } + else + { + x = PyFloat_AsDouble(arg); + if (x == -1.0 && PyErr_Occurred()) { + goto exit; } - else { - /* - Here absx is exactly halfway between two multiples of absy, - and we need to choose the even multiple. x now has the form - - absx = n * absy + m - - for some integer n (recalling that m = 0.5*absy at this point). - If n is even we want to return m; if n is odd, we need to - return -m. + } + return_value = math_degrees_impl(module, x); - So +exit: + return return_value; +} - 0.5 * (absx - m) = (n/2) * absy +PyDoc_STRVAR(math_radians__doc__, +"radians($module, x, /)\n" +"--\n" +"\n" +"Convert angle x from degrees to radians."); - and now reducing modulo absy gives us: +#define MATH_RADIANS_METHODDEF \ + {"radians", (PyCFunction)math_radians, METH_O, math_radians__doc__}, - | m, if n is odd - fmod(0.5 * (absx - m), absy) = | - | 0, if n is even +static PyObject * +math_radians_impl(PyObject *module, double x); - Now m - 2.0 * fmod(...) gives the desired result: m - if n is even, -m if m is odd. +static PyObject * +math_radians(PyObject *module, PyObject *arg) +{ + PyObject *return_value = NULL; + double x; - Note that all steps in fmod(0.5 * (absx - m), absy) - will be computed exactly, with no rounding error - introduced. - */ - assert(m == c); - r = m - 2.0 * fmod(0.5 * (absx - m), absy); + if (PyFloat_CheckExact(arg)) { + x = PyFloat_AS_DOUBLE(arg); + } + else + { + x = PyFloat_AsDouble(arg); + if (x == -1.0 && PyErr_Occurred()) { + goto exit; } - return copysign(1.0, x) * r; } + return_value = math_radians_impl(module, x); - /* Special values. */ - if (Py_IS_NAN(x)) { - return x; - } - if (Py_IS_NAN(y)) { - return y; - } - if (Py_IS_INFINITY(x)) { - return Py_NAN; - } - assert(Py_IS_INFINITY(y)); - return x; +exit: + return return_value; } +PyDoc_STRVAR(math_isfinite__doc__, +"isfinite($module, x, /)\n" +"--\n" +"\n" +"Return True if x is neither an infinity nor a NaN, and False otherwise."); -/* - Various platforms (Solaris, OpenBSD) do nonstandard things for log(0), - log(-ve), log(NaN). Here are wrappers for log and log10 that deal with - special values directly, passing positive non-special values through to - the system log/log10. - */ +#define MATH_ISFINITE_METHODDEF \ + {"isfinite", (PyCFunction)math_isfinite, METH_O, math_isfinite__doc__}, -static double -m_log(double x) +static PyObject * +math_isfinite_impl(PyObject *module, double x); + +static PyObject * +math_isfinite(PyObject *module, PyObject *arg) { - if (Py_IS_FINITE(x)) { - if (x > 0.0) - return log(x); - errno = EDOM; - if (x == 0.0) - return -Py_HUGE_VAL; /* log(0) = -inf */ - else - return Py_NAN; /* log(-ve) = nan */ + PyObject *return_value = NULL; + double x; + + if (PyFloat_CheckExact(arg)) { + x = PyFloat_AS_DOUBLE(arg); } - else if (Py_IS_NAN(x)) - return x; /* log(nan) = nan */ - else if (x > 0.0) - return x; /* log(inf) = inf */ - else { - errno = EDOM; - return Py_NAN; /* log(-inf) = nan */ + else + { + x = PyFloat_AsDouble(arg); + if (x == -1.0 && PyErr_Occurred()) { + goto exit; + } } + return_value = math_isfinite_impl(module, x); + +exit: + return return_value; } -/* - log2: log to base 2. +PyDoc_STRVAR(math_isnan__doc__, +"isnan($module, x, /)\n" +"--\n" +"\n" +"Return True if x is a NaN (not a number), and False otherwise."); - Uses an algorithm that should: +#define MATH_ISNAN_METHODDEF \ + {"isnan", (PyCFunction)math_isnan, METH_O, math_isnan__doc__}, - (a) produce exact results for powers of 2, and - (b) give a monotonic log2 (for positive finite floats), - assuming that the system log is monotonic. -*/ +static PyObject * +math_isnan_impl(PyObject *module, double x); -static double -m_log2(double x) +static PyObject * +math_isnan(PyObject *module, PyObject *arg) { - if (!Py_IS_FINITE(x)) { - if (Py_IS_NAN(x)) - return x; /* log2(nan) = nan */ - else if (x > 0.0) - return x; /* log2(+inf) = +inf */ - else { - errno = EDOM; - return Py_NAN; /* log2(-inf) = nan, invalid-operation */ - } - } + PyObject *return_value = NULL; + double x; - if (x > 0.0) { -#ifdef HAVE_LOG2 - return log2(x); -#else - double m; - int e; - m = frexp(x, &e); - /* We want log2(m * 2**e) == log(m) / log(2) + e. Care is needed when - * x is just greater than 1.0: in that case e is 1, log(m) is negative, - * and we get significant cancellation error from the addition of - * log(m) / log(2) to e. The slight rewrite of the expression below - * avoids this problem. - */ - if (x >= 1.0) { - return log(2.0 * m) / log(2.0) + (e - 1); - } - else { - return log(m) / log(2.0) + e; - } -#endif - } - else if (x == 0.0) { - errno = EDOM; - return -Py_HUGE_VAL; /* log2(0) = -inf, divide-by-zero */ + if (PyFloat_CheckExact(arg)) { + x = PyFloat_AS_DOUBLE(arg); } - else { - errno = EDOM; - return Py_NAN; /* log2(-inf) = nan, invalid-operation */ + else + { + x = PyFloat_AsDouble(arg); + if (x == -1.0 && PyErr_Occurred()) { + goto exit; + } } -} + return_value = math_isnan_impl(module, x); -static double -m_log10(double x) -{ - if (Py_IS_FINITE(x)) { - if (x > 0.0) - return log10(x); - errno = EDOM; - if (x == 0.0) - return -Py_HUGE_VAL; /* log10(0) = -inf */ - else - return Py_NAN; /* log10(-ve) = nan */ - } - else if (Py_IS_NAN(x)) - return x; /* log10(nan) = nan */ - else if (x > 0.0) - return x; /* log10(inf) = inf */ - else { - errno = EDOM; - return Py_NAN; /* log10(-inf) = nan */ - } +exit: + return return_value; } +PyDoc_STRVAR(math_isinf__doc__, +"isinf($module, x, /)\n" +"--\n" +"\n" +"Return True if x is a positive or negative infinity, and False otherwise."); -/*[clinic input] -math.gcd - - x as a: object - y as b: object - / +#define MATH_ISINF_METHODDEF \ + {"isinf", (PyCFunction)math_isinf, METH_O, math_isinf__doc__}, -greatest common divisor of x and y -[clinic start generated code]*/ +static PyObject * +math_isinf_impl(PyObject *module, double x); static PyObject * -math_gcd_impl(PyObject *module, PyObject *a, PyObject *b) -/*[clinic end generated code: output=7b2e0c151bd7a5d8 input=c2691e57fb2a98fa]*/ +math_isinf(PyObject *module, PyObject *arg) { - PyObject *g; - - a = PyNumber_Index(a); - if (a == NULL) - return NULL; - b = PyNumber_Index(b); - if (b == NULL) { - Py_DECREF(a); - return NULL; - } - g = _PyLong_GCD(a, b); - Py_DECREF(a); - Py_DECREF(b); - return g; -} - + PyObject *return_value = NULL; + double x; -/* Call is_error when errno != 0, and where x is the result libm - * returned. is_error will usually set up an exception and return - * true (1), but may return false (0) without setting up an exception. - */ -static int -is_error(double x) -{ - int result = 1; /* presumption of guilt */ - assert(errno); /* non-zero errno is a precondition for calling */ - if (errno == EDOM) - PyErr_SetString(PyExc_ValueError, "math domain error"); - - else if (errno == ERANGE) { - /* ANSI C generally requires libm functions to set ERANGE - * on overflow, but also generally *allows* them to set - * ERANGE on underflow too. There's no consistency about - * the latter across platforms. - * Alas, C99 never requires that errno be set. - * Here we suppress the underflow errors (libm functions - * should return a zero on underflow, and +- HUGE_VAL on - * overflow, so testing the result for zero suffices to - * distinguish the cases). - * - * On some platforms (Ubuntu/ia64) it seems that errno can be - * set to ERANGE for subnormal results that do *not* underflow - * to zero. So to be safe, we'll ignore ERANGE whenever the - * function result is less than one in absolute value. - */ - if (fabs(x) < 1.0) - result = 0; - else - PyErr_SetString(PyExc_OverflowError, - "math range error"); + if (PyFloat_CheckExact(arg)) { + x = PyFloat_AS_DOUBLE(arg); } else - /* Unexpected math error */ - PyErr_SetFromErrno(PyExc_ValueError); - return result; -} + { + x = PyFloat_AsDouble(arg); + if (x == -1.0 && PyErr_Occurred()) { + goto exit; + } + } + return_value = math_isinf_impl(module, x); + +exit: + return return_value; +} + +PyDoc_STRVAR(math_isclose__doc__, +"isclose($module, /, a, b, *, rel_tol=1e-09, abs_tol=0.0)\n" +"--\n" +"\n" +"Determine whether two floating point numbers are close in value.\n" +"\n" +" rel_tol\n" +" maximum difference for being considered \"close\", relative to the\n" +" magnitude of the input values\n" +" abs_tol\n" +" maximum difference for being considered \"close\", regardless of the\n" +" magnitude of the input values\n" +"\n" +"Return True if a is close in value to b, and False otherwise.\n" +"\n" +"For the values to be considered close, the difference between them\n" +"must be smaller than at least one of the tolerances.\n" +"\n" +"-inf, inf and NaN behave similarly to the IEEE 754 Standard. That\n" +"is, NaN is not close to anything, even itself. inf and -inf are\n" +"only close to themselves."); + +#define MATH_ISCLOSE_METHODDEF \ + {"isclose", (PyCFunction)(void(*)(void))math_isclose, METH_FASTCALL|METH_KEYWORDS, math_isclose__doc__}, -/* - math_1 is used to wrap a libm function f that takes a double - argument and returns a double. - - The error reporting follows these rules, which are designed to do - the right thing on C89/C99 platforms and IEEE 754/non IEEE 754 - platforms. - - - a NaN result from non-NaN inputs causes ValueError to be raised - - an infinite result from finite inputs causes OverflowError to be - raised if can_overflow is 1, or raises ValueError if can_overflow - is 0. - - if the result is finite and errno == EDOM then ValueError is - raised - - if the result is finite and nonzero and errno == ERANGE then - OverflowError is raised - - The last rule is used to catch overflow on platforms which follow - C89 but for which HUGE_VAL is not an infinity. - - For the majority of one-argument functions these rules are enough - to ensure that Python's functions behave as specified in 'Annex F' - of the C99 standard, with the 'invalid' and 'divide-by-zero' - floating-point exceptions mapping to Python's ValueError and the - 'overflow' floating-point exception mapping to OverflowError. - math_1 only works for functions that don't have singularities *and* - the possibility of overflow; fortunately, that covers everything we - care about right now. -*/ +static int +math_isclose_impl(PyObject *module, double a, double b, double rel_tol, + double abs_tol); static PyObject * -math_1_to_whatever(PyObject *arg, double (*func) (double), - PyObject *(*from_double_func) (double), - int can_overflow) +math_isclose(PyObject *module, PyObject *const *args, Py_ssize_t nargs, PyObject *kwnames) { - double x, r; - x = PyFloat_AsDouble(arg); - if (x == -1.0 && PyErr_Occurred()) - return NULL; - errno = 0; - r = (*func)(x); - if (Py_IS_NAN(r) && !Py_IS_NAN(x)) { - PyErr_SetString(PyExc_ValueError, - "math domain error"); /* invalid arg */ - return NULL; + PyObject *return_value = NULL; + static const char * const _keywords[] = {"a", "b", "rel_tol", "abs_tol", NULL}; + static _PyArg_Parser _parser = {NULL, _keywords, "isclose", 0}; + PyObject *argsbuf[4]; + Py_ssize_t noptargs = nargs + (kwnames ? PyTuple_GET_SIZE(kwnames) : 0) - 2; + double a; + double b; + double rel_tol = 1e-09; + double abs_tol = 0.0; + int _return_value; + + args = _PyArg_UnpackKeywords(args, nargs, NULL, kwnames, &_parser, 2, 2, 0, argsbuf); + if (!args) { + goto exit; + } + if (PyFloat_CheckExact(args[0])) { + a = PyFloat_AS_DOUBLE(args[0]); + } + else + { + a = PyFloat_AsDouble(args[0]); + if (a == -1.0 && PyErr_Occurred()) { + goto exit; + } } - if (Py_IS_INFINITY(r) && Py_IS_FINITE(x)) { - if (can_overflow) - PyErr_SetString(PyExc_OverflowError, - "math range error"); /* overflow */ + if (PyFloat_CheckExact(args[1])) { + b = PyFloat_AS_DOUBLE(args[1]); + } + else + { + b = PyFloat_AsDouble(args[1]); + if (b == -1.0 && PyErr_Occurred()) { + goto exit; + } + } + if (!noptargs) { + goto skip_optional_kwonly; + } + if (args[2]) { + if (PyFloat_CheckExact(args[2])) { + rel_tol = PyFloat_AS_DOUBLE(args[2]); + } else - PyErr_SetString(PyExc_ValueError, - "math domain error"); /* singularity */ - return NULL; + { + rel_tol = PyFloat_AsDouble(args[2]); + if (rel_tol == -1.0 && PyErr_Occurred()) { + goto exit; + } + } + if (!--noptargs) { + goto skip_optional_kwonly; + } + } + if (PyFloat_CheckExact(args[3])) { + abs_tol = PyFloat_AS_DOUBLE(args[3]); + } + else + { + abs_tol = PyFloat_AsDouble(args[3]); + if (abs_tol == -1.0 && PyErr_Occurred()) { + goto exit; + } } - if (Py_IS_FINITE(r) && errno && is_error(r)) - /* this branch unnecessary on most platforms */ - return NULL; +skip_optional_kwonly: + _return_value = math_isclose_impl(module, a, b, rel_tol, abs_tol); + if ((_return_value == -1) && PyErr_Occurred()) { + goto exit; + } + return_value = PyBool_FromLong((long)_return_value); - return (*from_double_func)(r); +exit: + return return_value; } -/* variant of math_1, to be used when the function being wrapped is known to - set errno properly (that is, errno = EDOM for invalid or divide-by-zero, - errno = ERANGE for overflow). */ - -static PyObject * -math_1a(PyObject *arg, double (*func) (double)) -{ - double x, r; - x = PyFloat_AsDouble(arg); - if (x == -1.0 && PyErr_Occurred()) - return NULL; - errno = 0; - r = (*func)(x); - if (errno && is_error(r)) - return NULL; - return PyFloat_FromDouble(r); -} +PyDoc_STRVAR(math_prod__doc__, +"prod($module, iterable, /, *, start=1)\n" +"--\n" +"\n" +"Calculate the product of all the elements in the input iterable.\n" +"\n" +"The default start value for the product is 1.\n" +"\n" +"When the iterable is empty, return the start value. This function is\n" +"intended specifically for use with numeric values and may reject\n" +"non-numeric types."); -/* - math_2 is used to wrap a libm function f that takes two double - arguments and returns a double. - - The error reporting follows these rules, which are designed to do - the right thing on C89/C99 platforms and IEEE 754/non IEEE 754 - platforms. - - - a NaN result from non-NaN inputs causes ValueError to be raised - - an infinite result from finite inputs causes OverflowError to be - raised. - - if the result is finite and errno == EDOM then ValueError is - raised - - if the result is finite and nonzero and errno == ERANGE then - OverflowError is raised - - The last rule is used to catch overflow on platforms which follow - C89 but for which HUGE_VAL is not an infinity. - - For most two-argument functions (copysign, fmod, hypot, atan2) - these rules are enough to ensure that Python's functions behave as - specified in 'Annex F' of the C99 standard, with the 'invalid' and - 'divide-by-zero' floating-point exceptions mapping to Python's - ValueError and the 'overflow' floating-point exception mapping to - OverflowError. -*/ +#define MATH_PROD_METHODDEF \ + {"prod", (PyCFunction)(void(*)(void))math_prod, METH_FASTCALL|METH_KEYWORDS, math_prod__doc__}, static PyObject * -math_1(PyObject *arg, double (*func) (double), int can_overflow) -{ - return math_1_to_whatever(arg, func, PyFloat_FromDouble, can_overflow); -} +math_prod_impl(PyObject *module, PyObject *iterable, PyObject *start); static PyObject * -math_2(PyObject *const *args, Py_ssize_t nargs, - double (*func) (double, double), const char *funcname) +math_prod(PyObject *module, PyObject *const *args, Py_ssize_t nargs, PyObject *kwnames) { - double x, y, r; - if (!_PyArg_CheckPositional(funcname, nargs, 2, 2)) - return NULL; - x = PyFloat_AsDouble(args[0]); - y = PyFloat_AsDouble(args[1]); - if ((x == -1.0 || y == -1.0) && PyErr_Occurred()) - return NULL; - errno = 0; - r = (*func)(x, y); - if (Py_IS_NAN(r)) { - if (!Py_IS_NAN(x) && !Py_IS_NAN(y)) - errno = EDOM; - else - errno = 0; + PyObject *return_value = NULL; + static const char * const _keywords[] = {"", "start", NULL}; + static _PyArg_Parser _parser = {NULL, _keywords, "prod", 0}; + PyObject *argsbuf[2]; + Py_ssize_t noptargs = nargs + (kwnames ? PyTuple_GET_SIZE(kwnames) : 0) - 1; + PyObject *iterable; + PyObject *start = NULL; + + args = _PyArg_UnpackKeywords(args, nargs, NULL, kwnames, &_parser, 1, 1, 0, argsbuf); + if (!args) { + goto exit; } - else if (Py_IS_INFINITY(r)) { - if (Py_IS_FINITE(x) && Py_IS_FINITE(y)) - errno = ERANGE; - else - errno = 0; + iterable = args[0]; + if (!noptargs) { + goto skip_optional_kwonly; } - if (errno && is_error(r)) - return NULL; - else - return PyFloat_FromDouble(r); -} + start = args[1]; +skip_optional_kwonly: + return_value = math_prod_impl(module, iterable, start); -#define FUNC1(funcname, func, can_overflow, docstring) \ - static PyObject * math_##funcname(PyObject *self, PyObject *args) { \ - return math_1(args, func, can_overflow); \ - }\ - PyDoc_STRVAR(math_##funcname##_doc, docstring); - -#define FUNC1A(funcname, func, docstring) \ - static PyObject * math_##funcname(PyObject *self, PyObject *args) { \ - return math_1a(args, func); \ - }\ - PyDoc_STRVAR(math_##funcname##_doc, docstring); - -#define FUNC2(funcname, func, docstring) \ - static PyObject * math_##funcname(PyObject *self, PyObject *const *args, Py_ssize_t nargs) { \ - return math_2(args, nargs, func, #funcname); \ - }\ - PyDoc_STRVAR(math_##funcname##_doc, docstring); - -FUNC1(acos, acos, 0, - "acos($module, x, /)\n--\n\n" - "Return the arc cosine (measured in radians) of x.\n\n" - "The result is between 0 and pi.") -FUNC1(acosh, m_acosh, 0, - "acosh($module, x, /)\n--\n\n" - "Return the inverse hyperbolic cosine of x.") -FUNC1(asin, asin, 0, - "asin($module, x, /)\n--\n\n" - "Return the arc sine (measured in radians) of x.\n\n" - "The result is between -pi/2 and pi/2.") -FUNC1(asinh, m_asinh, 0, - "asinh($module, x, /)\n--\n\n" - "Return the inverse hyperbolic sine of x.") -FUNC1(atan, atan, 0, - "atan($module, x, /)\n--\n\n" - "Return the arc tangent (measured in radians) of x.\n\n" - "The result is between -pi/2 and pi/2.") -FUNC2(atan2, m_atan2, - "atan2($module, y, x, /)\n--\n\n" - "Return the arc tangent (measured in radians) of y/x.\n\n" - "Unlike atan(y/x), the signs of both x and y are considered.") -FUNC1(atanh, m_atanh, 0, - "atanh($module, x, /)\n--\n\n" - "Return the inverse hyperbolic tangent of x.") - -/*[clinic input] -math.ceil +exit: + return return_value; +} - x as number: object - / +PyDoc_STRVAR(math_perm__doc__, +"perm($module, n, k=None, /)\n" +"--\n" +"\n" +"Number of ways to choose k items from n items without repetition and with order.\n" +"\n" +"Evaluates to n! / (n - k)! when k <= n and evaluates\n" +"to zero when k > n.\n" +"\n" +"If k is not specified or is None, then k defaults to n\n" +"and the function returns n!.\n" +"\n" +"Raises TypeError if either of the arguments are not integers.\n" +"Raises ValueError if either of the arguments are negative."); -Return the ceiling of x as an Integral. +#define MATH_PERM_METHODDEF \ + {"perm", (PyCFunction)(void(*)(void))math_perm, METH_FASTCALL, math_perm__doc__}, -This is the smallest integer >= x. -[clinic start generated code]*/ +static PyObject * +math_perm_impl(PyObject *module, PyObject *n, PyObject *k); static PyObject * -math_ceil(PyObject *module, PyObject *number) -/*[clinic end generated code: output=6c3b8a78bc201c67 input=2725352806399cab]*/ +math_perm(PyObject *module, PyObject *const *args, Py_ssize_t nargs) { - _Py_IDENTIFIER(__ceil__); - - if (!PyFloat_CheckExact(number)) { - PyObject *method = _PyObject_LookupSpecial(number, &PyId___ceil__); - if (method != NULL) { - PyObject *result = _PyObject_CallNoArg(method); - Py_DECREF(method); - return result; - } - if (PyErr_Occurred()) - return NULL; + PyObject *return_value = NULL; + PyObject *n; + PyObject *k = Py_None; + + if (!_PyArg_CheckPositional("perm", nargs, 1, 2)) { + goto exit; + } + n = args[0]; + if (nargs < 2) { + goto skip_optional; } - double x = PyFloat_AsDouble(number); - if (x == -1.0 && PyErr_Occurred()) - return NULL; + k = args[1]; +skip_optional: + return_value = math_perm_impl(module, n, k); - return PyLong_FromDouble(ceil(x)); +exit: + return return_value; } -FUNC2(copysign, copysign, - "copysign($module, x, y, /)\n--\n\n" - "Return a float with the magnitude (absolute value) of x but the sign of y.\n\n" - "On platforms that support signed zeros, copysign(1.0, -0.0)\n" - "returns -1.0.\n") -FUNC1(cos, cos, 0, - "cos($module, x, /)\n--\n\n" - "Return the cosine of x (measured in radians).") -FUNC1(cosh, cosh, 1, - "cosh($module, x, /)\n--\n\n" - "Return the hyperbolic cosine of x.") -FUNC1A(erf, m_erf, - "erf($module, x, /)\n--\n\n" - "Error function at x.") -FUNC1A(erfc, m_erfc, - "erfc($module, x, /)\n--\n\n" - "Complementary error function at x.") -FUNC1(exp, exp, 1, - "exp($module, x, /)\n--\n\n" - "Return e raised to the power of x.") -FUNC1(expm1, m_expm1, 1, - "expm1($module, x, /)\n--\n\n" - "Return exp(x)-1.\n\n" - "This function avoids the loss of precision involved in the direct " - "evaluation of exp(x)-1 for small x.") -FUNC1(fabs, fabs, 0, - "fabs($module, x, /)\n--\n\n" - "Return the absolute value of the float x.") - -/*[clinic input] -math.floor - - x as number: object - / +PyDoc_STRVAR(math_comb__doc__, +"comb($module, n, k, /)\n" +"--\n" +"\n" +"Number of ways to choose k items from n items without repetition and without order.\n" +"\n" +"Evaluates to n! / (k! * (n - k)!) when k <= n and evaluates\n" +"to zero when k > n.\n" +"\n" +"Also called the binomial coefficient because it is equivalent\n" +"to the coefficient of k-th term in polynomial expansion of the\n" +"expression (1 + x)**n.\n" +"\n" +"Raises TypeError if either of the arguments are not integers.\n" +"Raises ValueError if either of the arguments are negative."); -Return the floor of x as an Integral. +#define MATH_COMB_METHODDEF \ + {"comb", (PyCFunction)(void(*)(void))math_comb, METH_FASTCALL, math_comb__doc__}, -This is the largest integer <= x. -[clinic start generated code]*/ +static PyObject * +math_comb_impl(PyObject *module, PyObject *n, PyObject *k); static PyObject * -math_floor(PyObject *module, PyObject *number) -/*[clinic end generated code: output=c6a65c4884884b8a input=63af6b5d7ebcc3d6]*/ +math_comb(PyObject *module, PyObject *const *args, Py_ssize_t nargs) { - _Py_IDENTIFIER(__floor__); - - if (!PyFloat_CheckExact(number)) { - PyObject *method = _PyObject_LookupSpecial(number, &PyId___floor__); - if (method != NULL) { - PyObject *result = _PyObject_CallNoArg(method); - Py_DECREF(method); - return result; - } - if (PyErr_Occurred()) - return NULL; + PyObject *return_value = NULL; + PyObject *n; + PyObject *k; + + if (!_PyArg_CheckPositional("comb", nargs, 2, 2)) { + goto exit; } - double x = PyFloat_AsDouble(number); - if (x == -1.0 && PyErr_Occurred()) - return NULL; + n = args[0]; + k = args[1]; + return_value = math_comb_impl(module, n, k); - return PyLong_FromDouble(floor(x)); +exit: + return return_value; } -FUNC1A(gamma, m_tgamma, - "gamma($module, x, /)\n--\n\n" - "Gamma function at x.") -FUNC1A(lgamma, m_lgamma, - "lgamma($module, x, /)\n--\n\n" - "Natural logarithm of absolute value of Gamma function at x.") -FUNC1(log1p, m_log1p, 0, - "log1p($module, x, /)\n--\n\n" - "Return the natural logarithm of 1+x (base e).\n\n" - "The result is computed in a way which is accurate for x near zero.") -FUNC2(remainder, m_remainder, - "remainder($module, x, y, /)\n--\n\n" - "Difference between x and the closest integer multiple of y.\n\n" - "Return x - n*y where n*y is the closest integer multiple of y.\n" - "In the case where x is exactly halfway between two multiples of\n" - "y, the nearest even value of n is used. The result is always exact.") -FUNC1(sin, sin, 0, - "sin($module, x, /)\n--\n\n" - "Return the sine of x (measured in radians).") -FUNC1(sinh, sinh, 1, - "sinh($module, x, /)\n--\n\n" - "Return the hyperbolic sine of x.") -FUNC1(sqrt, sqrt, 0, - "sqrt($module, x, /)\n--\n\n" - "Return the square root of x.") -FUNC1(tan, tan, 0, - "tan($module, x, /)\n--\n\n" - "Return the tangent of x (measured in radians).") -FUNC1(tanh, tanh, 0, - "tanh($module, x, /)\n--\n\n" - "Return the hyperbolic tangent of x.") - -/* Precision summation function as msum() by Raymond Hettinger in - , - enhanced with the exact partials sum and roundoff from Mark - Dickinson's post at . - See those links for more details, proofs and other references. - - Note 1: IEEE 754R floating point semantics are assumed, - but the current implementation does not re-establish special - value semantics across iterations (i.e. handling -Inf + Inf). - - Note 2: No provision is made for intermediate overflow handling; - therefore, sum([1e+308, 1e-308, 1e+308]) returns 1e+308 while - sum([1e+308, 1e+308, 1e-308]) raises an OverflowError due to the - overflow of the first partial sum. - - Note 3: The intermediate values lo, yr, and hi are declared volatile so - aggressive compilers won't algebraically reduce lo to always be exactly 0.0. - Also, the volatile declaration forces the values to be stored in memory as - regular doubles instead of extended long precision (80-bit) values. This - prevents double rounding because any addition or subtraction of two doubles - can be resolved exactly into double-sized hi and lo values. As long as the - hi value gets forced into a double before yr and lo are computed, the extra - bits in downstream extended precision operations (x87 for example) will be - exactly zero and therefore can be losslessly stored back into a double, - thereby preventing double rounding. - - Note 4: A similar implementation is in Modules/cmathmodule.c. - Be sure to update both when making changes. - - Note 5: The signature of math.fsum() differs from builtins.sum() - because the start argument doesn't make sense in the context of - accurate summation. Since the partials table is collapsed before - returning a result, sum(seq2, start=sum(seq1)) may not equal the - accurate result returned by sum(itertools.chain(seq1, seq2)). -*/ - -#define NUM_PARTIALS 32 /* initial partials array size, on stack */ - -/* Extend the partials array p[] by doubling its size. */ -static int /* non-zero on error */ -_fsum_realloc(double **p_ptr, Py_ssize_t n, - double *ps, Py_ssize_t *m_ptr) +PyDoc_STRVAR(math_nextafter__doc__, +"nextafter($module, x, y, /)\n" +"--\n" +"\n" +"Return the next floating-point value after x towards y."); + +#define MATH_NEXTAFTER_METHODDEF \ + {"nextafter", (PyCFunction)(void(*)(void))math_nextafter, METH_FASTCALL, math_nextafter__doc__}, + +static PyObject * +math_nextafter_impl(PyObject *module, double x, double y); + +static PyObject * +math_nextafter(PyObject *module, PyObject *const *args, Py_ssize_t nargs) { - void *v = NULL; - Py_ssize_t m = *m_ptr; - - m += m; /* double */ - if (n < m && (size_t)m < ((size_t)PY_SSIZE_T_MAX / sizeof(double))) { - double *p = *p_ptr; - if (p == ps) { - v = PyMem_Malloc(sizeof(double) * m); - if (v != NULL) - memcpy(v, ps, sizeof(double) * n); + PyObject *return_value = NULL; + double x; + double y; + + if (!_PyArg_CheckPositional("nextafter", nargs, 2, 2)) { + goto exit; + } + if (PyFloat_CheckExact(args[0])) { + x = PyFloat_AS_DOUBLE(args[0]); + } + else + { + x = PyFloat_AsDouble(args[0]); + if (x == -1.0 && PyErr_Occurred()) { + goto exit; } - else - v = PyMem_Realloc(p, sizeof(double) * m); } - if (v == NULL) { /* size overflow or no memory */ - PyErr_SetString(PyExc_MemoryError, "math.fsum partials"); - return 1; + if (PyFloat_CheckExact(args[1])) { + y = PyFloat_AS_DOUBLE(args[1]); } - *p_ptr = (double*) v; - *m_ptr = m; - return 0; -} - -/* Full precision summation of a sequence of floats. - - def msum(iterable): - partials = [] # sorted, non-overlapping partial sums - for x in iterable: - i = 0 - for y in partials: - if abs(x) < abs(y): - x, y = y, x - hi = x + y - lo = y - (hi - x) - if lo: - partials[i] = lo - i += 1 - x = hi - partials[i:] = [x] - return sum_exact(partials) - - Rounded x+y stored in hi with the roundoff stored in lo. Together hi+lo - are exactly equal to x+y. The inner loop applies hi/lo summation to each - partial so that the list of partial sums remains exact. - - Sum_exact() adds the partial sums exactly and correctly rounds the final - result (using the round-half-to-even rule). The items in partials remain - non-zero, non-special, non-overlapping and strictly increasing in - magnitude, but possibly not all having the same sign. - - Depends on IEEE 754 arithmetic guarantees and half-even rounding. -*/ + else + { + y = PyFloat_AsDouble(args[1]); + if (y == -1.0 && PyErr_Occurred()) { + goto exit; + } + } + return_value = math_nextafter_impl(module, x, y); -/*[clinic input] -math.fsum +exit: + return return_value; +} - seq: object - / +PyDoc_STRVAR(math_ulp__doc__, +"ulp($module, x, /)\n" +"--\n" +"\n" +"Return the value of the least significant bit of the float x."); -Return an accurate floating point sum of values in the iterable seq. +#define MATH_ULP_METHODDEF \ + {"ulp", (PyCFunction)math_ulp, METH_O, math_ulp__doc__}, -Assumes IEEE-754 floating point arithmetic. -[clinic start generated code]*/ +static double +math_ulp_impl(PyObject *module, double x); static PyObject * -math_fsum(PyObject *module, PyObject *seq) -/*[clinic end generated code: output=ba5c672b87fe34fc input=c51b7d8caf6f6e82]*/ +math_ulp(PyObject *module, PyObject *arg) { - PyObject *item, *iter, *sum = NULL; - Py_ssize_t i, j, n = 0, m = NUM_PARTIALS; - double x, y, t, ps[NUM_PARTIALS], *p = ps; - double xsave, special_sum = 0.0, inf_sum = 0.0; - volatile double hi, yr, lo; - - iter = PyObject_GetIter(seq); - if (iter == NULL) - return NULL; - - for(;;) { /* for x in iterable */ - assert(0 <= n && n <= m); - assert((m == NUM_PARTIALS && p == ps) || - (m > NUM_PARTIALS && p != NULL)); - - item = PyIter_Next(iter); - if (item == NULL) { - if (PyErr_Occurred()) - goto _fsum_error; - break; - } - ASSIGN_DOUBLE(x, item, error_with_item); - Py_DECREF(item); - - xsave = x; - for (i = j = 0; j < n; j++) { /* for y in partials */ - y = p[j]; - if (fabs(x) < fabs(y)) { - t = x; x = y; y = t; - } - hi = x + y; - yr = hi - x; - lo = y - yr; - if (lo != 0.0) - p[i++] = lo; - x = hi; - } - - n = i; /* ps[i:] = [x] */ - if (x != 0.0) { - if (! Py_IS_FINITE(x)) { - /* a nonfinite x could arise either as - a result of intermediate overflow, or - as a result of a nan or inf in the - summands */ - if (Py_IS_FINITE(xsave)) { - PyErr_SetString(PyExc_OverflowError, - "intermediate overflow in fsum"); - goto _fsum_error; - } - if (Py_IS_INFINITY(xsave)) - inf_sum += xsave; - special_sum += xsave; - /* reset partials */ - n = 0; - } - else if (n >= m && _fsum_realloc(&p, n, ps, &m)) - goto _fsum_error; - else - p[n++] = x; - } - } - - if (special_sum != 0.0) { - if (Py_IS_NAN(inf_sum)) - PyErr_SetString(PyExc_ValueError, - "-inf + inf in fsum"); - else - sum = PyFloat_FromDouble(special_sum); - goto _fsum_error; - } - - hi = 0.0; - if (n > 0) { - hi = p[--n]; - /* sum_exact(ps, hi) from the top, stop when the sum becomes - inexact. */ - while (n > 0) { - x = hi; - y = p[--n]; - assert(fabs(y) < fabs(x)); - hi = x + y; - yr = hi - x; - lo = y - yr; - if (lo != 0.0) - break; - } - /* Make half-even rounding work across multiple partials. - Needed so that sum([1e-16, 1, 1e16]) will round-up the last - digit to two instead of down to zero (the 1e-16 makes the 1 - slightly closer to two). With a potential 1 ULP rounding - error fixed-up, math.fsum() can guarantee commutativity. */ - if (n > 0 && ((lo < 0.0 && p[n-1] < 0.0) || - (lo > 0.0 && p[n-1] > 0.0))) { - y = lo * 2.0; - x = hi + y; - yr = x - hi; - if (y == yr) - hi = x; - } - } - sum = PyFloat_FromDouble(hi); - - _fsum_error: - Py_DECREF(iter); - if (p != ps) - PyMem_Free(p); - return sum; - - error_with_item: - Py_DECREF(item); - goto _fsum_error; -} - -#undef NUM_PARTIALS - - -static unsigned long -count_set_bits(unsigned long n) -{ - unsigned long count = 0; - while (n != 0) { - ++count; - n &= n - 1; /* clear least significant bit */ - } - return count; -} - -/* Integer square root - -Given a nonnegative integer `n`, we want to compute the largest integer -`a` for which `a * a <= n`, or equivalently the integer part of the exact -square root of `n`. - -We use an adaptive-precision pure-integer version of Newton's iteration. Given -a positive integer `n`, the algorithm produces at each iteration an integer -approximation `a` to the square root of `n >> s` for some even integer `s`, -with `s` decreasing as the iterations progress. On the final iteration, `s` is -zero and we have an approximation to the square root of `n` itself. - -At every step, the approximation `a` is strictly within 1.0 of the true square -root, so we have - - (a - 1)**2 < (n >> s) < (a + 1)**2 - -After the final iteration, a check-and-correct step is needed to determine -whether `a` or `a - 1` gives the desired integer square root of `n`. - -The algorithm is remarkable in its simplicity. There's no need for a -per-iteration check-and-correct step, and termination is straightforward: the -number of iterations is known in advance (it's exactly `floor(log2(log2(n)))` -for `n > 1`). The only tricky part of the correctness proof is in establishing -that the bound `(a - 1)**2 < (n >> s) < (a + 1)**2` is maintained from one -iteration to the next. A sketch of the proof of this is given below. - -In addition to the proof sketch, a formal, computer-verified proof -of correctness (using Lean) of an equivalent recursive algorithm can be found -here: - - https://github.com/mdickinson/snippets/blob/master/proofs/isqrt/src/isqrt.lean - - -Here's Python code equivalent to the C implementation below: - - def isqrt(n): - """ - Return the integer part of the square root of the input. - """ - n = operator.index(n) - - if n < 0: - raise ValueError("isqrt() argument must be nonnegative") - if n == 0: - return 0 - - c = (n.bit_length() - 1) // 2 - a = 1 - d = 0 - for s in reversed(range(c.bit_length())): - # Loop invariant: (a-1)**2 < (n >> 2*(c - d)) < (a+1)**2 - e = d - d = c >> s - a = (a << d - e - 1) + (n >> 2*c - e - d + 1) // a - - return a - (a*a > n) - - -Sketch of proof of correctness ------------------------------- - -The delicate part of the correctness proof is showing that the loop invariant -is preserved from one iteration to the next. That is, just before the line - - a = (a << d - e - 1) + (n >> 2*c - e - d + 1) // a - -is executed in the above code, we know that - - (1) (a - 1)**2 < (n >> 2*(c - e)) < (a + 1)**2. - -(since `e` is always the value of `d` from the previous iteration). We must -prove that after that line is executed, we have - - (a - 1)**2 < (n >> 2*(c - d)) < (a + 1)**2 - -To facilitate the proof, we make some changes of notation. Write `m` for -`n >> 2*(c-d)`, and write `b` for the new value of `a`, so - - b = (a << d - e - 1) + (n >> 2*c - e - d + 1) // a - -or equivalently: - - (2) b = (a << d - e - 1) + (m >> d - e + 1) // a - -Then we can rewrite (1) as: - - (3) (a - 1)**2 < (m >> 2*(d - e)) < (a + 1)**2 - -and we must show that (b - 1)**2 < m < (b + 1)**2. - -From this point on, we switch to mathematical notation, so `/` means exact -division rather than integer division and `^` is used for exponentiation. We -use the `√` symbol for the exact square root. In (3), we can remove the -implicit floor operation to give: - - (4) (a - 1)^2 < m / 4^(d - e) < (a + 1)^2 - -Taking square roots throughout (4), scaling by `2^(d-e)`, and rearranging gives - - (5) 0 <= | 2^(d-e)a - √m | < 2^(d-e) - -Squaring and dividing through by `2^(d-e+1) a` gives - - (6) 0 <= 2^(d-e-1) a + m / (2^(d-e+1) a) - √m < 2^(d-e-1) / a - -We'll show below that `2^(d-e-1) <= a`. Given that, we can replace the -right-hand side of (6) with `1`, and now replacing the central -term `m / (2^(d-e+1) a)` with its floor in (6) gives - - (7) -1 < 2^(d-e-1) a + m // 2^(d-e+1) a - √m < 1 - -Or equivalently, from (2): - - (7) -1 < b - √m < 1 - -and rearranging gives that `(b-1)^2 < m < (b+1)^2`, which is what we needed -to prove. - -We're not quite done: we still have to prove the inequality `2^(d - e - 1) <= -a` that was used to get line (7) above. From the definition of `c`, we have -`4^c <= n`, which implies - - (8) 4^d <= m - -also, since `e == d >> 1`, `d` is at most `2e + 1`, from which it follows -that `2d - 2e - 1 <= d` and hence that - - (9) 4^(2d - 2e - 1) <= m - -Dividing both sides by `4^(d - e)` gives - - (10) 4^(d - e - 1) <= m / 4^(d - e) - -But we know from (4) that `m / 4^(d-e) < (a + 1)^2`, hence - - (11) 4^(d - e - 1) < (a + 1)^2 - -Now taking square roots of both sides and observing that both `2^(d-e-1)` and -`a` are integers gives `2^(d - e - 1) <= a`, which is what we needed. This -completes the proof sketch. - -*/ - - -/* Approximate square root of a large 64-bit integer. - - Given `n` satisfying `2**62 <= n < 2**64`, return `a` - satisfying `(a - 1)**2 < n < (a + 1)**2`. */ - -static uint64_t -_approximate_isqrt(uint64_t n) -{ - uint32_t u = 1U + (n >> 62); - u = (u << 1) + (n >> 59) / u; - u = (u << 3) + (n >> 53) / u; - u = (u << 7) + (n >> 41) / u; - return (u << 15) + (n >> 17) / u; -} - -/*[clinic input] -math.isqrt - - n: object - / - -Return the integer part of the square root of the input. -[clinic start generated code]*/ - -static PyObject * -math_isqrt(PyObject *module, PyObject *n) -/*[clinic end generated code: output=35a6f7f980beab26 input=5b6e7ae4fa6c43d6]*/ -{ - int a_too_large, c_bit_length; - size_t c, d; - uint64_t m, u; - PyObject *a = NULL, *b; - - n = PyNumber_Index(n); - if (n == NULL) { - return NULL; - } - - if (_PyLong_Sign(n) < 0) { - PyErr_SetString( - PyExc_ValueError, - "isqrt() argument must be nonnegative"); - goto error; - } - if (_PyLong_Sign(n) == 0) { - Py_DECREF(n); - return PyLong_FromLong(0); - } - - /* c = (n.bit_length() - 1) // 2 */ - c = _PyLong_NumBits(n); - if (c == (size_t)(-1)) { - goto error; - } - c = (c - 1U) / 2U; - - /* Fast path: if c <= 31 then n < 2**64 and we can compute directly with a - fast, almost branch-free algorithm. In the final correction, we use `u*u - - 1 >= m` instead of the simpler `u*u > m` in order to get the correct - result in the corner case where `u=2**32`. */ - if (c <= 31U) { - m = (uint64_t)PyLong_AsUnsignedLongLong(n); - Py_DECREF(n); - if (m == (uint64_t)(-1) && PyErr_Occurred()) { - return NULL; - } - u = _approximate_isqrt(m << (62U - 2U*c)) >> (31U - c); - u -= u * u - 1U >= m; - return PyLong_FromUnsignedLongLong((unsigned long long)u); - } - - /* Slow path: n >= 2**64. We perform the first five iterations in C integer - arithmetic, then switch to using Python long integers. */ - - /* From n >= 2**64 it follows that c.bit_length() >= 6. */ - c_bit_length = 6; - while ((c >> c_bit_length) > 0U) { - ++c_bit_length; - } - - /* Initialise d and a. */ - d = c >> (c_bit_length - 5); - b = _PyLong_Rshift(n, 2U*c - 62U); - if (b == NULL) { - goto error; - } - m = (uint64_t)PyLong_AsUnsignedLongLong(b); - Py_DECREF(b); - if (m == (uint64_t)(-1) && PyErr_Occurred()) { - goto error; - } - u = _approximate_isqrt(m) >> (31U - d); - a = PyLong_FromUnsignedLongLong((unsigned long long)u); - if (a == NULL) { - goto error; - } - - for (int s = c_bit_length - 6; s >= 0; --s) { - PyObject *q; - size_t e = d; - - d = c >> s; - - /* q = (n >> 2*c - e - d + 1) // a */ - q = _PyLong_Rshift(n, 2U*c - d - e + 1U); - if (q == NULL) { - goto error; - } - Py_SETREF(q, PyNumber_FloorDivide(q, a)); - if (q == NULL) { - goto error; - } - - /* a = (a << d - 1 - e) + q */ - Py_SETREF(a, _PyLong_Lshift(a, d - 1U - e)); - if (a == NULL) { - Py_DECREF(q); - goto error; - } - Py_SETREF(a, PyNumber_Add(a, q)); - Py_DECREF(q); - if (a == NULL) { - goto error; - } - } - - /* The correct result is either a or a - 1. Figure out which, and - decrement a if necessary. */ - - /* a_too_large = n < a * a */ - b = PyNumber_Multiply(a, a); - if (b == NULL) { - goto error; - } - a_too_large = PyObject_RichCompareBool(n, b, Py_LT); - Py_DECREF(b); - if (a_too_large == -1) { - goto error; - } - - if (a_too_large) { - Py_SETREF(a, PyNumber_Subtract(a, _PyLong_One)); - } - Py_DECREF(n); - return a; - - error: - Py_XDECREF(a); - Py_DECREF(n); - return NULL; -} - -/* Divide-and-conquer factorial algorithm - * - * Based on the formula and pseudo-code provided at: - * http://www.luschny.de/math/factorial/binarysplitfact.html - * - * Faster algorithms exist, but they're more complicated and depend on - * a fast prime factorization algorithm. - * - * Notes on the algorithm - * ---------------------- - * - * factorial(n) is written in the form 2**k * m, with m odd. k and m are - * computed separately, and then combined using a left shift. - * - * The function factorial_odd_part computes the odd part m (i.e., the greatest - * odd divisor) of factorial(n), using the formula: - * - * factorial_odd_part(n) = - * - * product_{i >= 0} product_{0 < j <= n / 2**i, j odd} j - * - * Example: factorial_odd_part(20) = - * - * (1) * - * (1) * - * (1 * 3 * 5) * - * (1 * 3 * 5 * 7 * 9) - * (1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19) - * - * Here i goes from large to small: the first term corresponds to i=4 (any - * larger i gives an empty product), and the last term corresponds to i=0. - * Each term can be computed from the last by multiplying by the extra odd - * numbers required: e.g., to get from the penultimate term to the last one, - * we multiply by (11 * 13 * 15 * 17 * 19). - * - * To see a hint of why this formula works, here are the same numbers as above - * but with the even parts (i.e., the appropriate powers of 2) included. For - * each subterm in the product for i, we multiply that subterm by 2**i: - * - * factorial(20) = - * - * (16) * - * (8) * - * (4 * 12 * 20) * - * (2 * 6 * 10 * 14 * 18) * - * (1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19) - * - * The factorial_partial_product function computes the product of all odd j in - * range(start, stop) for given start and stop. It's used to compute the - * partial products like (11 * 13 * 15 * 17 * 19) in the example above. It - * operates recursively, repeatedly splitting the range into two roughly equal - * pieces until the subranges are small enough to be computed using only C - * integer arithmetic. - * - * The two-valuation k (i.e., the exponent of the largest power of 2 dividing - * the factorial) is computed independently in the main math_factorial - * function. By standard results, its value is: - * - * two_valuation = n//2 + n//4 + n//8 + .... - * - * It can be shown (e.g., by complete induction on n) that two_valuation is - * equal to n - count_set_bits(n), where count_set_bits(n) gives the number of - * '1'-bits in the binary expansion of n. - */ - -/* factorial_partial_product: Compute product(range(start, stop, 2)) using - * divide and conquer. Assumes start and stop are odd and stop > start. - * max_bits must be >= bit_length(stop - 2). */ - -static PyObject * -factorial_partial_product(unsigned long start, unsigned long stop, - unsigned long max_bits) -{ - unsigned long midpoint, num_operands; - PyObject *left = NULL, *right = NULL, *result = NULL; - - /* If the return value will fit an unsigned long, then we can - * multiply in a tight, fast loop where each multiply is O(1). - * Compute an upper bound on the number of bits required to store - * the answer. - * - * Storing some integer z requires floor(lg(z))+1 bits, which is - * conveniently the value returned by bit_length(z). The - * product x*y will require at most - * bit_length(x) + bit_length(y) bits to store, based - * on the idea that lg product = lg x + lg y. - * - * We know that stop - 2 is the largest number to be multiplied. From - * there, we have: bit_length(answer) <= num_operands * - * bit_length(stop - 2) - */ - - num_operands = (stop - start) / 2; - /* The "num_operands <= 8 * SIZEOF_LONG" check guards against the - * unlikely case of an overflow in num_operands * max_bits. */ - if (num_operands <= 8 * SIZEOF_LONG && - num_operands * max_bits <= 8 * SIZEOF_LONG) { - unsigned long j, total; - for (total = start, j = start + 2; j < stop; j += 2) - total *= j; - return PyLong_FromUnsignedLong(total); - } - - /* find midpoint of range(start, stop), rounded up to next odd number. */ - midpoint = (start + num_operands) | 1; - left = factorial_partial_product(start, midpoint, - _Py_bit_length(midpoint - 2)); - if (left == NULL) - goto error; - right = factorial_partial_product(midpoint, stop, max_bits); - if (right == NULL) - goto error; - result = PyNumber_Multiply(left, right); - - error: - Py_XDECREF(left); - Py_XDECREF(right); - return result; -} - -/* factorial_odd_part: compute the odd part of factorial(n). */ - -static PyObject * -factorial_odd_part(unsigned long n) -{ - long i; - unsigned long v, lower, upper; - PyObject *partial, *tmp, *inner, *outer; - - inner = PyLong_FromLong(1); - if (inner == NULL) - return NULL; - outer = inner; - Py_INCREF(outer); - - upper = 3; - for (i = _Py_bit_length(n) - 2; i >= 0; i--) { - v = n >> i; - if (v <= 2) - continue; - lower = upper; - /* (v + 1) | 1 = least odd integer strictly larger than n / 2**i */ - upper = (v + 1) | 1; - /* Here inner is the product of all odd integers j in the range (0, - n/2**(i+1)]. The factorial_partial_product call below gives the - product of all odd integers j in the range (n/2**(i+1), n/2**i]. */ - partial = factorial_partial_product(lower, upper, _Py_bit_length(upper-2)); - /* inner *= partial */ - if (partial == NULL) - goto error; - tmp = PyNumber_Multiply(inner, partial); - Py_DECREF(partial); - if (tmp == NULL) - goto error; - Py_DECREF(inner); - inner = tmp; - /* Now inner is the product of all odd integers j in the range (0, - n/2**i], giving the inner product in the formula above. */ - - /* outer *= inner; */ - tmp = PyNumber_Multiply(outer, inner); - if (tmp == NULL) - goto error; - Py_DECREF(outer); - outer = tmp; - } - Py_DECREF(inner); - return outer; - - error: - Py_DECREF(outer); - Py_DECREF(inner); - return NULL; -} - - -/* Lookup table for small factorial values */ - -static const unsigned long SmallFactorials[] = { - 1, 1, 2, 6, 24, 120, 720, 5040, 40320, - 362880, 3628800, 39916800, 479001600, -#if SIZEOF_LONG >= 8 - 6227020800, 87178291200, 1307674368000, - 20922789888000, 355687428096000, 6402373705728000, - 121645100408832000, 2432902008176640000 -#endif -}; - -/*[clinic input] -math.factorial - - x as arg: object - / - -Find x!. - -Raise a ValueError if x is negative or non-integral. -[clinic start generated code]*/ - -static PyObject * -math_factorial(PyObject *module, PyObject *arg) -/*[clinic end generated code: output=6686f26fae00e9ca input=6d1c8105c0d91fb4]*/ -{ - long x, two_valuation; - int overflow; - PyObject *result, *odd_part, *pyint_form; - - if (PyFloat_Check(arg)) { - if (PyErr_WarnEx(PyExc_DeprecationWarning, - "Using factorial() with floats is deprecated", - 1) < 0) - { - return NULL; - } - PyObject *lx; - double dx = PyFloat_AS_DOUBLE((PyFloatObject *)arg); - if (!(Py_IS_FINITE(dx) && dx == floor(dx))) { - PyErr_SetString(PyExc_ValueError, - "factorial() only accepts integral values"); - return NULL; - } - lx = PyLong_FromDouble(dx); - if (lx == NULL) - return NULL; - x = PyLong_AsLongAndOverflow(lx, &overflow); - Py_DECREF(lx); - } - else { - pyint_form = PyNumber_Index(arg); - if (pyint_form == NULL) { - return NULL; - } - x = PyLong_AsLongAndOverflow(pyint_form, &overflow); - Py_DECREF(pyint_form); - } - - if (x == -1 && PyErr_Occurred()) { - return NULL; - } - else if (overflow == 1) { - PyErr_Format(PyExc_OverflowError, - "factorial() argument should not exceed %ld", - LONG_MAX); - return NULL; - } - else if (overflow == -1 || x < 0) { - PyErr_SetString(PyExc_ValueError, - "factorial() not defined for negative values"); - return NULL; - } - - /* use lookup table if x is small */ - if (x < (long)Py_ARRAY_LENGTH(SmallFactorials)) - return PyLong_FromUnsignedLong(SmallFactorials[x]); - - /* else express in the form odd_part * 2**two_valuation, and compute as - odd_part << two_valuation. */ - odd_part = factorial_odd_part(x); - if (odd_part == NULL) - return NULL; - two_valuation = x - count_set_bits(x); - result = _PyLong_Lshift(odd_part, two_valuation); - Py_DECREF(odd_part); - return result; -} - - -/*[clinic input] -math.lcm - x as a: object - y as b: object - / -least common multiple of x and y -[clinic start generated code]*/ - -static PyObject * -math_lcm_impl(PyObject *module, PyObject *a, PyObject *b) -/*[clinic end generated code: output=6f83fb6d671074ba input=efb3d7b7334b7118]*/ - - PyObject *g, *m, *f, *ab; - - a = PyNumber_Index(a); - if (a == NULL) { - return NULL; - } - b = PyNumber_Index(b); - if (b == NULL) { - Py_DECREF(a); - return NULL; - } - if (_PyLong_Sign(a) == 0 || _PyLong_Sign(b) == 0) { - Py_DECREF(a); - Py_DECREF(b); - return PyLong_FromLong(0); - } - g = _PyLong_GCD(a, b); - if (g == NULL) { - Py_DECREF(a); - Py_DECREF(b); - return NULL; - } - f = PyNumber_FloorDivide(a, g); - Py_DECREF(g); - Py_DECREF(a); - if (f == NULL) { - Py_DECREF(b); - return NULL; - } - m = PyNumber_Multiply(f, b); - Py_DECREF(f); - Py_DECREF(b); - if (m == NULL) { - return NULL; - } - ab = PyNumber_Absolute(m); - Py_DECREF(m); - return ab; -} - - -/*[clinic input] -math.trunc - - x: object - / - -Truncates the Real x to the nearest Integral toward 0. - -Uses the __trunc__ magic method. -[clinic start generated code]*/ - -static PyObject * -math_trunc(PyObject *module, PyObject *x) -/*[clinic end generated code: output=34b9697b707e1031 input=2168b34e0a09134d]*/ -{ - _Py_IDENTIFIER(__trunc__); - PyObject *trunc, *result; - - if (PyFloat_CheckExact(x)) { - return PyFloat_Type.tp_as_number->nb_int(x); - } - - if (Py_TYPE(x)->tp_dict == NULL) { - if (PyType_Ready(Py_TYPE(x)) < 0) - return NULL; - } - - trunc = _PyObject_LookupSpecial(x, &PyId___trunc__); - if (trunc == NULL) { - if (!PyErr_Occurred()) - PyErr_Format(PyExc_TypeError, - "type %.100s doesn't define __trunc__ method", - Py_TYPE(x)->tp_name); - return NULL; - } - result = _PyObject_CallNoArg(trunc); - Py_DECREF(trunc); - return result; -} - - -/*[clinic input] -math.frexp - - x: double - / - -Return the mantissa and exponent of x, as pair (m, e). - -m is a float and e is an int, such that x = m * 2.**e. -If x is 0, m and e are both 0. Else 0.5 <= abs(m) < 1.0. -[clinic start generated code]*/ - -static PyObject * -math_frexp_impl(PyObject *module, double x) -/*[clinic end generated code: output=03e30d252a15ad4a input=96251c9e208bc6e9]*/ -{ - int i; - /* deal with special cases directly, to sidestep platform - differences */ - if (Py_IS_NAN(x) || Py_IS_INFINITY(x) || !x) { - i = 0; - } - else { - x = frexp(x, &i); - } - return Py_BuildValue("(di)", x, i); -} - - -/*[clinic input] -math.ldexp - - x: double - i: object - / - -Return x * (2**i). - -This is essentially the inverse of frexp(). -[clinic start generated code]*/ - -static PyObject * -math_ldexp_impl(PyObject *module, double x, PyObject *i) -/*[clinic end generated code: output=b6892f3c2df9cc6a input=17d5970c1a40a8c1]*/ -{ - double r; - long exp; - int overflow; - - if (PyLong_Check(i)) { - /* on overflow, replace exponent with either LONG_MAX - or LONG_MIN, depending on the sign. */ - exp = PyLong_AsLongAndOverflow(i, &overflow); - if (exp == -1 && PyErr_Occurred()) - return NULL; - if (overflow) - exp = overflow < 0 ? LONG_MIN : LONG_MAX; - } - else { - PyErr_SetString(PyExc_TypeError, - "Expected an int as second argument to ldexp."); - return NULL; - } - - if (x == 0. || !Py_IS_FINITE(x)) { - /* NaNs, zeros and infinities are returned unchanged */ - r = x; - errno = 0; - } else if (exp > INT_MAX) { - /* overflow */ - r = copysign(Py_HUGE_VAL, x); - errno = ERANGE; - } else if (exp < INT_MIN) { - /* underflow to +-0 */ - r = copysign(0., x); - errno = 0; - } else { - errno = 0; - r = ldexp(x, (int)exp); - if (Py_IS_INFINITY(r)) - errno = ERANGE; - } - - if (errno && is_error(r)) - return NULL; - return PyFloat_FromDouble(r); -} - - -/*[clinic input] -math.modf - - x: double - / - -Return the fractional and integer parts of x. - -Both results carry the sign of x and are floats. -[clinic start generated code]*/ - -static PyObject * -math_modf_impl(PyObject *module, double x) -/*[clinic end generated code: output=90cee0260014c3c0 input=b4cfb6786afd9035]*/ -{ - double y; - /* some platforms don't do the right thing for NaNs and - infinities, so we take care of special cases directly. */ - if (!Py_IS_FINITE(x)) { - if (Py_IS_INFINITY(x)) - return Py_BuildValue("(dd)", copysign(0., x), x); - else if (Py_IS_NAN(x)) - return Py_BuildValue("(dd)", x, x); - } - - errno = 0; - x = modf(x, &y); - return Py_BuildValue("(dd)", x, y); -} - - -/* A decent logarithm is easy to compute even for huge ints, but libm can't - do that by itself -- loghelper can. func is log or log10, and name is - "log" or "log10". Note that overflow of the result isn't possible: an int - can contain no more than INT_MAX * SHIFT bits, so has value certainly less - than 2**(2**64 * 2**16) == 2**2**80, and log2 of that is 2**80, which is - small enough to fit in an IEEE single. log and log10 are even smaller. - However, intermediate overflow is possible for an int if the number of bits - in that int is larger than PY_SSIZE_T_MAX. */ - -static PyObject* -loghelper(PyObject* arg, double (*func)(double), const char *funcname) -{ - /* If it is int, do it ourselves. */ - if (PyLong_Check(arg)) { - double x, result; - Py_ssize_t e; - - /* Negative or zero inputs give a ValueError. */ - if (Py_SIZE(arg) <= 0) { - PyErr_SetString(PyExc_ValueError, - "math domain error"); - return NULL; - } - - x = PyLong_AsDouble(arg); - if (x == -1.0 && PyErr_Occurred()) { - if (!PyErr_ExceptionMatches(PyExc_OverflowError)) - return NULL; - /* Here the conversion to double overflowed, but it's possible - to compute the log anyway. Clear the exception and continue. */ - PyErr_Clear(); - x = _PyLong_Frexp((PyLongObject *)arg, &e); - if (x == -1.0 && PyErr_Occurred()) - return NULL; - /* Value is ~= x * 2**e, so the log ~= log(x) + log(2) * e. */ - result = func(x) + func(2.0) * e; - } - else - /* Successfully converted x to a double. */ - result = func(x); - return PyFloat_FromDouble(result); - } - - /* Else let libm handle it by itself. */ - return math_1(arg, func, 0); -} - - -/*[clinic input] -math.log - - x: object - [ - base: object(c_default="NULL") = math.e - ] - / - -Return the logarithm of x to the given base. - -If the base not specified, returns the natural logarithm (base e) of x. -[clinic start generated code]*/ - -static PyObject * -math_log_impl(PyObject *module, PyObject *x, int group_right_1, - PyObject *base) -/*[clinic end generated code: output=7b5a39e526b73fc9 input=0f62d5726cbfebbd]*/ -{ - PyObject *num, *den; - PyObject *ans; - - num = loghelper(x, m_log, "log"); - if (num == NULL || base == NULL) - return num; - - den = loghelper(base, m_log, "log"); - if (den == NULL) { - Py_DECREF(num); - return NULL; - } - - ans = PyNumber_TrueDivide(num, den); - Py_DECREF(num); - Py_DECREF(den); - return ans; -} - - -/*[clinic input] -math.log2 - - x: object - / - -Return the base 2 logarithm of x. -[clinic start generated code]*/ - -static PyObject * -math_log2(PyObject *module, PyObject *x) -/*[clinic end generated code: output=5425899a4d5d6acb input=08321262bae4f39b]*/ -{ - return loghelper(x, m_log2, "log2"); -} - - -/*[clinic input] -math.log10 - - x: object - / - -Return the base 10 logarithm of x. -[clinic start generated code]*/ - -static PyObject * -math_log10(PyObject *module, PyObject *x) -/*[clinic end generated code: output=be72a64617df9c6f input=b2469d02c6469e53]*/ -{ - return loghelper(x, m_log10, "log10"); -} - - -/*[clinic input] -math.fmod - - x: double - y: double - / - -Return fmod(x, y), according to platform C. - -x % y may differ. -[clinic start generated code]*/ - -static PyObject * -math_fmod_impl(PyObject *module, double x, double y) -/*[clinic end generated code: output=7559d794343a27b5 input=4f84caa8cfc26a03]*/ -{ - double r; - /* fmod(x, +/-Inf) returns x for finite x. */ - if (Py_IS_INFINITY(y) && Py_IS_FINITE(x)) - return PyFloat_FromDouble(x); - errno = 0; - r = fmod(x, y); - if (Py_IS_NAN(r)) { - if (!Py_IS_NAN(x) && !Py_IS_NAN(y)) - errno = EDOM; - else - errno = 0; - } - if (errno && is_error(r)) - return NULL; - else - return PyFloat_FromDouble(r); -} - -/* -Given an *n* length *vec* of values and a value *max*, compute: - - max * sqrt(sum((x / max) ** 2 for x in vec)) - -The value of the *max* variable must be non-negative and -equal to the absolute value of the largest magnitude -entry in the vector. If n==0, then *max* should be 0.0. -If an infinity is present in the vec, *max* should be INF. - -The *found_nan* variable indicates whether some member of -the *vec* is a NaN. - -To improve accuracy and to increase the number of cases where -vector_norm() is commutative, we use a variant of Neumaier -summation specialized to exploit that we always know that -|csum| >= |x|. - -The *csum* variable tracks the cumulative sum and *frac* tracks -the cumulative fractional errors at each step. Since this -variant assumes that |csum| >= |x| at each step, we establish -the precondition by starting the accumulation from 1.0 which -represents the largest possible value of (x/max)**2. - -After the loop is finished, the initial 1.0 is subtracted out -for a net zero effect on the final sum. Since *csum* will be -greater than 1.0, the subtraction of 1.0 will not cause -fractional digits to be dropped from *csum*. - -*/ - -static inline double -vector_norm(Py_ssize_t n, double *vec, double max, int found_nan) -{ - double x, csum = 1.0, oldcsum, frac = 0.0; - Py_ssize_t i; - - if (Py_IS_INFINITY(max)) { - return max; - } - if (found_nan) { - return Py_NAN; - } - if (max == 0.0 || n <= 1) { - return max; - } - for (i=0 ; i < n ; i++) { - x = vec[i]; - assert(Py_IS_FINITE(x) && fabs(x) <= max); - x /= max; - x = x*x; - oldcsum = csum; - csum += x; - assert(csum >= x); - frac += (oldcsum - csum) + x; - } - return max * sqrt(csum - 1.0 + frac); -} - -#define NUM_STACK_ELEMS 16 - -/*[clinic input] -math.dist - - p: object - q: object - / - -Return the Euclidean distance between two points p and q. - -The points should be specified as sequences (or iterables) of -coordinates. Both inputs must have the same dimension. - -Roughly equivalent to: - sqrt(sum((px - qx) ** 2.0 for px, qx in zip(p, q))) -[clinic start generated code]*/ - -static PyObject * -math_dist_impl(PyObject *module, PyObject *p, PyObject *q) -/*[clinic end generated code: output=56bd9538d06bbcfe input=74e85e1b6092e68e]*/ -{ - PyObject *item; - double max = 0.0; - double x, px, qx, result; - Py_ssize_t i, m, n; - int found_nan = 0, p_allocated = 0, q_allocated = 0; - double diffs_on_stack[NUM_STACK_ELEMS]; - double *diffs = diffs_on_stack; - - if (!PyTuple_Check(p)) { - p = PySequence_Tuple(p); - if (p == NULL) { - return NULL; - } - p_allocated = 1; - } - if (!PyTuple_Check(q)) { - q = PySequence_Tuple(q); - if (q == NULL) { - if (p_allocated) { - Py_DECREF(p); - } - return NULL; - } - q_allocated = 1; - } - - m = PyTuple_GET_SIZE(p); - n = PyTuple_GET_SIZE(q); - if (m != n) { - PyErr_SetString(PyExc_ValueError, - "both points must have the same number of dimensions"); - return NULL; - - } - if (n > NUM_STACK_ELEMS) { - diffs = (double *) PyObject_Malloc(n * sizeof(double)); - if (diffs == NULL) { - return PyErr_NoMemory(); - } - } - for (i=0 ; i max) { - max = x; - } - } - result = vector_norm(n, diffs, max, found_nan); - if (diffs != diffs_on_stack) { - PyObject_Free(diffs); - } - if (p_allocated) { - Py_DECREF(p); - } - if (q_allocated) { - Py_DECREF(q); - } - return PyFloat_FromDouble(result); - - error_exit: - if (diffs != diffs_on_stack) { - PyObject_Free(diffs); - } - if (p_allocated) { - Py_DECREF(p); - } - if (q_allocated) { - Py_DECREF(q); - } - return NULL; -} - -/* AC: cannot convert yet, waiting for *args support */ -static PyObject * -math_hypot(PyObject *self, PyObject *const *args, Py_ssize_t nargs) -{ - Py_ssize_t i; - PyObject *item; - double max = 0.0; - double x, result; - int found_nan = 0; - double coord_on_stack[NUM_STACK_ELEMS]; - double *coordinates = coord_on_stack; - - if (nargs > NUM_STACK_ELEMS) { - coordinates = (double *) PyObject_Malloc(nargs * sizeof(double)); - if (coordinates == NULL) { - return PyErr_NoMemory(); - } - } - for (i = 0; i < nargs; i++) { - item = args[i]; - ASSIGN_DOUBLE(x, item, error_exit); - x = fabs(x); - coordinates[i] = x; - found_nan |= Py_IS_NAN(x); - if (x > max) { - max = x; - } - } - result = vector_norm(nargs, coordinates, max, found_nan); - if (coordinates != coord_on_stack) { - PyObject_Free(coordinates); - } - return PyFloat_FromDouble(result); - - error_exit: - if (coordinates != coord_on_stack) { - PyObject_Free(coordinates); - } - return NULL; -} - -#undef NUM_STACK_ELEMS - -PyDoc_STRVAR(math_hypot_doc, - "hypot(*coordinates) -> value\n\n\ -Multidimensional Euclidean distance from the origin to a point.\n\ -\n\ -Roughly equivalent to:\n\ - sqrt(sum(x**2 for x in coordinates))\n\ -\n\ -For a two dimensional point (x, y), gives the hypotenuse\n\ -using the Pythagorean theorem: sqrt(x*x + y*y).\n\ -\n\ -For example, the hypotenuse of a 3/4/5 right triangle is:\n\ -\n\ - >>> hypot(3.0, 4.0)\n\ - 5.0\n\ -"); - -/* pow can't use math_2, but needs its own wrapper: the problem is - that an infinite result can arise either as a result of overflow - (in which case OverflowError should be raised) or as a result of - e.g. 0.**-5. (for which ValueError needs to be raised.) -*/ - -/*[clinic input] -math.pow - - x: double - y: double - / - -Return x**y (x to the power of y). -[clinic start generated code]*/ + PyObject *return_value = NULL; + double x; + double _return_value; -static PyObject * -math_pow_impl(PyObject *module, double x, double y) -/*[clinic end generated code: output=fff93e65abccd6b0 input=c26f1f6075088bfd]*/ -{ - double r; - int odd_y; - - /* deal directly with IEEE specials, to cope with problems on various - platforms whose semantics don't exactly match C99 */ - r = 0.; /* silence compiler warning */ - if (!Py_IS_FINITE(x) || !Py_IS_FINITE(y)) { - errno = 0; - if (Py_IS_NAN(x)) - r = y == 0. ? 1. : x; /* NaN**0 = 1 */ - else if (Py_IS_NAN(y)) - r = x == 1. ? 1. : y; /* 1**NaN = 1 */ - else if (Py_IS_INFINITY(x)) { - odd_y = Py_IS_FINITE(y) && fmod(fabs(y), 2.0) == 1.0; - if (y > 0.) - r = odd_y ? x : fabs(x); - else if (y == 0.) - r = 1.; - else /* y < 0. */ - r = odd_y ? copysign(0., x) : 0.; - } - else if (Py_IS_INFINITY(y)) { - if (fabs(x) == 1.0) - r = 1.; - else if (y > 0. && fabs(x) > 1.0) - r = y; - else if (y < 0. && fabs(x) < 1.0) { - r = -y; /* result is +inf */ - if (x == 0.) /* 0**-inf: divide-by-zero */ - errno = EDOM; - } - else - r = 0.; - } - } - else { - /* let libm handle finite**finite */ - errno = 0; - r = pow(x, y); - /* a NaN result should arise only from (-ve)**(finite - non-integer); in this case we want to raise ValueError. */ - if (!Py_IS_FINITE(r)) { - if (Py_IS_NAN(r)) { - errno = EDOM; - } - /* - an infinite result here arises either from: - (A) (+/-0.)**negative (-> divide-by-zero) - (B) overflow of x**y with x and y finite - */ - else if (Py_IS_INFINITY(r)) { - if (x == 0.) - errno = EDOM; - else - errno = ERANGE; - } - } + if (PyFloat_CheckExact(arg)) { + x = PyFloat_AS_DOUBLE(arg); } - - if (errno && is_error(r)) - return NULL; else - return PyFloat_FromDouble(r); -} - - -static const double degToRad = Py_MATH_PI / 180.0; -static const double radToDeg = 180.0 / Py_MATH_PI; - -/*[clinic input] -math.degrees - - x: double - / - -Convert angle x from radians to degrees. -[clinic start generated code]*/ - -static PyObject * -math_degrees_impl(PyObject *module, double x) -/*[clinic end generated code: output=7fea78b294acd12f input=81e016555d6e3660]*/ -{ - return PyFloat_FromDouble(x * radToDeg); -} - - -/*[clinic input] -math.radians - - x: double - / - -Convert angle x from degrees to radians. -[clinic start generated code]*/ - -static PyObject * -math_radians_impl(PyObject *module, double x) -/*[clinic end generated code: output=34daa47caf9b1590 input=91626fc489fe3d63]*/ -{ - return PyFloat_FromDouble(x * degToRad); -} - - -/*[clinic input] -math.isfinite - - x: double - / - -Return True if x is neither an infinity nor a NaN, and False otherwise. -[clinic start generated code]*/ - -static PyObject * -math_isfinite_impl(PyObject *module, double x) -/*[clinic end generated code: output=8ba1f396440c9901 input=46967d254812e54a]*/ -{ - return PyBool_FromLong((long)Py_IS_FINITE(x)); -} - - -/*[clinic input] -math.isnan - - x: double - / - -Return True if x is a NaN (not a number), and False otherwise. -[clinic start generated code]*/ - -static PyObject * -math_isnan_impl(PyObject *module, double x) -/*[clinic end generated code: output=f537b4d6df878c3e input=935891e66083f46a]*/ -{ - return PyBool_FromLong((long)Py_IS_NAN(x)); -} - - -/*[clinic input] -math.isinf - - x: double - / - -Return True if x is a positive or negative infinity, and False otherwise. -[clinic start generated code]*/ - -static PyObject * -math_isinf_impl(PyObject *module, double x) -/*[clinic end generated code: output=9f00cbec4de7b06b input=32630e4212cf961f]*/ -{ - return PyBool_FromLong((long)Py_IS_INFINITY(x)); -} - - -/*[clinic input] -math.isclose -> bool - - a: double - b: double - * - rel_tol: double = 1e-09 - maximum difference for being considered "close", relative to the - magnitude of the input values - abs_tol: double = 0.0 - maximum difference for being considered "close", regardless of the - magnitude of the input values - -Determine whether two floating point numbers are close in value. - -Return True if a is close in value to b, and False otherwise. - -For the values to be considered close, the difference between them -must be smaller than at least one of the tolerances. - --inf, inf and NaN behave similarly to the IEEE 754 Standard. That -is, NaN is not close to anything, even itself. inf and -inf are -only close to themselves. -[clinic start generated code]*/ - -static int -math_isclose_impl(PyObject *module, double a, double b, double rel_tol, - double abs_tol) -/*[clinic end generated code: output=b73070207511952d input=f28671871ea5bfba]*/ -{ - double diff = 0.0; - - /* sanity check on the inputs */ - if (rel_tol < 0.0 || abs_tol < 0.0 ) { - PyErr_SetString(PyExc_ValueError, - "tolerances must be non-negative"); - return -1; - } - - if ( a == b ) { - /* short circuit exact equality -- needed to catch two infinities of - the same sign. And perhaps speeds things up a bit sometimes. - */ - return 1; - } - - /* This catches the case of two infinities of opposite sign, or - one infinity and one finite number. Two infinities of opposite - sign would otherwise have an infinite relative tolerance. - Two infinities of the same sign are caught by the equality check - above. - */ - - if (Py_IS_INFINITY(a) || Py_IS_INFINITY(b)) { - return 0; - } - - /* now do the regular computation - this is essentially the "weak" test from the Boost library - */ - - diff = fabs(b - a); - - return (((diff <= fabs(rel_tol * b)) || - (diff <= fabs(rel_tol * a))) || - (diff <= abs_tol)); -} - -static inline int -_check_long_mult_overflow(long a, long b) { - - /* From Python2's int_mul code: - - Integer overflow checking for * is painful: Python tried a couple ways, but - they didn't work on all platforms, or failed in endcases (a product of - -sys.maxint-1 has been a particular pain). - - Here's another way: - - The native long product x*y is either exactly right or *way* off, being - just the last n bits of the true product, where n is the number of bits - in a long (the delivered product is the true product plus i*2**n for - some integer i). - - The native double product (double)x * (double)y is subject to three - rounding errors: on a sizeof(long)==8 box, each cast to double can lose - info, and even on a sizeof(long)==4 box, the multiplication can lose info. - But, unlike the native long product, it's not in *range* trouble: even - if sizeof(long)==32 (256-bit longs), the product easily fits in the - dynamic range of a double. So the leading 50 (or so) bits of the double - product are correct. - - We check these two ways against each other, and declare victory if they're - approximately the same. Else, because the native long product is the only - one that can lose catastrophic amounts of information, it's the native long - product that must have overflowed. - - */ - - long longprod = (long)((unsigned long)a * b); - double doubleprod = (double)a * (double)b; - double doubled_longprod = (double)longprod; - - if (doubled_longprod == doubleprod) { - return 0; - } - - const double diff = doubled_longprod - doubleprod; - const double absdiff = diff >= 0.0 ? diff : -diff; - const double absprod = doubleprod >= 0.0 ? doubleprod : -doubleprod; - - if (32.0 * absdiff <= absprod) { - return 0; - } - - return 1; -} - -/*[clinic input] -math.prod - - iterable: object - / - * - start: object(c_default="NULL") = 1 - -Calculate the product of all the elements in the input iterable. - -The default start value for the product is 1. - -When the iterable is empty, return the start value. This function is -intended specifically for use with numeric values and may reject -non-numeric types. -[clinic start generated code]*/ - -static PyObject * -math_prod_impl(PyObject *module, PyObject *iterable, PyObject *start) -/*[clinic end generated code: output=36153bedac74a198 input=4c5ab0682782ed54]*/ -{ - PyObject *result = start; - PyObject *temp, *item, *iter; - - iter = PyObject_GetIter(iterable); - if (iter == NULL) { - return NULL; - } - - if (result == NULL) { - result = PyLong_FromLong(1); - if (result == NULL) { - Py_DECREF(iter); - return NULL; - } - } else { - Py_INCREF(result); - } -#ifndef SLOW_PROD - /* Fast paths for integers keeping temporary products in C. - * Assumes all inputs are the same type. - * If the assumption fails, default to use PyObjects instead. - */ - if (PyLong_CheckExact(result)) { - int overflow; - long i_result = PyLong_AsLongAndOverflow(result, &overflow); - /* If this already overflowed, don't even enter the loop. */ - if (overflow == 0) { - Py_DECREF(result); - result = NULL; - } - /* Loop over all the items in the iterable until we finish, we overflow - * or we found a non integer element */ - while(result == NULL) { - item = PyIter_Next(iter); - if (item == NULL) { - Py_DECREF(iter); - if (PyErr_Occurred()) { - return NULL; - } - return PyLong_FromLong(i_result); - } - if (PyLong_CheckExact(item)) { - long b = PyLong_AsLongAndOverflow(item, &overflow); - if (overflow == 0 && !_check_long_mult_overflow(i_result, b)) { - long x = i_result * b; - i_result = x; - Py_DECREF(item); - continue; - } - } - /* Either overflowed or is not an int. - * Restore real objects and process normally */ - result = PyLong_FromLong(i_result); - if (result == NULL) { - Py_DECREF(item); - Py_DECREF(iter); - return NULL; - } - temp = PyNumber_Multiply(result, item); - Py_DECREF(result); - Py_DECREF(item); - result = temp; - if (result == NULL) { - Py_DECREF(iter); - return NULL; - } - } - } - - /* Fast paths for floats keeping temporary products in C. - * Assumes all inputs are the same type. - * If the assumption fails, default to use PyObjects instead. - */ - if (PyFloat_CheckExact(result)) { - double f_result = PyFloat_AS_DOUBLE(result); - Py_DECREF(result); - result = NULL; - while(result == NULL) { - item = PyIter_Next(iter); - if (item == NULL) { - Py_DECREF(iter); - if (PyErr_Occurred()) { - return NULL; - } - return PyFloat_FromDouble(f_result); - } - if (PyFloat_CheckExact(item)) { - f_result *= PyFloat_AS_DOUBLE(item); - Py_DECREF(item); - continue; - } - if (PyLong_CheckExact(item)) { - long value; - int overflow; - value = PyLong_AsLongAndOverflow(item, &overflow); - if (!overflow) { - f_result *= (double)value; - Py_DECREF(item); - continue; - } - } - result = PyFloat_FromDouble(f_result); - if (result == NULL) { - Py_DECREF(item); - Py_DECREF(iter); - return NULL; - } - temp = PyNumber_Multiply(result, item); - Py_DECREF(result); - Py_DECREF(item); - result = temp; - if (result == NULL) { - Py_DECREF(iter); - return NULL; - } - } - } -#endif - /* Consume rest of the iterable (if any) that could not be handled - * by specialized functions above.*/ - for(;;) { - item = PyIter_Next(iter); - if (item == NULL) { - /* error, or end-of-sequence */ - if (PyErr_Occurred()) { - Py_DECREF(result); - result = NULL; - } - break; - } - temp = PyNumber_Multiply(result, item); - Py_DECREF(result); - Py_DECREF(item); - result = temp; - if (result == NULL) - break; - } - Py_DECREF(iter); - return result; -} - - -/*[clinic input] -math.perm - - n: object - k: object = None - / - -Number of ways to choose k items from n items without repetition and with order. - -Evaluates to n! / (n - k)! when k <= n and evaluates -to zero when k > n. - -If k is not specified or is None, then k defaults to n -and the function returns n!. - -Raises TypeError if either of the arguments are not integers. -Raises ValueError if either of the arguments are negative. -[clinic start generated code]*/ - -static PyObject * -math_perm_impl(PyObject *module, PyObject *n, PyObject *k) -/*[clinic end generated code: output=e021a25469653e23 input=5311c5a00f359b53]*/ -{ - PyObject *result = NULL, *factor = NULL; - int overflow, cmp; - long long i, factors; - - if (k == Py_None) { - return math_factorial(module, n); - } - n = PyNumber_Index(n); - if (n == NULL) { - return NULL; - } - if (!PyLong_CheckExact(n)) { - Py_SETREF(n, _PyLong_Copy((PyLongObject *)n)); - if (n == NULL) { - return NULL; - } - } - k = PyNumber_Index(k); - if (k == NULL) { - Py_DECREF(n); - return NULL; - } - if (!PyLong_CheckExact(k)) { - Py_SETREF(k, _PyLong_Copy((PyLongObject *)k)); - if (k == NULL) { - Py_DECREF(n); - return NULL; - } - } - - if (Py_SIZE(n) < 0) { - PyErr_SetString(PyExc_ValueError, - "n must be a non-negative integer"); - goto error; - } - if (Py_SIZE(k) < 0) { - PyErr_SetString(PyExc_ValueError, - "k must be a non-negative integer"); - goto error; - } - - cmp = PyObject_RichCompareBool(n, k, Py_LT); - if (cmp != 0) { - if (cmp > 0) { - result = PyLong_FromLong(0); - goto done; - } - goto error; - } - - factors = PyLong_AsLongLongAndOverflow(k, &overflow); - if (overflow > 0) { - PyErr_Format(PyExc_OverflowError, - "k must not exceed %lld", - LLONG_MAX); - goto error; - } - else if (factors == -1) { - /* k is nonnegative, so a return value of -1 can only indicate error */ - goto error; - } - - if (factors == 0) { - result = PyLong_FromLong(1); - goto done; - } - - result = n; - Py_INCREF(result); - if (factors == 1) { - goto done; - } - - factor = n; - Py_INCREF(factor); - for (i = 1; i < factors; ++i) { - Py_SETREF(factor, PyNumber_Subtract(factor, _PyLong_One)); - if (factor == NULL) { - goto error; - } - Py_SETREF(result, PyNumber_Multiply(result, factor)); - if (result == NULL) { - goto error; - } - } - Py_DECREF(factor); - -done: - Py_DECREF(n); - Py_DECREF(k); - return result; - -error: - Py_XDECREF(factor); - Py_XDECREF(result); - Py_DECREF(n); - Py_DECREF(k); - return NULL; -} - - -/*[clinic input] -math.comb - - n: object - k: object - / - -Number of ways to choose k items from n items without repetition and without order. - -Evaluates to n! / (k! * (n - k)!) when k <= n and evaluates -to zero when k > n. - -Also called the binomial coefficient because it is equivalent -to the coefficient of k-th term in polynomial expansion of the -expression (1 + x)**n. - -Raises TypeError if either of the arguments are not integers. -Raises ValueError if either of the arguments are negative. - -[clinic start generated code]*/ - -static PyObject * -math_comb_impl(PyObject *module, PyObject *n, PyObject *k) -/*[clinic end generated code: output=bd2cec8d854f3493 input=9a05315af2518709]*/ -{ - PyObject *result = NULL, *factor = NULL, *temp; - int overflow, cmp; - long long i, factors; - - n = PyNumber_Index(n); - if (n == NULL) { - return NULL; - } - if (!PyLong_CheckExact(n)) { - Py_SETREF(n, _PyLong_Copy((PyLongObject *)n)); - if (n == NULL) { - return NULL; - } - } - k = PyNumber_Index(k); - if (k == NULL) { - Py_DECREF(n); - return NULL; - } - if (!PyLong_CheckExact(k)) { - Py_SETREF(k, _PyLong_Copy((PyLongObject *)k)); - if (k == NULL) { - Py_DECREF(n); - return NULL; - } - } - - if (Py_SIZE(n) < 0) { - PyErr_SetString(PyExc_ValueError, - "n must be a non-negative integer"); - goto error; - } - if (Py_SIZE(k) < 0) { - PyErr_SetString(PyExc_ValueError, - "k must be a non-negative integer"); - goto error; - } - - /* k = min(k, n - k) */ - temp = PyNumber_Subtract(n, k); - if (temp == NULL) { - goto error; - } - if (Py_SIZE(temp) < 0) { - Py_DECREF(temp); - result = PyLong_FromLong(0); - goto done; - } - cmp = PyObject_RichCompareBool(temp, k, Py_LT); - if (cmp > 0) { - Py_SETREF(k, temp); - } - else { - Py_DECREF(temp); - if (cmp < 0) { - goto error; - } - } - - factors = PyLong_AsLongLongAndOverflow(k, &overflow); - if (overflow > 0) { - PyErr_Format(PyExc_OverflowError, - "min(n - k, k) must not exceed %lld", - LLONG_MAX); - goto error; - } - if (factors == -1) { - /* k is nonnegative, so a return value of -1 can only indicate error */ - goto error; - } - - if (factors == 0) { - result = PyLong_FromLong(1); - goto done; - } - - result = n; - Py_INCREF(result); - if (factors == 1) { - goto done; - } - - factor = n; - Py_INCREF(factor); - for (i = 1; i < factors; ++i) { - Py_SETREF(factor, PyNumber_Subtract(factor, _PyLong_One)); - if (factor == NULL) { - goto error; - } - Py_SETREF(result, PyNumber_Multiply(result, factor)); - if (result == NULL) { - goto error; - } - - temp = PyLong_FromUnsignedLongLong((unsigned long long)i + 1); - if (temp == NULL) { - goto error; - } - Py_SETREF(result, PyNumber_FloorDivide(result, temp)); - Py_DECREF(temp); - if (result == NULL) { - goto error; + { + x = PyFloat_AsDouble(arg); + if (x == -1.0 && PyErr_Occurred()) { + goto exit; } } - Py_DECREF(factor); - -done: - Py_DECREF(n); - Py_DECREF(k); - return result; - -error: - Py_XDECREF(factor); - Py_XDECREF(result); - Py_DECREF(n); - Py_DECREF(k); - return NULL; -} - - -/*[clinic input] -math.nextafter - - x: double - y: double - / - -Return the next floating-point value after x towards y. -[clinic start generated code]*/ - -static PyObject * -math_nextafter_impl(PyObject *module, double x, double y) -/*[clinic end generated code: output=750c8266c1c540ce input=02b2d50cd1d9f9b6]*/ -{ -#if defined(_AIX) - if (x == y) { - /* On AIX 7.1, libm nextafter(-0.0, +0.0) returns -0.0. - Bug fixed in bos.adt.libm 7.2.2.0 by APAR IV95512. */ - return PyFloat_FromDouble(y); - } -#endif - return PyFloat_FromDouble(nextafter(x, y)); -} - - -/*[clinic input] -math.ulp -> double - - x: double - / - -Return the value of the least significant bit of the float x. -[clinic start generated code]*/ - -static double -math_ulp_impl(PyObject *module, double x) -/*[clinic end generated code: output=f5207867a9384dd4 input=31f9bfbbe373fcaa]*/ -{ - if (Py_IS_NAN(x)) { - return x; - } - x = fabs(x); - if (Py_IS_INFINITY(x)) { - return x; + _return_value = math_ulp_impl(module, x); + if ((_return_value == -1.0) && PyErr_Occurred()) { + goto exit; } - double inf = m_inf(); - double x2 = nextafter(x, inf); - if (Py_IS_INFINITY(x2)) { - /* special case: x is the largest positive representable float */ - x2 = nextafter(x, -inf); - return x - x2; - } - return x2 - x; -} + return_value = PyFloat_FromDouble(_return_value); - -static PyMethodDef math_methods[] = { - {"acos", math_acos, METH_O, math_acos_doc}, - {"acosh", math_acosh, METH_O, math_acosh_doc}, - {"asin", math_asin, METH_O, math_asin_doc}, - {"asinh", math_asinh, METH_O, math_asinh_doc}, - {"atan", math_atan, METH_O, math_atan_doc}, - {"atan2", (PyCFunction)(void(*)(void))math_atan2, METH_FASTCALL, math_atan2_doc}, - {"atanh", math_atanh, METH_O, math_atanh_doc}, - MATH_CEIL_METHODDEF - {"copysign", (PyCFunction)(void(*)(void))math_copysign, METH_FASTCALL, math_copysign_doc}, - {"cos", math_cos, METH_O, math_cos_doc}, - {"cosh", math_cosh, METH_O, math_cosh_doc}, - MATH_DEGREES_METHODDEF - MATH_DIST_METHODDEF - {"erf", math_erf, METH_O, math_erf_doc}, - {"erfc", math_erfc, METH_O, math_erfc_doc}, - {"exp", math_exp, METH_O, math_exp_doc}, - {"expm1", math_expm1, METH_O, math_expm1_doc}, - {"fabs", math_fabs, METH_O, math_fabs_doc}, - MATH_FACTORIAL_METHODDEF - MATH_FLOOR_METHODDEF - MATH_FMOD_METHODDEF - MATH_FREXP_METHODDEF - MATH_FSUM_METHODDEF - {"gamma", math_gamma, METH_O, math_gamma_doc}, - MATH_GCD_METHODDEF - {"hypot", (PyCFunction)(void(*)(void))math_hypot, METH_FASTCALL, math_hypot_doc}, - MATH_ISCLOSE_METHODDEF - MATH_ISFINITE_METHODDEF - MATH_ISINF_METHODDEF - MATH_ISNAN_METHODDEF - MATH_ISQRT_METHODDEF - MATH_LDEXP_METHODDEF - {"lgamma", math_lgamma, METH_O, math_lgamma_doc}, - MATH_LOG_METHODDEF - {"log1p", math_log1p, METH_O, math_log1p_doc}, - MATH_LOG10_METHODDEF - MATH_LOG2_METHODDEF - MATH_MODF_METHODDEF - MATH_POW_METHODDEF - MATH_RADIANS_METHODDEF - {"remainder", (PyCFunction)(void(*)(void))math_remainder, METH_FASTCALL, math_remainder_doc}, - {"sin", math_sin, METH_O, math_sin_doc}, - {"sinh", math_sinh, METH_O, math_sinh_doc}, - {"sqrt", math_sqrt, METH_O, math_sqrt_doc}, - {"tan", math_tan, METH_O, math_tan_doc}, - {"tanh", math_tanh, METH_O, math_tanh_doc}, - MATH_TRUNC_METHODDEF - MATH_PROD_METHODDEF - MATH_PERM_METHODDEF - MATH_COMB_METHODDEF - MATH_NEXTAFTER_METHODDEF - MATH_ULP_METHODDEF - {NULL, NULL} /* sentinel */ -}; - - -PyDoc_STRVAR(module_doc, -"This module provides access to the mathematical functions\n" -"defined by the C standard."); - - -static struct PyModuleDef mathmodule = { - PyModuleDef_HEAD_INIT, - "math", - module_doc, - -1, - math_methods, - NULL, - NULL, - NULL, - NULL -}; - -PyMODINIT_FUNC -PyInit_math(void) -{ - PyObject *m; - - m = PyModule_Create(&mathmodule); - if (m == NULL) - goto finally; - - PyModule_AddObject(m, "pi", PyFloat_FromDouble(Py_MATH_PI)); - PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E)); - PyModule_AddObject(m, "tau", PyFloat_FromDouble(Py_MATH_TAU)); /* 2pi */ - PyModule_AddObject(m, "inf", PyFloat_FromDouble(m_inf())); -#if !defined(PY_NO_SHORT_FLOAT_REPR) || defined(Py_NAN) - PyModule_AddObject(m, "nan", PyFloat_FromDouble(m_nan())); -#endif - - finally: - return m; +exit: + return return_value; } +/*[clinic end generated code: output=f8daa185c043a7b7 input=a9049054013a1b77]*/