Index: Misc/NEWS =================================================================== --- Misc/NEWS (revision 81174) +++ Misc/NEWS (working copy) @@ -1120,6 +1120,12 @@ Extension Modules ----------------- +- Issue #8692: Optimize math.factorial: replace the previous naive + algorithm with an improved 'binary-split' algorithm that uses fewer + multiplications and allows many of the multiplications to be + performed using plain C integer arithmetic instead of PyLong + arithmetic. Also uses a lookup table for small arguments. + - Issue #8674: Fixed a number of incorrect or undefined-behaviour-inducing overflow checks in the audioop module. Index: Lib/test/test_math.py =================================================================== --- Lib/test/test_math.py (revision 81174) +++ Lib/test/test_math.py (working copy) @@ -365,18 +365,18 @@ self.ftest('fabs(1)', math.fabs(1), 1) def testFactorial(self): - def fact(n): - result = 1 - for i in range(1, int(n)+1): - result *= i - return result - values = list(range(10)) + [50, 100, 500] - random.shuffle(values) - for x in values: - for cast in (int, float): - self.assertEqual(math.factorial(cast(x)), fact(x), (x, fact(x), math.factorial(x))) + self.assertEqual(math.factorial(0), 1) + self.assertEqual(math.factorial(0.0), 1) + total = 1 + for i in range(1, 1000): + total *= i + self.assertEqual(math.factorial(i), total) + self.assertEqual(math.factorial(float(i)), total) self.assertRaises(ValueError, math.factorial, -1) + self.assertRaises(ValueError, math.factorial, -1.0) self.assertRaises(ValueError, math.factorial, math.pi) + self.assertRaises(OverflowError, math.factorial, sys.maxsize+1) + self.assertRaises(OverflowError, math.factorial, 10e100) def testFloor(self): self.assertRaises(TypeError, math.floor) Index: Modules/mathmodule.c =================================================================== --- Modules/mathmodule.c (revision 81174) +++ Modules/mathmodule.c (working copy) @@ -1129,11 +1129,211 @@ Return an accurate floating point sum of values in the iterable.\n\ Assumes IEEE-754 floating point arithmetic."); +/* Find the index of the highest set bit. Equivalent to floor(lg(x))+1. + * Also equivalent to: bitwidth_of_type - count_leading_zero_bits(x) + */ + +/* XXX: This routine does more or less the same thing as + * bits_in_digit() in Objects/longobject.c.Someday it would be nice to + * consolidate them. On BSD, there's a library function called fls() + * that we could use, and GCC provides __builtin_clz(). + */ + +static unsigned long +find_last_set_bit(unsigned long n) +{ + unsigned long len = 0; + while (n != 0) { + ++len; + n >>= 1; + } + return len; +} + +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; +} + +/* 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. + * + * Notes on the algorithm + * ---------------------- + * + * The function factorial_loop computes the odd part (i.e., the greatest odd + * divisor) of factorial(n), using the formula: + * + * odd-part-of-factorial(n) = + * + * product_{i >= 0} product_{0 < j <= n / 2**i, j odd} j + * + * Example: odd-part-of-factorial(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(n, m) for given n and m. 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 even part of the factorial is computed independently in the main + * math_factorial function. Its value is 2**k where k = n//2 + n//4 + n//8 + + * .... It can be shown (e.g., by complete induction on n) that k 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(n, m, 2)) using divide and + * conquer. Assumes n and m are odd and m > n. max_bits must be >= + * find_last_set_bit(m-2). */ + static PyObject * +factorial_partial_product(unsigned long n, unsigned long m, + unsigned long max_bits) +{ + unsigned long k, 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 find_last_set_bit(z). The + * product x*y will require at most + * find_last_set_bit(x)+find_last_set_bit(y) bits to store, based + * on the idea that lg product = lg x + lg y. + * + * We know that m is the largest number to be multiplied. From + * there, we have: + * find_last_set_bit(answer) <= num_operands * find_last_set_bit(m) + */ + + num_operands = (m - n) / 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 total = n; + for (n += 2; n < m; n += 2) + total *= n; + return PyLong_FromUnsignedLong(total); + } + + /* k = midpoint of range(n, m), rounded up to next odd number. */ + k = (n + num_operands) | 1; + left = factorial_partial_product(n, k, find_last_set_bit(k - 2)); + if (left == NULL) + goto error; + right = factorial_partial_product(k, m, max_bits); + if (right == NULL) + goto error; + result = PyNumber_Multiply(left, right); + + error: + Py_XDECREF(left); + Py_XDECREF(right); + return result; +} + +/* factorial_loop: compute the odd part of factorial(n). */ + +static PyObject * +factorial_loop(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 = find_last_set_bit(n) - 2; i >= 0; i--) { + v = n >> i; + if (v <= 2) + continue; + lower = upper; + 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, + find_last_set_bit(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. */ + + /* r *= inner; */ + tmp = PyNumber_Multiply(outer, inner); + if (tmp == NULL) + goto error; + Py_DECREF(outer); + outer = tmp; + } + + goto done; + + error: + Py_DECREF(outer); + done: + Py_DECREF(inner); + return outer; +} + +static PyObject * math_factorial(PyObject *self, PyObject *arg) { - long i, x; - PyObject *result, *iobj, *newresult; + long x; + PyObject *result = NULL, *r = NULL, *nminusnumbits = NULL; if (PyFloat_Check(arg)) { PyObject *lx; @@ -1160,25 +1360,28 @@ return NULL; } - result = (PyObject *)PyLong_FromLong(1); - if (result == NULL) + if (x <= 12) { + static const unsigned lookup[] = { + 1, 1, 2, 6, 24, 120, 720, 5040, 40320, + 362880, 3628800, 39916800, 479001600 + }; + result = PyLong_FromLong(lookup[x]); + return result; + } + + r = factorial_loop(x); + 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; - } + + nminusnumbits = PyLong_FromLong(x - count_set_bits(x)); + if (nminusnumbits == NULL) + goto error; + result = PyNumber_Lshift(r, nminusnumbits); + Py_DECREF(nminusnumbits); + + error: + Py_DECREF(r); return result; - -error: - Py_DECREF(result); - return NULL; } PyDoc_STRVAR(math_factorial_doc,