diff --git Modules/mathmodule.c Modules/mathmodule.c
index 76d7906..08cb981 100644
--- Modules/mathmodule.c
+++ Modules/mathmodule.c
@@ -1129,56 +1129,148 @@ PyDoc_STRVAR(math_fsum_doc,
 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.
+ */
+
+/* Maximum value such that multiply_cutoff * multiply_cutoff fits in a long */
+static long multiply_cutoff;
+
+/* Compute product(k for k in range(n, m+1) if k is odd) using divide
+ * and conquer. */
+static PyObject *
+factorial_part_product(long n, long m)
+{
+    long k;
+    PyObject *left = NULL, *right = NULL, *result = NULL;
+
+    if (m <= (n + 1))
+        return PyLong_FromLong(n);
+    if (m == (n + 2)) {
+        if (m < multiply_cutoff && n < multiply_cutoff)
+            return PyLong_FromLong(m*n);
+        left = PyLong_FromLong(n);
+        if (left == NULL) goto done;
+        right = PyLong_FromLong(m);
+        if (right == NULL) goto done;
+    } else {
+        k = (n + m) / 2;
+        if ((k & 1) != 1)
+            k = k - 1;
+        left = factorial_part_product(n, k);
+        if (left == NULL) goto done;
+        right = factorial_part_product(k + 2, m);
+        if (right == NULL) goto done;
+    }
+    result = PyNumber_Multiply(left, right);
+done:
+    Py_XDECREF(left);
+    Py_XDECREF(right);
+    return result;
+}
+
+static PyObject *
+factorial_loop(long n)
+{
+    long i, v;
+    PyObject *part, *tmp, *p = NULL, *r = NULL, *result = NULL;
+
+    p = PyLong_FromLong(1);
+    if (p == NULL)
+        goto done;
+    r = PyLong_FromLong(1);
+    if (r == NULL)
+        goto done;
+
+    for (i = log2(n)-1; i >= 0; i--) {
+        v = n >> i;
+        if (v <= 2) continue;
+
+        part = factorial_part_product(v / 2 + 1 + ((v / 2) & 1),
+                                      v - 1 + (v & 1));
+        if (part == NULL)
+            goto done;
+        tmp = PyNumber_Multiply(p, part);
+        if (tmp == NULL)
+            goto done;
+        Py_DECREF(p);
+        p = tmp;
+
+        tmp = PyNumber_Multiply(r, p);
+        if (tmp == NULL)
+            goto done;
+        Py_DECREF(r);
+        r = tmp;
+    }
+
+    result = r;
+    r = NULL;
+
+done:
+    Py_XDECREF(p);
+    Py_XDECREF(r);
+    return result;
+}
+
 static PyObject *
 math_factorial(PyObject *self, PyObject *arg)
 {
-    long i, x;
-    PyObject *result, *iobj, *newresult;
+    long n, tmp;
+    PyObject *result = NULL, *r = NULL;
+    PyObject *nminusnumbits_ob = NULL;
+    size_t num_bits;
 
     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");
-        return NULL;
+                        "factorial() not defined for negative values");
+        goto done;
     }
 
-    result = (PyObject *)PyLong_FromLong(1);
-    if (result == 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;
+    if (n <= 12) {
+        static const long lookup[] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320,
+                                      362880, 3628800, 39916800, 479001600};
+        result = PyLong_FromLong(lookup[n]);
+        goto done;
     }
+
+    /* Count the number of set bits in n */
+    num_bits = 0;
+    for (tmp = n; tmp; tmp >>= 1)
+        num_bits += tmp & 1;
+
+    r = factorial_loop(n);
+    if (r == NULL)
+        goto done;
+
+    tmp = n - num_bits;
+    nminusnumbits_ob = PyLong_FromLong(tmp);
+    if (nminusnumbits_ob == NULL)
+        goto done;
+    result = PyNumber_Lshift(r, nminusnumbits_ob);
+
+done:
+    Py_XDECREF(nminusnumbits_ob);
+    Py_XDECREF(r);
     return result;
-
-error:
-    Py_DECREF(result);
-    return NULL;
 }
 
 PyDoc_STRVAR(math_factorial_doc,
@@ -1695,6 +1787,8 @@ PyInit_math(void)
     if (m == NULL)
         goto finally;
 
+    multiply_cutoff = (long) sqrt(LONG_MAX);
+
     PyModule_AddObject(m, "pi", PyFloat_FromDouble(Py_MATH_PI));
     PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E));