scipy: BUG: handle inconsistencies in factorial functions and their extensions

Thanks for the incredibly detailed report @h-vetinari 🤯! I suggest adding a TLDR; something like this if I understand the issue correctly:

1. factorial, factorial2, and factorialk behave inconsistently for several different combinations of inputs. (for example, factorial2 returning 0-dim arrays for scalars and, factorial2 and factorialk aren’t vectorized). See the first table for more details. This issue can be tacked by adding a new parameter like legacy which, by default, preserves the old behaviour but when false, handles the current inconsistent behavior correctly.
2. Support for float and complex inputs can easily be added using this extension. This enhancement adds a new boolean parameter extend that specifies whether to use the extension formula or the default formula. See the second table for more details.

Can this issue be split into two: one for addressing inconsistencies and one for adding the extension?

• factorial2(<scalar>, exact=False) returns a zero-dimensional array instead of a scalar.

This is somewhat common in stats too. Although it’s an inconsistency, such outputs mostly behave like scalars (and also get coerced into one if operated with scalars). I think there are more important things we should tackle first: for example, vectorizing factorial2 and factorialk (these functions throwing for array inputs is very unintuitive).

• As such, the exact keyword only really matters to positive scalar integers; in all other cases, it could be retired, resp. raise.

I agree that raising would be nice but that would again be a backward-incompatible change. Rather, we can add a keyword, similar to nan_policy, that asks the user whether to coerce the inputs to integers, raise, or fill NaNs for invalid cases in array inputs (e.g. exact_policy='raise' | 'coerce' | 'fill-nan'). For factorial, the default could be 'coerce'. For factorial2 and factorialk, raise.

• Complex extensions (together with non-positive floats) should be behind a keyword

This sounds good.

I would propose to introduce another keyword (sorry…) legacy, defaulting to True, but implement legacy=False with all the corrected behaviour. Where there’s any planned change in behaviour, we raises a deprecation warning recommending to switch to legacy=False. After 2 releases, we switch the default from legacy=True to legacy=False, and after two more releases we remove the keyword.

I am not against this but maybe we can also consider an alternative approach like adding exact_policy?

I consider this issue solved by #15841. For the open points regarding the extensions, I’ve opened a follow-up issue: https://github.com/scipy/scipy/issues/18409

Almost in every software I used, if we have a multiple function doing similar but different versions of a functionality 5 years down the line we struggle deprecating/removing them.

I would have liked this function to be in NumPy namespace since it is a very common math operation or from standard math library but it is now living in special functions territory because reasons.

So at least can we combine them into a single signature since they are special cases of k=1, 2 and get rid of the others eventually in the future?

Next problem - the “universal” approximations for n!! is very bad for even numbers (the odd ones work well). For example, the formula (1) from the OP, calculated to high precision with mpmath gives: 10!! = 3063.8767134830029.

However, the true value is very different:

>>> factorial2(10)
3840.

and the current approximation actually does an extremely good job (down to a single eps even for n=100)

>>> factorial2(10) / float(factorial2(10, exact=True)) -1
0.0
>>> factorial2(100) / float(factorial2(100, exact=True)) - 1
2.220446049250313e-16
>>> float(factorial2(100, exact=True))
3.4243224702511973e+79

However, the formula picked up by @j-bowhay in #15349 seems to do a much better job also for even integers.

While working on this, I just stumbled over test_factorial_0d_return_type, and thus https://github.com/scipy/scipy/issues/7400 & https://github.com/scipy/scipy/pull/10954. This shows that the 0d array-output for factorial2 was already undesired a long time ago, and was merely overlooked in #10954.

