Index: Modules/mathmodule.c =================================================================== --- Modules/mathmodule.c (revision 81150) +++ Modules/mathmodule.c (working copy) @@ -1129,56 +1129,173 @@ Return an accurate floating point sum of values in the iterable.\n\ Assumes IEEE-754 floating point arithmetic."); +/* Divide-and-conquer factorial algorithm + * + * Based on the formula and psuedo-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 factoriazation algorithm. + */ + +/* i-th odd number */ +#define ODD(i) (i << 1 | 1) +/* Compute product(ODD(i) for i in range(first, last+1)) */ static PyObject * +factorial_partial_product(unsigned long first, unsigned long last) +{ + PyObject *result; + PyObject *a, *x, *y; + Py_ssize_t n, i; + + n = last - first + 1; + if (n <= 0) + return PyLong_FromLong(1L); + if (n == 1) + return PyLong_FromUnsignedLong(ODD(first)); + + /* a = [ODD(i) for i in range(first, last + 1)] */ + a = PyList_New(last - first + 1); + if (a == NULL) + return NULL; + for (i = 0; i < n; ++i) { + result = PyLong_FromUnsignedLong(ODD(first + i)); + if (result == NULL) + goto done; + PyList_SET_ITEM(a, i, result); + } + while (n > 1) { + for (i = (n >> 1) - 1; i >= 0; --i) { + x = PyList_GET_ITEM(a, i); + y = PyList_GET_ITEM(a, n - i - 1); + result = PyNumber_Multiply(x, y); + if (result == NULL) + goto done; + PyList_SET_ITEM(a, i, result); + PyList_SET_ITEM(a, n - i - 1, NULL); + Py_DECREF(x); + Py_DECREF(y); + } + n = (n >> 1) + (n & 1); + } + /* Release reference held by the list */ + PyList_SET_ITEM(a, 0, NULL); + done: + Py_DECREF(a); + return result; +} + +static unsigned long +bit_length(unsigned long n) +{ + unsigned long len = 0; + while (n != 0) { + ++len; + n >>= 1; + } + return len; +} + +static unsigned long +bit_count(unsigned long n) +{ + unsigned long count = 0; + while (n != 0) { + ++count; + n &= n - 1; /* clear least significant bit */ + } + return count; +} + +static PyObject * +factorial_loop(unsigned long n) +{ + long i, v; + PyObject *p, *r; + PyObject *part, *tmp; + + p = PyLong_FromLong(1L); + if (p == NULL) + return NULL; + Py_INCREF(p); + r = p; + + for (i = bit_length(n) - 2; i >= 0; --i) { + v = n >> i; + if (v > 2) { + part = factorial_partial_product(((v >> 1) + 1) >> 1, (v - 1) >> 1); + if (part == NULL) + goto error; + + tmp = PyNumber_Multiply(p, part); + Py_DECREF(part); + if (tmp == NULL) + goto error; + Py_DECREF(p); + p = tmp; + + tmp = PyNumber_Multiply(r, p); + if (tmp == NULL) + goto error; + Py_DECREF(r); + r = tmp; + } + } + Py_DECREF(p); + return r; + error: + Py_DECREF(p); + Py_DECREF(r); + return NULL; +} + +static PyObject * math_factorial(PyObject *self, PyObject *arg) { - long i, x; - PyObject *result, *iobj, *newresult; + PyObject *result; + PyObject *r; /* result without trailing zeros */ + PyObject *ntz; /* number of trailing zeros */ + long n; + if (PyFloat_Check(arg)) { - 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) + n = (long) dx; + } else { + n = PyLong_AsLong(arg); + if (n == -1 && PyErr_Occurred()) return NULL; - x = PyLong_AsLong(lx); - Py_DECREF(lx); } - else - x = PyLong_AsLong(arg); - if (x == -1 && PyErr_Occurred()) - return NULL; - if (x < 0) { + if (n < 0) { PyErr_SetString(PyExc_ValueError, - "factorial() not defined for negative values"); + "factorial() not defined for negative values"); return NULL; } - result = (PyObject *)PyLong_FromLong(1); - if (result == NULL) + if (n <= 12) { + static const long lookup[] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, + 362880, 3628800, 39916800, 479001600}; + return PyLong_FromLong(lookup[n]); + } + + r = factorial_loop(n); + if (r == NULL) return NULL; - for (i=1 ; i<=x ; i++) { - iobj = (PyObject *)PyLong_FromLong(i); - if (iobj == NULL) - goto error; - newresult = PyNumber_Multiply(result, iobj); - Py_DECREF(iobj); - if (newresult == NULL) - goto error; - Py_DECREF(result); - result = newresult; + + ntz = PyLong_FromLong(n - bit_count(n)); + if (ntz != NULL) { + result = PyNumber_Lshift(r, ntz); + Py_DECREF(ntz); } + + Py_DECREF(r); return result; - -error: - Py_DECREF(result); - return NULL; } PyDoc_STRVAR(math_factorial_doc,