ENH: extensions of `factorial{,2,k}` for complex domains, recurrences, ...
Following the clean-up of #15600 in #15841, we can now discuss the extensions separately.

## Design Space

There are several dimensions that come into play:
- extensions to negative integers
- extensions to floats
- extensions to complex numbers
- extensions to poles through recurrence (possible e.g. for `factorial2`)

## Extension Issues

Some people would like to have the float/complex extensions of these functions (e..g. https://github.com/scipy/scipy/pull/15349). I tend to agree, but that raises another problem:

#### Complex extensions are incompatible with current inputs

First off, everything that's 0 currently on the negative axis gets a different value. More importantly, all even numbers get smaller by a factor ~0.79=sqrt(2/pi) in the [extension](https://en.wikipedia.org/wiki/Double_factorial#Complex_arguments).

For non-complex floats, there doesn't seem to be another relevant extension, especially not one that matches on all the integers and isn't just a gratuitous hack (this is what #15841 ended up stalling on for a long while). Compared to #15600 I now propose that floats continue to be treated as errors for `extend='zero'` in `factorial{2,k}`.

#### Extensions based on recurrence are incompatible with current status as well as complex extension

This came up for example in #15299. It's [possible](https://en.wikipedia.org/wiki/Double_factorial#Negative_arguments) (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

- 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 negative & even integers), 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 enum value for a potential `extend` keyword.
- Use the most general approximation, except when alternatives with better accuracy / performance exist.

## Possible Future Behaviour (note `extend=`, not `exact=`)

| Inputs | `factorial(`<br/>`extend=`<br/>`'zero')` | `factorial(`<br/>`extend=`<br/>`'complex')` | `factorial2(`<br/>`extend=`<br/>`'zero')` | `factorial2(`<br/>`extend=`<br/>`'complex')` | `factorialk(`<br/>`extend=`<br/>`'zero')` | `factorialk(`<br/>`extend=`<br/>`'complex')` |
|-|-|-|-|-|-|-|
| Scalar positive integer&nbsp;w/ `exact=True` | n! | n! | n!! | n!! | n!(k) | n!(k) |
| Scalar positive integer | Γ(n+1) | Γ(n+1) | [1a] for odd,<br/>[1b] for even<br/>(matches n!!) | [1a] (_doesn't_<br/>match n!!<br/>for even n) | n!(k)<br/>(no approx.) | [2] |
| Scalar negative integer / poles | 0.0 | NaN | 0.0 | NaN | 0.0 | NaN |
| Scalar positive float | Γ(n+1) | Γ(n+1) | `ValueError` | [1a] | `ValueError` | [2] |
| Scalar negative float | 0.0 | Γ(n+1) | `ValueError` | [1a] | `ValueError` | [2] |
| Scalar complex | `ValueError` | Γ(n+1) | `ValueError` | [1a] | `ValueError` | [2] |
| Scalar NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| Positive integer in array | Γ(n+1) | Γ(n+1) |[1a] for odd,<br/>[1b] for even<br/>(matches n!!) | [1a] (_doesn't_<br/>match n!!<br/>for even n) | n!(k)<br/>(no approx.) | [2] |
| Negative integer in array | 0.0 | NaN | 0.0 | NaN | 0.0 | NaN |
| Positive float in array | Γ(n+1) | Γ(n+1) | `ValueError` | [1a] | `ValueError` | [2] |
| Negative float in array | 0.0 | Γ(n+1) | `ValueError` | [1a] | `ValueError` | [2] |
| Complex in array | `ValueError` | Γ(n+1) | `ValueError` | [1a] | `ValueError` | [2] |
| NaN in array | NaN | NaN | NaN | NaN | NaN | NaN |
| Empty array | `array([])` | `array([])` | `array([])` | `array([])` | `array([])` | `array([])` |

The formulas in question are:
* [1a]: `z!! = 2^((z-1)/2) Γ(z/2 + 1) / Γ(1/2 + 1) == 2^(z/2) Γ(z/2 + 1) sqrt(2 / π)` (holds exactly for odd integers and is the basis for the extension)
* [1b]: `z!! = 2^(z/2) Γ(z/2 + 1)` (holds only for even integers)
* [2]: `z!(α) =  α^((z-1)/α) Γ(z/α + 1) / Γ(1/α + 1)` (extension based on α'th integers)

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'`, etc.)