The reason this is relevant is that if we treat it as a bugfix (similar to #10954), then we’d have no deprecations to do for the other points in this issue.

However, I wanted to tag the participants of those old discussion, because I think #10954 went too far in fixing #7400 - I don’t think the factorial-functions should be changing array-inputs (even zero-dimensional ones) to scalars. IMO it should be exclusively:

• scalar->scalar
• array->array (with identical dimensions).

CC @WarrenWeckesser @person142 @rgommers @larsoner

Maybe a good first step would be to add some benchmarks?

Absolutely! 😃

IMO we should aim basically for all the cases in the first table that don’t raise (except NaN/empty array), and some mixed inputs in the array case (though the negative inputs returning 0 are presumably less interesting).

But to do that, we need some benchmarks against the current most important code paths (haven’t verified current status of factorial in benchmarks/), and ideally some as-tight-as-possible accuracy tests.

Agreed, this was my hesitation in making them thin wrappers as there is a potential loss in performance and accuracy as some cases have more optimal ways of performing the calculation that don’t generalise.

In benchmarks/benchmarks/special.py there are currently no benchmarks for any of the factorial stuff. Maybe a good first step would be to add some benchmarks?

@tirthasheshpatel

Thanks for the incredibly detailed report @h-vetinari 🤯!

🙃

I suggest adding a TLDR; something like this if I understand the issue correctly:

Done, thanks!

Can this issue be split into two: one for addressing inconsistencies and one for adding the extension?

Fine by me, and effectively adopted with your TL;DR.

• factorial2(<scalar>, exact=False) returns a zero-dimensional array instead of a scalar.

This is somewhat common in stats too. Although it’s an inconsistency, such outputs mostly behave like scalars (and also get coerced into one if operated with scalars). I think there are more important things we should tackle first: for example, vectorizing factorial2 and factorialk (these functions throwing for array inputs is very unintuitive).

Prioritization aside, I think we should clean this up. 0-dim arrays used to be a bit of a wild west thing, but it’s getting cleaned up in lots of APIs (not just numpy/scipy, also other downstreams like cvxpy).

• As such, the exact keyword only really matters to positive scalar integers; in all other cases, it could be retired, resp. raise.

I agree that raising would be nice but that would again be a backward-incompatible change. Rather, we can add a keyword, similar to nan_policy, that asks the user whether to coerce the inputs to integers, raise, or fill NaNs for invalid cases in array inputs (e.g. exact_policy='raise' | 'coerce' | 'fill-nan'). For factorial, the default could be 'coerce'. For factorial2 and factorialk, raise.

I’m not a fan, to be honest. It has the same high price (another new keyword) for very little gain IMO. In particular, I think coerce is basically never the right thing, and if so, then people should get the warning/error and explicitly cast.

• Complex extensions (together with non-positive floats) should be behind a keyword

This sounds good.

We’ll have to do some bikeshedding on this. I proposed extend=True|False (more or less as a stand-in, because extend_complex seemed too long), but if we want to leave the recurrence thing on the table, we could actually do extend='zero'|'complex'|'recurrence' (with 'zero' being the default and 'recurrence' potentially being added later), which to me looks even better than boolean switches (even if 'recurrence' never materializes)

I would propose to introduce another keyword (sorry…) legacy, defaulting to True, but implement legacy=False with all the corrected behaviour. Where there’s any planned change in behaviour, we raises a deprecation warning recommending to switch to legacy=False. After 2 releases, we switch the default from legacy=True to legacy=False, and after two more releases we remove the keyword.

I am not against this but maybe we can also consider an alternative approach like adding exact_policy?

I understand that it’s not attractive, but I had a closer look and added this to the OP:

Looking more closely at relevant changes in (currently non-raising) behaviour, it really only affects the zero-dimensional arrays for factorial2(scalar_val, exact=False). In particular, we could choose to let exact=True only raise if extend=True as well, keeping essentially all current non-raising behaviour intact.

In other words, that’s the one behavioural change that we cannot avoid anyway, and using it for that aspect is IMO a reasonable transition (and only one affected function). We could still deprecate certain combinations that are currently non-raising by deprecating them, also with IMO acceptably little impact.