Index: Include/longobject.h
===================================================================
--- Include/longobject.h	(revision 60941)
+++ Include/longobject.h	(working copy)
@@ -119,6 +119,9 @@
    a leading "0", instead of the prefix "0o" */
 PyAPI_FUNC(PyObject *) _PyLong_Format(PyObject *aa, int base, int addL, int newstyle);
 
+/* For use by the gcd function in mathmodule.c */
+PyAPI_FUNC(PyLongObject *) _PyLong_GCD(PyLongObject *, PyLongObject *);
+
 #ifdef __cplusplus
 }
 #endif
Index: Objects/longobject.c
===================================================================
--- Objects/longobject.c	(revision 60941)
+++ Objects/longobject.c	(working copy)
@@ -3265,6 +3265,188 @@
 	return PyInt_FromLong(x);
 }
 
+static int
+bits_in_digit(digit d)
+{
+	int nbits;
+
+	if (d==0)
+		return 0;
+	nbits = 1;
+	if (d & ~255) {
+		d >>= 8;
+		nbits += 8;
+	}
+	if (d & ~15) {
+		d >>= 4;
+		nbits += 4;
+	}
+	if (d & ~3) {
+		d >>= 2;
+		nbits += 2;
+	}
+	if (d & ~1) {
+		d >>= 1;
+		nbits += 1;
+	}
+	return nbits;
+}
+
+PyLongObject *
+_PyLong_GCD(PyLongObject *a, PyLongObject *b)
+{
+	PyLongObject *c, *d;
+        stwodigits x, y, q, s, t, c_carry, d_carry;
+        digit A, B, C, D;
+	int nbits, k;
+	Py_ssize_t size_a, size_b;
+        digit *a_digit, *b_digit, *c_digit, *d_digit, *a_end, *b_end;
+
+	/* Initial reduction: make sure that 0 <= b <= a. */
+	a = (PyLongObject *)long_abs(a);
+	b = (PyLongObject *)long_abs(b);
+	if (long_compare(a, b) < 0) {
+		d = a;
+		a = b;
+		b = d;
+	}
+	/* We now own references to a and b */
+
+	/* reduce until a fits into 2 digits */
+	while ((size_a = Py_SIZE(a)) > 2) {
+		nbits = bits_in_digit(a->ob_digit[size_a-1]);
+		/* extract top 2*PyLong_SHIFT bits of a into x, along with
+		   corresponding bits of b into y */
+		size_b = Py_SIZE(b);
+		x = ((a->ob_digit[size_a-1] << (2*PyLong_SHIFT-nbits)) |
+		     (a->ob_digit[size_a-2] << (PyLong_SHIFT-nbits)) |
+		     (a->ob_digit[size_a-3] >> nbits));
+
+		y = ((size_b >= size_a - 2 ? b->ob_digit[size_a-3] >> nbits : 0) |
+		     (size_b >= size_a - 1 ? b->ob_digit[size_a-2] << (PyLong_SHIFT-nbits) : 0) |
+		     (size_b >= size_a ? b->ob_digit[size_a-1] << (2*PyLong_SHIFT-nbits) : 0));
+
+		/* inner loop of Lehmer's algorithm; A, B, C, D never grow
+		larger than PyLong_MASK during the algorithm. */
+		A = 1; B = 0; C = 0; D = 1;
+		for (k=0;; k++) {
+			if (y-C == 0)
+				break;
+			q = (x+(A-1))/(y-C);
+			s = B+q*D;
+			t = x-q*y;
+			if (s > t)
+				break;
+			x = y; y = t;
+			t = A+q*C; A = D; B = C; C = (digit)s; D = (digit)t;
+		}
+
+		if (k == 0) {
+			/* no progress; do a Euclidean step */
+			if (Py_SIZE(b) == 0) {
+				Py_DECREF(b);
+				return a;
+			}
+			if (l_divmod(a, b, NULL, &d) < 0) {
+				Py_DECREF(a);
+				Py_DECREF(b);
+				return NULL;
+			}
+			Py_DECREF(a);
+			a = b;
+			b = d;
+			continue;
+		}
+
+		/*
+		  a, b = A*b-B*a, D*a-C*b if k is odd
+		  a, b = A*a-B*b, D*b-C*a if k is even
+		*/
+		c = _PyLong_New(size_a);
+		if (c == NULL) {
+			Py_DECREF(a);
+			Py_DECREF(b);
+			return NULL;
+		}
+
+		d = _PyLong_New(size_a);
+		if (d == NULL) {
+			Py_DECREF(a);
+			Py_DECREF(b);
+			Py_DECREF(c);
+			return NULL;
+		}
+		a_end = a->ob_digit + size_a;
+		b_end = b->ob_digit + size_b;
+
+		/* compute new a and new b in parallel */
+		a_digit = a->ob_digit;
+		b_digit = b->ob_digit;
+		c_digit = c->ob_digit;
+		d_digit = d->ob_digit;
+		c_carry = 0;
+		d_carry = 0;
+		if (k&1) {
+			while (b_digit < b_end) {
+				c_carry += (A * *b_digit) - (B * *a_digit);
+				d_carry += (D * *a_digit++) - (C * *b_digit++);
+				*c_digit++ = (digit)(c_carry & PyLong_MASK);
+				*d_digit++ = (digit)(d_carry & PyLong_MASK);
+				c_carry >>= PyLong_SHIFT;
+				d_carry >>= PyLong_SHIFT;
+			}
+			while (a_digit < a_end) {
+				c_carry -= B * *a_digit;
+				d_carry += D * *a_digit++;
+				*c_digit++ = (digit)(c_carry & PyLong_MASK);
+				*d_digit++ = (digit)(d_carry & PyLong_MASK);
+				c_carry >>= PyLong_SHIFT;
+				d_carry >>= PyLong_SHIFT;
+			}
+		}
+		else {
+			while (b_digit < b_end) {
+				c_carry += (A * *a_digit) - (B * *b_digit);
+				d_carry += (D * *b_digit++) - (C * *a_digit++);
+				*c_digit++ = (digit)(c_carry & PyLong_MASK);
+				*d_digit++ = (digit)(d_carry & PyLong_MASK);
+				c_carry >>= PyLong_SHIFT;
+				d_carry >>= PyLong_SHIFT;
+			}
+			while (a_digit < a_end) {
+				c_carry += A * *a_digit;
+				d_carry -= C * *a_digit++;
+				*c_digit++ = (digit)(c_carry & PyLong_MASK);
+				*d_digit++ = (digit)(d_carry & PyLong_MASK);
+				c_carry >>= PyLong_SHIFT;
+				d_carry >>= PyLong_SHIFT;
+			}
+		}
+		assert(c_carry == 0);
+		assert(d_carry == 0);
+
+		Py_DECREF(a);
+		Py_DECREF(b);
+		a = long_normalize(c);
+		b = long_normalize(d);
+	}
+
+	/* a fits into a long, so b must too */
+	x = PyLong_AsLong((PyObject *)a);
+	y = PyLong_AsLong((PyObject *)b);
+	Py_DECREF(a);
+	Py_DECREF(b);
+
+	/* usual Euclidean algorithm for longs */
+	while (y != 0) {
+		t = y;
+		y = x % y;
+		x = t;
+	}
+	return (PyLongObject *)PyLong_FromLong(x);
+
+}
+
 static PyObject *
 long_float(PyObject *v)
 {
Index: Modules/mathmodule.c
===================================================================
--- Modules/mathmodule.c	(revision 60941)
+++ Modules/mathmodule.c	(working copy)
@@ -357,6 +357,82 @@
 Checks if float x is infinite (positive or negative)");
 
 
+static PyObject *
+math_gcd(PyObject *self, PyObject *args)
+{
+	PyObject *a, *b, *c;
+	long x, y, t;
+
+	if (!PyArg_UnpackTuple(args, "gcd", 2, 2, &a, &b))
+		return NULL;
+	if (PyInt_Check(a)) {
+		x = PyInt_AS_LONG(a);
+		if (PyInt_Check(b)) {
+			/* both a and b are integers;
+			   use the usual gcd algorithm */
+			y = PyInt_AS_LONG(b);
+			while (y != 0) {
+				t = x%y;
+				x = y;
+				y = t;
+			}
+			if (x < 0) {
+				/* need to take abs value; check whether -x overflows */
+				if ((unsigned long)(x) == 0-(unsigned long)(x)) {
+					PyObject *o = PyLong_FromLong(x);
+					if (o == NULL)
+						return NULL;
+					PyObject *result = PyNumber_Negative(o);
+					Py_DECREF(o);
+					return result;
+				}
+				else {
+					return PyInt_FromLong(-x);
+				}
+			}
+			else {
+				return PyInt_FromLong(x);
+			}
+		}
+		else if (PyLong_Check(b)) {
+			a = PyLong_FromLong(x);
+			c = (PyObject *)_PyLong_GCD((PyLongObject *)a, (PyLongObject *)b);
+			Py_DECREF(a);
+			return c;
+		}
+		else {
+			PyErr_SetString(PyExc_TypeError, "expected an integer");
+			return NULL;
+		}
+	}
+	else if (PyLong_Check(a)) {
+		if (PyInt_Check(b)) {
+			b = PyLong_FromLong(PyInt_AS_LONG(b));
+			c = (PyObject *)_PyLong_GCD((PyLongObject *)a, (PyLongObject *)b);
+			Py_DECREF(b);
+			return c;
+		}
+		else if (PyLong_Check(b)) {
+			c = (PyObject *)_PyLong_GCD((PyLongObject *)a, (PyLongObject *)b);
+			return c;
+		}
+		else {
+			PyErr_SetString(PyExc_TypeError, "expected an integer");
+			return NULL;
+		}
+	}
+	else {
+		PyErr_SetString(PyExc_TypeError, "expected an integer");
+		return NULL;
+	}		
+}
+
+PyDoc_STRVAR(math_gcd_doc,
+"gcd(x, y) -> int\n\
+greatest common divisor of x and y");
+
+
+
 static PyMethodDef math_methods[] = {
 	{"acos",	math_acos,	METH_O,		math_acos_doc},
 	{"asin",	math_asin,	METH_O,		math_asin_doc},
@@ -374,6 +450,7 @@
 	{"floor",	math_floor,	METH_O,		math_floor_doc},
 	{"fmod",	math_fmod,	METH_VARARGS,	math_fmod_doc},
 	{"frexp",	math_frexp,	METH_O,		math_frexp_doc},
+	{"gcd",         math_gcd,       METH_VARARGS,   math_gcd_doc},
 	{"hypot",	math_hypot,	METH_VARARGS,	math_hypot_doc},
 	{"isinf",	math_isinf,	METH_O,		math_isinf_doc},
 	{"isnan",	math_isnan,	METH_O,		math_isnan_doc},