Issue 1:
It's currently not possible to compute the full eigendecomposition with eigsh.

>>> A = rand(5, 4)
>>> AA = np.dot(A.T, A)
>>> from scipy.sparse.linalg import eigsh
>>> eigsh(AA, AA.shape[0])
[...]
ValueError: k must be between 1 and the order of the square input matrix.

https://github.com/scipy/scipy/blob/master/scipy/sparse/linalg/eigen/arpack/arpack.py#L1506

I think k >= n should be replaced by k > n.

Issue 2:
In the context of #15351 and reviewing #15299, I noticed a bunch of inconsistencies that should be fixed IMO.

TL;DR

Courtesy of @tirthasheshpatel:

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.

Support for float and complex inputs can easily be added using this extension. This enhancement adds a new boolean parameter extend='zero'|'complex' (and perhaps 'recurrence' down the line) that specifies whether to use the extension formula or the default formula. See the second table for more details.

Design Space

There are several dimensions that come into play:

scalar vs. array inputs

for integers: approximation or exact computation (exact=)

extensions to negative integers

extensions to floats

extensions to complex numbers

extensions to poles (possible e.g. for factorial2)

for the array-case, failing soft (e.g. inserting NaNs) or hard (raising error)

Current Behaviour (as of 1.8.0)

Inputs
factorial(
exact=False)
factorial(
exact=True)
factorial2(
exact=False)
factorial2(
exact=True)
factorialk(
exact=False)
factorialk(
exact=True)

Scalar positive integer
float
int
!! 0-dim
float array !!
int
NotImpl Error
int

Scalar negative integer
0.0
0
!! array(0.) !!
0
NotImpl Error
0

Scalar positive float
float
ValueError (but text misleading)
TypeError
TypeError
NotImpl Error
TypeError

Scalar negative float
0.0
0
TypeError
TypeError
NotImpl Error
TypeError

Scalar complex
TypeError
TypeError
TypeError
TypeError
NotImpl Error
TypeError

Scalar NaN
NaN
NaN
TypeError
TypeError
NotImpl Error
TypeError

Positive integer in array
float
int
float
!! TypeError !!
NotImpl Error
TypeError

Negative integer in array
0.0
0
0.0
!! TypeError !!
NotImpl Error
TypeError

Positive float in array
float
ValueError (but text misleading)
TypeError
TypeError
NotImpl Error
TypeError

Negative float in array
0.0
0
TypeError
TypeError
NotImpl Error
TypeError

Complex in array
TypeError
TypeError
TypeError
TypeError
NotImpl Error
TypeError

NaN in array
NaN
NaN
TypeError
TypeError
NotImpl Error
TypeError

Empty array
array([])
IndexError
TypeError
TypeError
NotImpl Error
TypeError

Note that:

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

All the TypeErrors are IMO implementation oversights (and in case of "!!", against documentation) - i.e. we should provide better error messages.

There is a full generalization of factorialk which is basically just a fraction of two gamma functions, which would work for all complex n & k (up to poles), as well as arrays.

Extension Issues

In addition to this hodgepodge of behaviour, some people would like to have the float/complex extensions of these functions (e..g. #15349). I tend to agree, but that raises another problem:

Complex extensions are incompatible with current inputs ≤0

First off, everything that's 0 currently on the negative axis gets a different value. But also note the "≤"; namely, 0!! is 1 in the regular definition, but ~0.79=sqrt(2/pi) in the extension.

Extensions based on recurrence are incompatible with current status as well as complex extension

This came up for example in #15299. It's possible (for example) to invert the recurrence relation n!! = n * (n-2)!!, and so extend the double-factorial to negative integers. The resulting value would be neither 0 (currently) nor the pole in the complex extension.

Design Considerations

If we (or at least I) were designing from scratch:

Changing from the exact=False to exact=True should never fundamentally change the result (e.g. going from 0 to non-zero, or casting inputs to integers as seen in #15351).

As such, the exact keyword only really matters to positive scalar integers; in all other cases, it could be retired, resp. raise, or simply be turned into a no-op (cf. #15359)

Open question whether it's worth to properly implement exact=True for array-case (currently, only factorial does).

Complex extensions (together with non-positive floats) should be behind a keyword, roughly because:

0-for-negative-values is a valid approach, and has been in scipy ~forever.

Since the values being produced are incompatible (for certain inputs), the complex extension should be explicitly opt-in.

Using recurrence relations for extending into poles is an open question; not sure if worth it, but if so, probably with another keyword value for the extend keyword.

Use the most general approximation, except when alternatives with better accuracy / performance exist.

Probably need to think about dtype casting rules for empty arrays (keep dtype? always float? ...)

Possible Future Behaviour (note extend=, not exact=)

Inputs
factorial(
extend=
'zero')
factorial(
extend=
'complex')
factorial2(
extend=
'zero')
factorial2(
extend=
'complex')
factorialk(
extend=
'zero')
factorialk(
extend=
'complex')

Scalar positive integer w/ exact=True
n!
n!
n!!
n!!
n!(k)
n!(k)

Scalar positive integer
Γ(n+1)
Γ(n+1)
(1)
(1)
(2)
(2)

Scalar negative integer / poles
0.0
NaN
0.0
NaN
0.0
NaN

Scalar positive float
Γ(n+1)
Γ(n+1)
(1)
(1)
(2)
(2)

Scalar negative float
0.0
Γ(n+1)
0.0
(1)
0.0
(2)

Scalar complex
ValueError
Γ(n+1)
ValueError
(1)
ValueError
(2)

Scalar NaN
NaN
NaN
NaN
NaN
NaN
NaN

Positive integer in array
Γ(n+1)
Γ(n+1)
(1)
(1)
(2)
(2)

Negative integer in array
0.0
NaN
0.0
NaN
0.0
NaN

Positive float in array
Γ(n+1)
Γ(n+1)
(1)
(1)
(2)
(2)

Negative float in array
0.0
Γ(n+1)
0.0
(1)
0.0
(2)

Complex in array
ValueError
Γ(n+1)
ValueError
(1)
ValueError
(2)

NaN in array
NaN
NaN
NaN
NaN
NaN
NaN

Empty array
array([])
array([])
array([])
array([])
array([])
array([])

The formulas in question are:

(1): z!! = 2^((z-1)/2) Γ(z/2 + 1) / Γ(1/2 + 1)

(2): z!(α) = α^((z-1)/α) Γ(z/α + 1) / Γ(1/α + 1)

I'd very much welcome feedback on this - there are some trade-offs to be made (e.g. whether negative floats should also raises without extend='complex' or whether the errors in the array-case should just map to NaN; also there are different approximations for n!! in use currently), but I think the above would be reasonable. In particular, the scalar and array case behave the same way.

Implementation Plan

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, after two more releases we raise on legacy=True, and after two more we remove the keyword.

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, if we let exact=True be a no-op for non-integers (the less "noisy" options from the point of UX), we would keep essentially all current non-raising behaviour intact.

CC @steppi @tirthasheshpatel @j-bowhay @tupui

Edit History:

v4 (2022-03-10): Make exact=True a no-op for non-integers (less noisy), following #15359. Extend deprecation period of legacy-keyword along the same lines as proposed in #15722.

v3 (2022-02-16; up to rev 20): Added (current & proposed) behaviour for NaN and empty arrays; good catch @tirthasheshpatel!

v2 (2022-02-16; up to rev 15): Added TL;DR, clarification on impact of legacy, changed extend from boolean argument to take an enum.

v1 (2022-02-16; up to rev 10): Original