sympy/sympy — pull request #1001 by ness01 (2012-01-19 14:42:45+00:00)
Title: Fix issue 3016.

1) Certain substitutions may result in exceptions (akin to 0**0).
Wrap such a subs in an exception handler and fail gracefully in meijerint
code _rewrite_single.

2) Add tests for beta distribution.

3) Make hyperexpand turn hyper([..], [..], 0) into 0.
(We don't want to autoevaluate, but simplify should do this.)

4) Add a test for this as well.

Some discussion
==============

Consider the integrand `betadist = x**(a - 1)*(-x + 1)**(b - 1)*gamma(a + b)/(gamma(a)*gamma(b))`, where `a, b` are declared positive. We wish to compute (e.g.) `integrate(betadist, (x, 0, 1))`. There are three approaches, tried in this order:

1. Find an antiderivative with heurisch, evalute it over the limits.
2. Find an antiderivative with meijerint, evaluate it over the limits.
3. Compute the definite integral with meijerint.

As explained in the `integrate()` docstring, the reason for trying things in this order is heuristical, but it usually does a good job. However, in this case, the order is exactly wrong!

(3) Works like a charm.

(2) Works, with the following caveat: the antiderivative is in terms of a hypergeometric function which cannot be expanded. When evaluating it over the limits `(0, 1)` it turns into two terms which *can* be computed (because it vanishes at zero, and at one there is a formula due to gauss). The `integrate()` code does not call hyperexpand, and so the result is expressed in terms of hypergeometric functions even though it is elementary.

(1) Hangs.

There are thus two problems: the cosmetic one that even when calling `integrate(betadist, (x,0,1), meijerg=True)` the result will *not* look elementary (however calling simplify resolves this); and the somewhat more important problem that without the parameter meijerg the integral does not seem to finish in any sensible amount of time.

I would argue that both of these behaviours are OK. Namely, having to call simplify is not a real problem. My code tries to hide hypergeometric terms, but if the most obvious thing to get rid of them works, that's fine. Regarding the hanging integral: I see no better heuristic. A hanging integral is probably something every user of sympy has seen, we just have to educate them that passing `meijerg=True` may sometimes (rather seldomly) alleviate the problem.

So much for the woffle. If you have better suggestions, I'm of course happy to implement them.