diff -r b6d771d318fc src/org/python/modules/math.java
--- a/src/org/python/modules/math.java	Tue Apr 24 15:35:45 2012 -0700
+++ b/src/org/python/modules/math.java	Tue Apr 24 15:37:33 2012 -0700
@@ -12,6 +12,7 @@
 import org.python.core.PySystemState;
 import org.python.core.PyTuple;
 import org.python.modules.math_erf;
+import org.python.modules.math_gamma;
 
 public class math implements ClassDictInit {
     public static PyFloat pi = new PyFloat(Math.PI);
@@ -565,6 +566,14 @@
     }
 
 
+    public static double gamma(double x) {
+        return math_gamma.gamma(x);
+    }
+
+    public static double lgamma(double x) {
+        return math_gamma.lgamma(x);
+    }
+
     private static double calculateLongLog(PyLong v) {
         int exp[] = new int[1];
         double x = v.scaledDoubleValue(exp);
diff -r b6d771d318fc src/org/python/modules/math_gamma.java
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/org/python/modules/math_gamma.java	Tue Apr 24 15:37:33 2012 -0700
@@ -0,0 +1,454 @@
+/*
+	Floating-point logarithm of the gamma and lgamma function.
+
+        This is a simple translation of the code in Google's go math/lgamma.go
+*/
+
+// The original C code and the long comment below are
+// from FreeBSD's /usr/src/lib/msun/src/e_lgamma_r.c and
+// came with this notice.  The go code is a simplified
+// version of the original C.
+//
+// ====================================================
+// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+//
+// Developed at SunPro, a Sun Microsystems, Inc. business.
+// Permission to use, copy, modify, and distribute this
+// software is freely granted, provided that this notice
+// is preserved.
+// ====================================================
+
+package org.python.modules;
+
+import org.python.core.Py;
+
+public class math_gamma {
+    static final double[] _gamP = {
+            1.60119522476751861407e-04,
+            1.19135147006586384913e-03,
+            1.04213797561761569935e-02,
+            4.76367800457137231464e-02,
+            2.07448227648435975150e-01,
+            4.94214826801497100753e-01,
+            9.99999999999999996796e-01,
+    };
+    static final double[] _gamQ = {
+            -2.31581873324120129819e-05,
+            5.39605580493303397842e-04,
+            -4.45641913851797240494e-03,
+            1.18139785222060435552e-02,
+            3.58236398605498653373e-02,
+            -2.34591795718243348568e-01,
+            7.14304917030273074085e-02,
+            1.00000000000000000320e+00,
+    };
+    static final double[] _gamS = {
+            7.87311395793093628397e-04,
+            -2.29549961613378126380e-04,
+            -2.68132617805781232825e-03,
+            3.47222221605458667310e-03,
+            8.33333333333482257126e-02,
+    };
+
+    private static double stirling(double x) {
+        final double sqrtTwoPi = 2.506628274631000502417;
+        final double maxStirling = 143.01608;
+	double w = 1 / x;
+	w = 1 + w*((((_gamS[0]*w+_gamS[1])*w+_gamS[2])*w+_gamS[3])*w+_gamS[4]);
+	double y = Math.exp(x);
+        if (x > maxStirling) { // avoid Pow() overflow
+            double v = Math.pow(x, 0.5*x-0.25);
+            y = v * (v / y);
+        } else {
+            y = Math.pow(x, x-0.5) / y;
+        }
+	y = sqrtTwoPi * y * w;
+	return y;
+    }
+
+    private static double small(double x, double z) {
+	final double euler = 0.57721566490153286060651209008240243104215933593992; // A001620
+	if (x == 0) {
+            return Double.POSITIVE_INFINITY;
+	}
+	return z / ((1 + euler*x) * x);
+    }
+
+    public static double gamma(double x) {
+        if (Double.POSITIVE_INFINITY == x) {
+            return x;
+        }
+        if ((Double.NEGATIVE_INFINITY == x) ||
+            (0.0 == x) ||
+            (-0.0 == x) ||
+            (-1.0 == x) || // FIXME: ???
+            (-2.0 == x)) { // FIXME: ???
+            throw Py.ValueError(Double.toString(x));
+        }
+        double result = _gamma(x);
+        if (Double.isInfinite(result)) {
+            throw Py.OverflowError(Double.toString(x));
+        }
+
+        return result;
+    }
+
+    public static double _gamma(double x) {
+        // special cases
+        if ((Double.NEGATIVE_INFINITY == x) || Double.isNaN(x)) {
+            return x;
+        }
+
+        if ((x < -170.5674972726612) || (x > 171.61447887182298)) {
+            return Double.POSITIVE_INFINITY;
+        }
+
+        double q = Math.abs(x);
+        double p = Math.floor(q);
+        if (q > 33) {
+            if (x >= 0) {
+                return stirling(x);
+            }
+            int signgam = 1;
+            int ip = (int)p;
+            if ((ip&1) == 0) {
+                signgam = -1;
+            }
+            double z = q - p;
+            if (z > 0.5) {
+                p = p + 1;
+                z = q - p;
+            }
+            z = q * Math.sin(Math.PI*z);
+            if (z == 0) {
+                return (signgam == 1) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
+            }
+            z = Math.PI / (Math.abs(z) * stirling(q));
+            return (double)(signgam) * z;
+        }
+
+        // Reduce argument
+        double z = 1.0;
+        while (x >= 3) {
+            x = x - 1;
+            z = z * x;
+        }
+        while (x < 0) {
+            if (x > -1e-09) {
+                return small(x, z);
+            }
+            z = z / x;
+            x = x + 1;
+        }
+        while (x < 2) {
+            if (x < 1e-09) {
+                return small(x, z);
+            }
+            z = z / x;
+            x = x + 1;
+        }
+
+        if (x == 2) {
+            return z;
+        }
+
+        x = x - 2;
+        p = (((((x*_gamP[0]+_gamP[1])*x+_gamP[2])*x+_gamP[3])*x+_gamP[4])*x+_gamP[5])*x + _gamP[6];
+        q = ((((((x*_gamQ[0]+_gamQ[1])*x+_gamQ[2])*x+_gamQ[3])*x+_gamQ[4])*x+_gamQ[5])*x+_gamQ[6])*x + _gamQ[7];
+        return z * p / q;
+    }
+
+
+
+    static final double[] _lgamA = {
+	7.72156649015328655494e-02, // 0x3FB3C467E37DB0C8
+	3.22467033424113591611e-01, // 0x3FD4A34CC4A60FAD
+	6.73523010531292681824e-02, // 0x3FB13E001A5562A7
+	2.05808084325167332806e-02, // 0x3F951322AC92547B
+	7.38555086081402883957e-03, // 0x3F7E404FB68FEFE8
+	2.89051383673415629091e-03, // 0x3F67ADD8CCB7926B
+	1.19270763183362067845e-03, // 0x3F538A94116F3F5D
+	5.10069792153511336608e-04, // 0x3F40B6C689B99C00
+	2.20862790713908385557e-04, // 0x3F2CF2ECED10E54D
+	1.08011567247583939954e-04, // 0x3F1C5088987DFB07
+	2.52144565451257326939e-05, // 0x3EFA7074428CFA52
+	4.48640949618915160150e-05, // 0x3F07858E90A45837
+    };
+
+    static final double[] _lgamR ={
+            1.0, // placeholder
+            1.39200533467621045958e+00, // 0x3FF645A762C4AB74
+            7.21935547567138069525e-01, // 0x3FE71A1893D3DCDC
+            1.71933865632803078993e-01, // 0x3FC601EDCCFBDF27
+            1.86459191715652901344e-02, // 0x3F9317EA742ED475
+            7.77942496381893596434e-04, // 0x3F497DDACA41A95B
+            7.32668430744625636189e-06, // 0x3EDEBAF7A5B38140
+    };
+    static final double[] _lgamS ={
+            -7.72156649015328655494e-02, // 0xBFB3C467E37DB0C8
+            2.14982415960608852501e-01,  // 0x3FCB848B36E20878
+            3.25778796408930981787e-01,  // 0x3FD4D98F4F139F59
+            1.46350472652464452805e-01,  // 0x3FC2BB9CBEE5F2F7
+            2.66422703033638609560e-02,  // 0x3F9B481C7E939961
+            1.84028451407337715652e-03,  // 0x3F5E26B67368F239
+            3.19475326584100867617e-05,  // 0x3F00BFECDD17E945
+    };
+    static final double[] _lgamT ={
+            4.83836122723810047042e-01,  // 0x3FDEF72BC8EE38A2
+            -1.47587722994593911752e-01, // 0xBFC2E4278DC6C509
+            6.46249402391333854778e-02,  // 0x3FB08B4294D5419B
+            -3.27885410759859649565e-02, // 0xBFA0C9A8DF35B713
+            1.79706750811820387126e-02,  // 0x3F9266E7970AF9EC
+            -1.03142241298341437450e-02, // 0xBF851F9FBA91EC6A
+            6.10053870246291332635e-03,  // 0x3F78FCE0E370E344
+            -3.68452016781138256760e-03, // 0xBF6E2EFFB3E914D7
+            2.25964780900612472250e-03,  // 0x3F6282D32E15C915
+            -1.40346469989232843813e-03, // 0xBF56FE8EBF2D1AF1
+            8.81081882437654011382e-04,  // 0x3F4CDF0CEF61A8E9
+            -5.38595305356740546715e-04, // 0xBF41A6109C73E0EC
+            3.15632070903625950361e-04,  // 0x3F34AF6D6C0EBBF7
+            -3.12754168375120860518e-04, // 0xBF347F24ECC38C38
+            3.35529192635519073543e-04,  // 0x3F35FD3EE8C2D3F4
+    };
+    static final double[] _lgamU ={
+            -7.72156649015328655494e-02, // 0xBFB3C467E37DB0C8
+            6.32827064025093366517e-01,  // 0x3FE4401E8B005DFF
+            1.45492250137234768737e+00,  // 0x3FF7475CD119BD6F
+            9.77717527963372745603e-01,  // 0x3FEF497644EA8450
+            2.28963728064692451092e-01,  // 0x3FCD4EAEF6010924
+            1.33810918536787660377e-02,  // 0x3F8B678BBF2BAB09
+    };
+    static final double[] _lgamV ={
+            1.0,
+            2.45597793713041134822e+00, // 0x4003A5D7C2BD619C
+            2.12848976379893395361e+00, // 0x40010725A42B18F5
+            7.69285150456672783825e-01, // 0x3FE89DFBE45050AF
+            1.04222645593369134254e-01, // 0x3FBAAE55D6537C88
+            3.21709242282423911810e-03, // 0x3F6A5ABB57D0CF61
+    };
+    static final double[] _lgamW ={
+            4.18938533204672725052e-01,  // 0x3FDACFE390C97D69
+            8.33333333333329678849e-02,  // 0x3FB555555555553B
+            -2.77777777728775536470e-03, // 0xBF66C16C16B02E5C
+            7.93650558643019558500e-04,  // 0x3F4A019F98CF38B6
+            -5.95187557450339963135e-04, // 0xBF4380CB8C0FE741
+            8.36339918996282139126e-04,  // 0x3F4B67BA4CDAD5D1
+            -1.63092934096575273989e-03, // 0xBF5AB89D0B9E43E4
+    };
+
+    public static double lgamma(double x) {
+        if (Double.isInfinite(x)) {
+            return Double.POSITIVE_INFINITY;
+        }
+
+        if ((0.0 == x) ||
+            (-0.0 == x) ||
+            (-1.0 == x) ||
+            (-2.0 == x)) {
+            throw Py.ValueError(Double.toString(x));
+        }
+
+        double result = _lgamma(x);
+        if (Double.isInfinite(result)) {
+            throw Py.OverflowError(Double.toString(x));
+        }
+
+        return result;
+    }
+
+    public static double _lgamma(double x) {
+	// special cases
+        if (Double.isNaN(x)) {
+            return x;
+        }
+
+        if (Double.isInfinite(x)) {
+            return x;
+        }
+
+        if (0 == x) {
+            return Double.POSITIVE_INFINITY;
+        }
+
+        final double ymin  = 1.461632144968362245;
+        final double two52 = 1 << 52; // 0x4330000000000000 ~4.5036e+15
+        final double two53 = 1 << 53; // 0x4340000000000000 ~9.0072e+15
+        final double two58 = 1 << 58; // 0x4390000000000000 ~2.8823e+17
+        final double tiny  = 1.0 / (1 << 70); // 0x3b90000000000000 ~8.47033e-22
+        final double tc    = 1.46163214496836224576e+00; // 0x3FF762D86356BE3F
+        final double tf    = -1.21486290535849611461e-01; // 0xBFBF19B9BCC38A42
+        // tt = -(tail of tf)
+        final double tt = -3.63867699703950536541e-18; // 0xBC50C7CAA48A971F
+
+
+        int sign = 1;
+        boolean neg = false;
+        if (x < 0) {
+            x = -x;
+            neg = true;
+        }
+
+        if (x < tiny) { // if |x| < 2**-70, return -log(|x|)
+            if (neg) {
+                sign = -1;
+            }
+            return -Math.log(x) * sign;
+        }
+
+        double nadj = 0.0;
+        if (neg) {
+            if (x >= two52) { // |x| >= 2**52, must be -integer
+                return Double.POSITIVE_INFINITY;
+            }
+            double t = sinPi(x);
+            if (t == 0) {
+                return Double.POSITIVE_INFINITY;
+            }
+            nadj = Math.log(Math.PI / Math.abs(t*x));
+            if (t < 0) {
+                sign = -1;
+            }
+        }
+
+        if ((x == 1) || (x == 2)) {
+            return 0.0;
+        }
+
+        double lgamma;
+        double z, p, p1, p2, p3, w, t, y;
+        int i;
+        if (x < 2) {
+            if (x <= 0.9) {
+                lgamma = -Math.log(x);
+                if (x >= (ymin - 1 + 0.27)) { // 0.7316 <= x <=  0.9
+                    y = 1 - x;
+                    i = 0;
+                }
+                else if (x >= (ymin - 1 - 0.27)) { // 0.2316 <= x < 0.7316
+                    y = x - (tc - 1);
+                    i = 1;
+                } else { // 0 < x < 0.2316
+                    y = x;
+                    i = 2;
+                }
+            } else {
+                lgamma = 0;
+                if (x >= (ymin + 0.27)) { // 1.7316 <= x < 2
+                    y = 2 - x;
+                    i = 0;
+                } else if (x >= (ymin - 0.27)) { // 1.2316 <= x < 1.7316
+                    y = x - tc;
+                    i = 1;
+                } else { // 0.9 < x < 1.2316
+                    y = x - 1;
+                    i = 2;
+                }
+            }
+            switch (i) {
+                case 0:
+                    z = y * y;
+                    p1 = _lgamA[0] + z*(_lgamA[2]+z*(_lgamA[4]+z*(_lgamA[6]+z*(_lgamA[8]+z*_lgamA[10]))));
+                    p2 = z * (_lgamA[1] + z*(+_lgamA[3]+z*(_lgamA[5]+z*(_lgamA[7]+z*(_lgamA[9]+z*_lgamA[11])))));
+                    p = y*p1 + p2;
+                    lgamma += (p - 0.5*y);
+                    break;
+                case 1:
+                    z = y * y;
+                    w = z * y;
+                    p1 = _lgamT[0] + w*(_lgamT[3]+w*(_lgamT[6]+w*(_lgamT[9]+w*_lgamT[12]))); // parallel comp
+                    p2 = _lgamT[1] + w*(_lgamT[4]+w*(_lgamT[7]+w*(_lgamT[10]+w*_lgamT[13])));
+                    p3 = _lgamT[2] + w*(_lgamT[5]+w*(_lgamT[8]+w*(_lgamT[11]+w*_lgamT[14])));
+                    p = z*p1 - (tt - w*(p2+y*p3));
+                    lgamma += (tf + p);
+                    break;
+                case 2:
+                    p1 = y * (_lgamU[0] + y*(_lgamU[1]+y*(_lgamU[2]+y*(_lgamU[3]+y*(_lgamU[4]+y*_lgamU[5])))));
+                    p2 = 1 + y*(_lgamV[1]+y*(_lgamV[2]+y*(_lgamV[3]+y*(_lgamV[4]+y*_lgamV[5]))));
+                    lgamma += (-0.5*y + p1/p2);
+                    break;
+            }
+        } else if (x < 8) { // 2 <= x < 8
+		i = (int)x;
+		y = x - (double)i;
+		p = y * (_lgamS[0] + y*(_lgamS[1]+y*(_lgamS[2]+y*(_lgamS[3]+y*(_lgamS[4]+y*(_lgamS[5]+y*_lgamS[6]))))));
+		double q = 1 + y*(_lgamR[1]+y*(_lgamR[2]+y*(_lgamR[3]+y*(_lgamR[4]+y*(_lgamR[5]+y*_lgamR[6])))));
+		lgamma = 0.5*y + p/q;
+		z = 1.0; // Lgamma(1+s) = Log(s) + Lgamma(s)
+                switch (i) {
+                    case 7:
+                        z *= (y + 6);
+                    case 6:
+                        z *= (y + 5);
+                    case 5:
+                        z *= (y + 4);
+                    case 4:
+                        z *= (y + 3);
+                    case 3:
+                        z *= (y + 2);
+                }
+        } else if (x < two58) { // 8 <= x < 2**58
+		t = Math.log(x);
+		z = 1 / x;
+		y = z * z;
+		w = _lgamW[0] + z*(_lgamW[1]+y*(_lgamW[2]+y*(_lgamW[3]+y*(_lgamW[4]+y*(_lgamW[5]+y*_lgamW[6])))));
+		lgamma = (x-0.5)*(t-1) + w;
+        } else { // 2**58 <= x <= Inf
+            lgamma = x * (Math.log(x) - 1);
+        }
+
+        if (neg) {
+            lgamma = nadj - lgamma;
+        }
+
+        return lgamma * sign;
+    }
+
+    public static double sinPi(double x) {
+        final double two52 = 1 << 52; // 0x4330000000000000 ~4.5036e+15
+        final double two53 = 1 << 53; // 0x4340000000000000 ~9.0072e+15
+        if (x < 0.25) {
+            return -Math.sin(Math.PI * x);
+        }
+
+	// argument reduction
+	double z = Math.floor(x);
+	int n;
+        if (z != x) { // inexact
+            x = x % 2.0;
+            n = (int)(x * 4);
+        } else {
+            if (x >= two53) { // x must be even
+                x = 0;
+                n = 0;
+            } else {
+                if (x < two52) {
+                    z = x + two52; // exact
+                }
+                n = (int)(1 & Double.doubleToLongBits(z));
+                x = (double)n;
+                n <<= 2;
+            }
+        }
+        switch (n) {
+            case 0:
+                x = Math.sin(Math.PI * x);
+                break;
+            case 1:
+            case 2:
+                x = Math.cos(Math.PI * (0.5 - x));
+                break;
+            case 3:
+            case 4:
+                x = Math.sin(Math.PI * (1 - x));
+                break;
+            case 5:
+            case 6:
+                x = -Math.cos(Math.PI * (x - 1.5));
+                break;
+            default:
+                x = Math.sin(Math.PI * (x - 2));
+                break;
+        }
+        return -x;
+    } 
+}