Move automatic simplifications of gamma into a calculate function - #2123

This PR defers calculations of the gamma function from the automatic simplification to an explicit function call. It ties in to the discussion in #2056, but instead of arguing why this is a good idea in general I am here trying to optimize the use of the gamma function in some real applications.

Apologies for the long text, but since this has been a point of discussion before, I would like to make sure that the reasoning behind this is clear to everyone.

Before the PR

Before this PR the gamma function had automatic simplifications for:

• Integers
• Rationals of the form p / 2
• Reals

What this PR does

• Removes all automatic simplifications for the gamma function
• Adds the calculate function to calculate gamma for these cases (and possible future cases)

calculate (or if we prefer some other name) would be like a doit function that performs "deferred" calculations. This is then distinct from eval that evaluates expressions onto some floating point number and refine that leverage assumptions to simplify expressions.

Benchmarks

The gamma function can be very slow to calculate. I did some benchmarking of the current situation using symengine.py on my laptop:

Integers

It looks like at least O(n) in the number of digits for integers.

Rationals

This looks like O(n) for number of digits in the numerator.

RealDouble

Without auto simplification

Couldn't get good measurements, but should be somewhere around 20ns - 100ns.

Why?

There are arguments against automatic simplification in general (see #2056). I will give arguments for this particular case.

Result might not be needed

Always calculating the gamma for integers takes time and space to store the often large result and the result might not be needed. Some examples:

• We have a Piecewise like: Piecewise(x > 0, gamma(8977887798), 0). If doing subs with x=-1 the calculation of the gamma was useless.
• We have gamma(8387889) * n where I do subs with n=0
• We want to know if an expression containing gamma is positive. We do not need to calculate the gamma to know that.
• We want to know if an expression containing gamma is a real value.
• For high enough arguments to gamma the result might not fit in RAM and a crash will occur. If the results are not really needed we can avoid a crash in these situations.

There might be faster ways

If we have the expression 1001*gamma(1001) instead of calculating gamma(1001) and do the multiplication we could leverage that $\Gamma(z + 1) = z\Gamma(z)$ and instead calculate gamma(1002) removing the need for the extra multiplication.

If we have the expression gamma(1001)*gamma(-1000) we could use that $\Gamma(z) \Gamma(1 - z) = \frac{\pi}{\sin{\pi z}}$ directly.

If we have gamma in an inequality for example Ge(gamma(2392), 2) we could potentially evaluate the inequality faster without having to calculate the gamma function.

All these cases might arise from subs situations.

We might want higher precision

If we have gamma(2.35) it will be automatically calculated as RealDouble making it impossible to evaluate with a higher precision later.

Smaller memory footprint

For gamma(237) we would instead of directly calculating and storing the resulting integer we will now instead keep gamma(237) which has a slightly more complex expression tree, but takes less memory since the resulting integer would be much larger (460 decimal digits).

No drawbacks?

One drawback could be if the only thing you want is to calculate gamma for some number now you would need to do calculate(gamma(2393)).

Is this a slippery slope?

I the discussion of #2056 we could see pros and cons of automatic simplification in general. I think this PR still fits under the current system, since as I have argued, there are few benefits of automatic simplification in this particular case, but this must be taken on a case by case basis. I also think this way of thinking can be extended to more cases in functions.cpp, but not into more core operators like mul, add etc.