Since we're not in a hurry for Py2.7 and 3.1, I would like to this 
kicked around a bit on the newsgroup and in numpy forums (personally, I 
would also post a pure python equivalent to the ASPN cookbook for 
further commentary).

There are several different approximations to choose from.  Each of 
them has their own implications for speed and accuracy.  IIRC, the one 
I used in test.test_random.gamma() was accompanied by a proof that its 
maximum error never exceeded a certain amount.  I think there were some 
formulas that made guarantees only over a certain interval and others 
that had nice properties in the first and second derivatives (one that 
don't have those properties can throw newtonian solvers wildly off the 
mark). 

Let's let the community use its collective wisdom to select the best 
approach and not immediately commit ourselves to the one in this patch.

At one point, Tim was reluctant to add any of these functions because 
it is non-trivial to do well and because it would open a can of worms 
about why python gets a different answer (possibly not as close to 
true) and some other favorite tool (matlab, excel, numpy, and 
engineering calculator, etc).

FWIW, I changed the 2.6 version of test_random.gamma() to take 
advantage of msum() to increase its accuracy when some of the summands 
have nearly cancelling values.