Repository: scipy/scipy
Issue category: Feature Request

Issue:
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.)



Hints:
Not a general comment, but just a remark on one specific case.

In SciPy-1.10.*, the result of `factorial2(-1)` was `1` and this changed to `0` in version 1.11, which was an unexpected change that broke a library I'm involved in. We relied on the old behavior for the implementation of integrals involving products of Gaussian functions. In the recurrence relations for these integrals, `factorial2(2*n-1)` often appears in denominators, where `n` can be zero or a positive integer.

In case of the double factorial, it is common to define a meaningful result for negative odd integer arguments, see e.g. https://en.wikipedia.org/wiki/Double_factorial#Negative_arguments. Is there any chance of supporting odd negative arguments in `factorial2`? (Note that this is not related to complex arguments, but rather comes from a generalization of the standard recurrence relation.)




> Is there any chance of supporting odd negative arguments in `factorial2`?

Yes it's possible; I mentioned this in the OP already:🙃 

> 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.

Personally I was waiting for 1.12 to branch before touching the factorial functions again, but feel free to raise a PR!
Sorry for having completely missed your point. I should have gone for that coffee first. The reason I didn't get it, is that in the rows "Scalar negative integer / poles" and "Negative integer in array" contain 0.0 or NaN. Shouldn't these follow the same rules? Letting the outcome depend on the dtype of the argument may cause a lot of confusion.

Do you suggest a PR for a 1.11.x version in which support for negative odd arguments is added to `factorial2` without changing the API? Or do you prefer to go for the API change right away?

P.S. I'm not convinced by the backward compatibility argument to treat integers differently. With 1.11, backward compatibility is already broken.
What I'm suggesting - as described in the OP (so you picked the right issue!) - is to add a keyword argument `extend=` to `factorial2` that allows _opting into_ a different behaviour for `n<0`. In your case, using the values that arise from the recurrence relation. Independently of that, we'll also want to eventually allow using the complex extension, which the design allows without further overhauls.

> P.S. I'm not convinced by the backward compatibility argument to treat integers differently. With 1.11, backward compatibility is already broken.

Yes, it was breakage we accepted (perhaps wrongly without a deprecation period), see https://github.com/scipy/scipy/issues/18813
Thanks for clarifying, and sorry for having missed some points. Just to avoid confusion, I'd rather ask a few more questions:
- The proposal is to have `extend="zero"` as the default value, because this is the current behavior?
- For our use case, the option `extend="recurrence"` would need to be implemented in `factorial2`? (It's not in the table, but this seems to make sense.)
- The `extend="recurrence"` is only meaningful for (some) integer arguments of `factorial2`. Would you return `NaN` when it is combined with a `float` or `complex` argument, or when a negative even integer argument is given?
> The proposal is to have `extend="zero"` as the default value, because this is the current behavior?

Yes

> For our use case, the option `extend="recurrence"` would need to be implemented in `factorial2`?

Yes

> The `extend="recurrence"` is only meaningful for (some) integer arguments of `factorial2`. Would you return `NaN` when it is combined with a `float` or `complex` argument, or when a negative even integer argument is given?

Raising an error for known-to-be-incompatible values is the best approach IMO.

Patch:
diff --git a/scipy/special/_basic.py b/scipy/special/_basic.py
index 103f2c867632..5fe7f17e2784 100644
--- a/scipy/special/_basic.py
+++ b/scipy/special/_basic.py
@@ -11,7 +11,7 @@
                    sin, place, issubdtype, extract, inexact, nan, zeros, sinc)
 
 from . import _ufuncs
-from ._ufuncs import (mathieu_a, mathieu_b, iv, jv, gamma,
+from ._ufuncs import (mathieu_a, mathieu_b, iv, jv, gamma, rgamma,
                       psi, hankel1, hankel2, yv, kv, poch, binom,
                       _stirling2_inexact)
 
@@ -2911,34 +2911,68 @@ def _factorialx_array_exact(n, k=1):
     return out
 
 
-def _factorialx_array_approx(n, k):
+def _factorialx_array_approx(n, k, extend):
     """
     Calculate approximation to multifactorial for array n and integer k.
 
-    Ensure we only call _factorialx_approx_core where necessary/required.
+    Ensure that values aren't calculated unnecessarily.
     """
+    if extend == "complex":
+        return _factorialx_approx_core(n, k=k, extend=extend)
+
+    # at this point we are guaranteed that extend='zero' and that k>0 is an integer
     result = zeros(n.shape)
     # keep nans as nans
     place(result, np.isnan(n), np.nan)
     # only compute where n >= 0 (excludes nans), everything else is 0
     cond = (n >= 0)
     n_to_compute = extract(cond, n)
-    place(result, cond, _factorialx_approx_core(n_to_compute, k=k))
+    place(result, cond, _factorialx_approx_core(n_to_compute, k=k, extend=extend))
     return result
 
 
-def _factorialx_approx_core(n, k):
+def _gamma1p(vals):
+    """
+    returns gamma(n+1), though with NaN at -1 instead of inf, c.f. #21827
+    """
+    res = gamma(vals + 1)
+    # replace infinities at -1 (from gamma function at 0) with nan
+    # gamma only returns inf for real inputs; can ignore complex case
+    if isinstance(res, np.ndarray):
+        if not _is_subdtype(vals.dtype, "c"):
+            res[vals == -1] = np.nan
+    elif np.isinf(res) and vals == -1:
+        res = np.float64("nan")
+    return res
+
+
+def _factorialx_approx_core(n, k, extend):
     """
     Core approximation to multifactorial for array n and integer k.
     """
     if k == 1:
-        # shortcut for k=1
-        result = gamma(n + 1)
+        # shortcut for k=1; same for both extensions, because we assume the
+        # handling of extend == 'zero' happens in _factorialx_array_approx
+        result = _gamma1p(n)
         if isinstance(n, np.ndarray):
-            # gamma does not maintain 0-dim arrays
+            # gamma does not maintain 0-dim arrays; fix it
             result = np.array(result)
         return result
 
+    if extend == "complex":
+        # see https://numpy.org/doc/stable/reference/generated/numpy.power.html
+        p_dtype = complex if (_is_subdtype(type(k), "c") or k < 0) else None
+        with warnings.catch_warnings():
+            # do not warn about 0 * inf, nan / nan etc.; the results are correct
+            warnings.simplefilter("ignore", RuntimeWarning)
+            result = np.power(k, (n - 1) / k, dtype=p_dtype) * _gamma1p(n / k)
+            result *= rgamma(1 / k + 1)
+        if isinstance(n, np.ndarray):
+            # ensure we keep array-ness for 0-dim inputs; already n/k above loses it
+            result = np.array(result)
+        return result
+
+    # at this point we are guaranteed that extend='zero' and that k>0 is an integer
     n_mod_k = n % k
     # scalar case separately, unified handling would be inefficient for arrays;
     # don't use isscalar due to numpy/numpy#23574; 0-dim arrays treated below
@@ -2988,7 +3022,7 @@ def _is_subdtype(dtype, dtypes):
     return any(np.issubdtype(dtype, dt) for dt in dtypes)
 
 
-def factorial(n, exact=False):
+def factorial(n, exact=False, extend="zero"):
     """
     The factorial of a number or array of numbers.
 
@@ -3006,6 +3040,10 @@ def factorial(n, exact=False):
         If False, result is approximated in floating point rapidly using the
         `gamma` function.
         Default is False.
+    extend : string, optional
+        One of ``'zero'`` or ``'complex'``; this determines how values ``n<0``
+        are handled - by default they are 0, but it is possible to opt into the
+        complex extension of the factorial (the Gamma function).
 
     Returns
     -------
@@ -3036,47 +3074,64 @@ def factorial(n, exact=False):
     120
 
     """
+    if extend not in ("zero", "complex"):
+        raise ValueError(
+            f"argument `extend` must be either 'zero' or 'complex', received: {extend}"
+        )
+    if exact and extend == "complex":
+        raise ValueError("Incompatible options: `exact=True` and `extend='complex'`")
+
+    msg_needs_complex = (
+        "In order to use non-integer arguments, you must opt into this by passing "
+        "`extend='complex'`. Note that this changes the result for all negative "
+        "arguments (which by default return 0)."
+    )
     msg_wrong_dtype = (
         "Unsupported data type for factorial: {dtype}\n"
-        "Permitted data types are integers and floating point numbers."
+        "Permitted data types are integers and floating point numbers, "
+        "as well as complex numbers if `extend='complex'` is passed."
+    )
+    msg_exact_not_possible = (
+        "`exact=True` only supports integers, cannot use data type {dtype}"
     )
 
     # don't use isscalar due to numpy/numpy#23574; 0-dim arrays treated below
     if np.ndim(n) == 0 and not isinstance(n, np.ndarray):
         # scalar cases
-        if not _is_subdtype(type(n), ["i", "f", type(None)]):
+        if not _is_subdtype(type(n), ["i", "f", "c", type(None)]):
             raise ValueError(msg_wrong_dtype.format(dtype=type(n)))
+        elif _is_subdtype(type(n), "c") and extend != "complex":
+            raise ValueError(msg_needs_complex)
         elif n is None or np.isnan(n):
-            return np.float64("nan")
-        elif n < 0:
+            complexify = (extend == "complex") and _is_subdtype(type(n), "c")
+            return np.complex128("nan+nanj") if complexify else np.float64("nan")
+        elif extend == "zero" and n < 0:
             return 0 if exact else np.float64(0)
         elif exact and _is_subdtype(type(n), "i"):
             return math.factorial(n)
         elif exact:
-            msg = ("Non-integer values of `n` together with `exact=True` are "
-                   "not supported. Either ensure integer `n` or use `exact=False`.")
-            raise ValueError(msg)
-        return _factorialx_approx_core(n, k=1)
+            raise ValueError(msg_exact_not_possible.format(dtype=type(n)))
+        return _factorialx_approx_core(n, k=1, extend=extend)
 
     # arrays & array-likes
     n = asarray(n)
 
-    if not _is_subdtype(n.dtype, ["i", "f"]):
+    if not _is_subdtype(n.dtype, ["i", "f", "c"]):
         raise ValueError(msg_wrong_dtype.format(dtype=n.dtype))
-    if exact and not _is_subdtype(n.dtype, "i"):
-        msg = ("factorial with `exact=True` does not "
-               "support non-integral arrays")
-        raise ValueError(msg)
+    elif _is_subdtype(n.dtype, "c") and extend != "complex":
+        raise ValueError(msg_needs_complex)
+    elif exact and _is_subdtype(n.dtype, ["f", "c"]):
+        raise ValueError(msg_exact_not_possible.format(dtype=n.dtype))
 
     if n.size == 0:
         # return empty arrays unchanged
         return n
     elif exact:
         return _factorialx_array_exact(n, k=1)
-    return _factorialx_array_approx(n, k=1)
+    return _factorialx_array_approx(n, k=1, extend=extend)
 
 
-def factorial2(n, exact=False):
+def factorial2(n, exact=False, extend="zero"):
     """Double factorial.
 
     This is the factorial with every second value skipped.  E.g., ``7!! = 7 * 5
@@ -3086,6 +3141,8 @@ def factorial2(n, exact=False):
           = 2 ** (n / 2) * gamma(n / 2 + 1)                 n even
           = 2 ** (n / 2) * (n / 2)!                         n even
 
+    The formula for `n odd` is the basis for the complex extension.
+
     Parameters
     ----------
     n : int or array_like
@@ -3094,6 +3151,17 @@ def factorial2(n, exact=False):
         The result can be approximated rapidly using the gamma-formula
         above (default).  If `exact` is set to True, calculate the
         answer exactly using integer arithmetic.
+    extend : string, optional
+        One of ``'zero'`` or ``'complex'``; this determines how values ``n<0``
+        are handled - by default they are 0, but it is possible to opt into the
+        complex extension of the double factorial. This also enables passing
+        complex values to ``n``.
+
+        .. warning::
+
+           Using the ``'complex'`` extension also changes the values of the
+           double factorial for even integers, reducing them by a factor of
+           ``sqrt(2/pi) ~= 0.79``, see [1].
 
     Returns
     -------
@@ -3109,41 +3177,68 @@ def factorial2(n, exact=False):
     >>> factorial2(7, exact=True)
     105
 
+    References
+    ----------
+    .. [1] Complex extension to double factorial
+            https://en.wikipedia.org/wiki/Double_factorial#Complex_arguments
     """
+    if extend not in ("zero", "complex"):
+        raise ValueError(
+            f"argument `extend` must be either 'zero' or 'complex', received: {extend}"
+        )
+    if exact and extend == "complex":
+        raise ValueError("Incompatible options: `exact=True` and `extend='complex'`")
+
+    msg_needs_complex = (
+        "In order to use non-integer arguments, you must opt into this by passing "
+        "`extend='complex'`. Note that this changes the result for all negative "
+        "arguments (which by default return 0). Additionally, it will rescale the "
+        "values of the double factorial at even integers by a factor of sqrt(2/pi)."
+    )
     msg_wrong_dtype = (
         "Unsupported data type for factorial2: {dtype}\n"
-        "Only integers are permitted."
+        "Only integers are permitted by default, though floating point "
+        "and complex numbers can be used if `extend='complex'` is passed."
     )
 
     # don't use isscalar due to numpy/numpy#23574; 0-dim arrays treated below
     if np.ndim(n) == 0 and not isinstance(n, np.ndarray):
         # scalar cases
-        if not _is_subdtype(type(n), "i"):
+        if not _is_subdtype(type(n), ["i", "f", "c", type(None)]):
             raise ValueError(msg_wrong_dtype.format(dtype=type(n)))
-        elif n < 0:
-            return 0
+        elif _is_subdtype(type(n), ["f", "c"]) and extend != "complex":
+            raise ValueError(msg_needs_complex)
+        elif n is None or np.isnan(n):
+            complexify = (extend == "complex") and _is_subdtype(type(n), "c")
+            return np.complex128("nan+nanj") if complexify else np.float64("nan")
+        elif extend == "zero" and n < 0:
+            return 0 if exact else np.float64(0)
         elif n in {0, 1}:
-            return 1
-        # general integer case
+            return 1 if exact else np.float64(1)
+
         if exact:
+            # general integer case
             return _range_prod(1, n, k=2)
-        return _factorialx_approx_core(n, k=2)
+        # approximation
+        return _factorialx_approx_core(n, k=2, extend=extend)
 
     # arrays & array-likes
     n = asarray(n)
 
-    if not _is_subdtype(n.dtype, "i"):
+    if not _is_subdtype(n.dtype, ["i", "f", "c"]):
         raise ValueError(msg_wrong_dtype.format(dtype=n.dtype))
+    elif _is_subdtype(n.dtype, ["f", "c"]) and extend != "complex":
+        raise ValueError(msg_needs_complex)
 
     if n.size == 0:
         # return empty arrays unchanged
         return n
     elif exact:
         return _factorialx_array_exact(n, k=2)
-    return _factorialx_array_approx(n, k=2)
+    return _factorialx_array_approx(n, k=2, extend=extend)
 
 
-def factorialk(n, k, exact=None):
+def factorialk(n, k, exact=None, extend="zero"):
     """Multifactorial of n of order k, n(!!...!).
 
     This is the multifactorial of n skipping k values.  For example,
@@ -3170,6 +3265,17 @@ def factorialk(n, k, exact=None):
         .. warning::
            The default value for ``exact`` will be changed to
            ``False`` in SciPy 1.15.0.
+    extend : string, optional
+        One of ``'zero'`` or ``'complex'``; this determines how values `n<0` are
+        handled - by default they are 0, but it is possible to opt into the complex
+        extension of the multifactorial. This enables passing complex values,
+        not only to ``n`` but also to ``k``.
+
+        .. warning::
+
+           Using the ``'complex'`` extension also changes the values of the
+           multifactorial at integers ``n != 1 (mod k)`` by a factor depending
+           on both ``k`` and ``n % k``, see below or [1].
 
     Returns
     -------
@@ -3197,15 +3303,23 @@ def factorialk(n, k, exact=None):
 
       n!(k) = k ** ((n - r)/k) * gamma(n/k + 1) / gamma(r/k + 1) * max(r, 1)
 
-    This is the basis of the approximation when ``exact=False``. Compare also [1].
+    This is the basis of the approximation when ``exact=False``.
+
+    In principle, any fixed choice of ``r`` (ignoring its relation ``r = n%k``
+    to ``n``) would provide a suitable analytic continuation from integer ``n``
+    to complex ``z`` (not only satisfying the functional equation but also
+    being logarithmically convex, c.f. Bohr-Mollerup theorem) -- in fact, the
+    choice of ``r`` above only changes the function by a constant factor. The
+    final constraint that determines the canonical continuation is ``f(1) = 1``,
+    which forces ``r = 1`` (see also [1]).::
+
+      z!(k) = k ** ((z - 1)/k) * gamma(z/k + 1) / gamma(1/k + 1)
 
     References
     ----------
     .. [1] Complex extension to multifactorial
             https://en.wikipedia.org/wiki/Double_factorial#Alternative_extension_of_the_multifactorial
     """
-    if not _is_subdtype(type(k), "i") or k < 1:
-        raise ValueError(f"k must be a positive integer, received: {k}")
     if exact is None:
         msg = (
             "factorialk will default to `exact=False` starting from SciPy "
@@ -3216,42 +3330,72 @@ def factorialk(n, k, exact=None):
         warnings.warn(msg, DeprecationWarning, stacklevel=2)
         exact = True
 
+    if extend not in ("zero", "complex"):
+        raise ValueError(
+            f"argument `extend` must be either 'zero' or 'complex', received: {extend}"
+        )
+    if exact and extend == "complex":
+        raise ValueError("Incompatible options: `exact=True` and `extend='complex'`")
+
+    msg_needs_complex = (
+        "In order to use non-integer arguments, you must opt into this by passing "
+        "`extend='complex'`. Note that this changes the result for all negative "
+        "arguments (which by default return 0). Additionally, it will perturb "
+        "the values of the multifactorial at most positive integers `n`."
+    )
     msg_wrong_dtype = (
-        "Unsupported data type for factorialk: {dtype}\n"
-        "Only integers are permitted."
+        "Unsupported data type for factorialk in {varname}: {dtype}\n"
+        "Only integers are permitted by default, though floating point "
+        "and complex numbers can be used if `extend='complex'` is passed."
     )
 
-    helpmsg = ""
-    if k in {1, 2}:
-        func = "factorial" if k == 1 else "factorial2"
-        helpmsg = f"\nYou can try to use {func} instead"
+    # check type of k
+    if not _is_subdtype(type(k), ["i", "f", "c"]):
+        raise ValueError(msg_wrong_dtype.format(varname="k", dtype=type(k)))
+    elif _is_subdtype(type(k), ["f", "c"]) and extend != "complex":
+        raise ValueError(msg_needs_complex)
+    # check value of k
+    if extend == "zero" and k < 1:
+        msg = f"For `extend='zero'`, k must be a positive integer, received: {k}"
+        raise ValueError(msg)
+    elif k == 0:
+        raise ValueError("Parameter k cannot be zero!")
 
     # don't use isscalar due to numpy/numpy#23574; 0-dim arrays treated below
     if np.ndim(n) == 0 and not isinstance(n, np.ndarray):
         # scalar cases
-        if not _is_subdtype(type(n), "i"):
-            raise ValueError(msg_wrong_dtype.format(dtype=type(n)) + helpmsg)
-        elif n < 0:
-            return 0
+        if not _is_subdtype(type(n), ["i", "f", "c", type(None)]):
+            raise ValueError(msg_wrong_dtype.format(varname="n", dtype=type(n)))
+        elif _is_subdtype(type(n), ["f", "c"]) and extend != "complex":
+            raise ValueError(msg_needs_complex)
+        elif n is None or np.isnan(n):
+            complexify = (extend == "complex") and _is_subdtype(type(n), "c")
+            return np.complex128("nan+nanj") if complexify else np.float64("nan")
+        elif extend == "zero" and n < 0:
+            return 0 if exact else np.float64(0)
         elif n in {0, 1}:
-            return 1
-        # general integer case
+            return 1 if exact else np.float64(1)
+
         if exact:
+            # general integer case
             return _range_prod(1, n, k=k)
-        return _factorialx_approx_core(n, k=k)
+        # approximation
+        return _factorialx_approx_core(n, k=k, extend=extend)
 
     # arrays & array-likes
     n = asarray(n)
 
-    if not _is_subdtype(n.dtype, "i"):
-        raise ValueError(msg_wrong_dtype.format(dtype=n.dtype) + helpmsg)
+    if not _is_subdtype(n.dtype, ["i", "f", "c"]):
+        raise ValueError(msg_wrong_dtype.format(varname="n", dtype=n.dtype))
+    elif _is_subdtype(n.dtype, ["f", "c"]) and extend != "complex":
+        raise ValueError(msg_needs_complex)
 
     if n.size == 0:
         # return empty arrays unchanged
         return n
     elif exact:
         return _factorialx_array_exact(n, k=k)
-    return _factorialx_array_approx(n, k=k)
+    return _factorialx_array_approx(n, k=k, extend=extend)
 
 
 def stirling2(N, K, *, exact=False):


Test patch:
diff --git a/scipy/_lib/tests/test_warnings.py b/scipy/_lib/tests/test_warnings.py
index f18bc20cb033..f200b1a6e975 100644
--- a/scipy/_lib/tests/test_warnings.py
+++ b/scipy/_lib/tests/test_warnings.py
@@ -115,6 +115,7 @@ def test_warning_calls_filters(warning_calls):
         os.path.join('optimize', '_nnls.py'),
         os.path.join('signal', '_ltisys.py'),
         os.path.join('sparse', '__init__.py'),  # np.matrix pending-deprecation
+        os.path.join('special', '_basic.py'),  # gh-21801
         os.path.join('stats', '_discrete_distns.py'),  # gh-14901
         os.path.join('stats', '_continuous_distns.py'),
         os.path.join('stats', '_binned_statistic.py'),  # gh-19345
diff --git a/scipy/special/tests/test_basic.py b/scipy/special/tests/test_basic.py
index 259a6fd35055..54d2bb0570f5 100644
--- a/scipy/special/tests/test_basic.py
+++ b/scipy/special/tests/test_basic.py
@@ -2133,28 +2133,39 @@ def assert_complex_nan(x):
             assert_really_equal(elem_x, elem_y, rtol=rtol)
     elif np.isnan(x) and np.isnan(y) and _is_subdtype(type(x), "c"):
         assert_complex_nan(x) and assert_complex_nan(y)
+    # no need to consider complex infinities due to numpy/numpy#25493
     else:
         assert_func(x, y)
 
 
 class TestFactorialFunctions:
-    @pytest.mark.parametrize("exact", [True, False])
-    def test_factorialx_scalar_return_type(self, exact):
-        kw = {"exact": exact}
+    def factorialk_ref(self, n, k, exact, extend):
+        if exact:
+            return special.factorialk(n, k=k, exact=True)
+        # for details / explanation see factorialk-docstring
+        r = np.mod(n, k) if extend == "zero" else 1
+        vals = np.power(k, (n - r)/k) * special.gamma(n/k + 1) * special.rgamma(r/k + 1)
+        # np.maximum is element-wise, which is what we want
+        return vals * np.maximum(r, 1)
+
+    @pytest.mark.parametrize("exact,extend",
+                             [(True, "zero"), (False, "zero"), (False, "complex")])
+    def test_factorialx_scalar_return_type(self, exact, extend):
+        kw = {"exact": exact, "extend": extend}
         assert np.isscalar(special.factorial(1, **kw))
         assert np.isscalar(special.factorial2(1, **kw))
         assert np.isscalar(special.factorialk(1, k=3, **kw))
 
     @pytest.mark.parametrize("n", [-1, -2, -3])
     @pytest.mark.parametrize("exact", [True, False])
-    def test_factorialx_negative(self, exact, n):
+    def test_factorialx_negative_extend_zero(self, exact, n):
         kw = {"exact": exact}
         assert_equal(special.factorial(n, **kw), 0)
         assert_equal(special.factorial2(n, **kw), 0)
         assert_equal(special.factorialk(n, k=3, **kw), 0)
 
     @pytest.mark.parametrize("exact", [True, False])
-    def test_factorialx_negative_array(self, exact):
+    def test_factorialx_negative_extend_zero_array(self, exact):
         kw = {"exact": exact}
         rtol = 1e-15
         n = [-5, -4, 0, 1]
@@ -2164,42 +2175,115 @@ def test_factorialx_negative_array(self, exact):
         assert_really_equal(special.factorial2(n, **kw), expected, rtol=rtol)
         assert_really_equal(special.factorialk(n, k=3, **kw), expected, rtol=rtol)
 
+    @pytest.mark.parametrize("n", [-1.1, -2.2, -3.3])
+    def test_factorialx_negative_extend_complex(self, n):
+        # factorialk defaults to exact=True, need to explicitly set to False
+        kw = {"exact": False, "extend": "complex"}
+        exp_1 = {-1.1: -10.686287021193184771,
+                 -2.2:   4.8509571405220931958,
+                 -3.3:  -1.4471073942559181166}
+        exp_2 = {-1.1:  1.0725776858167496309,
+                 -2.2: -3.9777171783768419874,
+                 -3.3: -0.99588841846200555977}
+        exp_k = {-1.1:  0.73565345382163025659,
+                 -2.2:  1.1749163167190809498,
+                 -3.3: -2.4780584257450583713}
+        rtol = 3e-15
+        assert_allclose(special.factorial(n, **kw), exp_1[n], rtol=rtol)
+        assert_allclose(special.factorial2(n, **kw), exp_2[n], rtol=rtol)
+        assert_allclose(special.factorialk(n, k=3, **kw), exp_k[n], rtol=rtol)
+        assert_allclose(special.factorial([n], **kw)[0], exp_1[n], rtol=rtol)
+        assert_allclose(special.factorial2([n], **kw)[0], exp_2[n], rtol=rtol)
+        assert_allclose(special.factorialk([n], k=3, **kw)[0], exp_k[n], rtol=rtol)
+
+    @pytest.mark.parametrize("imag", [0, 0j])
+    @pytest.mark.parametrize("n_outer", [-1, -2, -3])
+    def test_factorialx_negative_extend_complex_poles(self, n_outer, imag):
+        def _check(n):
+            complexify = _is_subdtype(type(n), "c")
+            # like for gamma, we expect complex nans for complex inputs
+            complex_nan = np.complex128("nan+nanj")
+            exp = np.complex128("nan+nanj") if complexify else np.float64("nan")
+            # factorialk defaults to incompatible exact=True, explicitly set to False
+            kw = {"exact": False, "extend": "complex"}
+            # poles are at negative integers that are multiples of k
+            assert_really_equal(special.factorial(n, **kw), exp)
+            assert_really_equal(special.factorial2(n * 2, **kw), exp)
+            assert_really_equal(special.factorialk(n * 3, k=3, **kw), exp)
+            # also test complex k for factorialk
+            c = 1.5 - 2j
+            assert_really_equal(special.factorialk(n * c, k=c, **kw), complex_nan)
+            # same for array case
+            assert_really_equal(special.factorial([n], **kw)[0], exp)
+            assert_really_equal(special.factorial2([n * 2], **kw)[0], exp)
+            assert_really_equal(special.factorialk([n * 3], k=3, **kw)[0], exp)
+            assert_really_equal(special.factorialk([n * c], k=c, **kw)[0], complex_nan)
+            # more specific tests in test_factorial{,2,k}_complex_reference
+
+        # imag ensures we test both real and complex representations of the poles
+        _check(n_outer + imag)
+        # check for large multiple of period
+        _check(100_000 * n_outer + imag)
+
     @pytest.mark.parametrize("boxed", [True, False])
+    @pytest.mark.parametrize("extend", ["zero", "complex"])
     @pytest.mark.parametrize(
         "n",
         [
             np.nan, np.float64("nan"), np.nan + np.nan*1j, np.complex128("nan+nanj"),
-            None, np.datetime64("nat")
+            np.inf, np.inf + 0j, -np.inf, -np.inf + 0j, None, np.datetime64("nat")
         ],
         ids=[
             "NaN", "np.float64('nan')", "NaN+i*NaN", "np.complex128('nan+nanj')",
-            "None", "NaT"
+            "inf", "inf+0i", "-inf", "-inf+0i", "None", "NaT"
         ]
     )
     @pytest.mark.parametrize(
         "factorialx",
         [special.factorial, special.factorial2, special.factorialk]
     )
-    def test_factorialx_nan(self, factorialx, n, boxed):
+    def test_factorialx_inf_nan(self, factorialx, n, extend, boxed):
         # NaNs not allowed (by dtype) for exact=True
-        kw = {"exact": False}
+        kw = {"exact": False, "extend": extend}
         if factorialx == special.factorialk:
             kw["k"] = 3
 
-        permissible_types = ["i"]
-        # factorial also allows floats
-        if factorialx == special.factorial:
-            # None is allowed for scalars, but would cause object type in array case
-            permissible_types = ["i", "f"] if boxed else ["i", "f", type(None)]
+        # None is allowed for scalars, but would cause object type in array case
+        permissible_types = ["i", "f", "c"] if boxed else ["i", "f", "c", type(None)]
+        # factorial allows floats also for extend="zero"
+        types_need_complex_ext = "c" if factorialx == special.factorial else ["f", "c"]
 
         if not _is_subdtype(type(n), permissible_types):
             with pytest.raises(ValueError, match="Unsupported data type.*"):
                 factorialx([n] if boxed else n, **kw)
+        elif _is_subdtype(type(n), types_need_complex_ext) and extend != "complex":
+            with pytest.raises(ValueError, match="In order to use non-integer.*"):
+                factorialx([n] if boxed else n, **kw)
         else:
+            # account for type and whether extend="complex"
+            complexify = (extend == "complex") and _is_subdtype(type(n), "c")
+            # note that the type of the naïve `np.nan + np.nan * 1j` is `complex`
+            # instead of `numpy.complex128`, which trips up assert_really_equal
+            expected = np.complex128("nan+nanj") if complexify else np.float64("nan")
+            # the only exception are real infinities
+            if _is_subdtype(type(n), "f") and np.isinf(n):
+                # unchanged for positive infinity; negative one depends on extension
+                neg_inf_result = np.float64(0 if (extend == "zero") else "nan")
+                expected = np.float64("inf") if (n > 0) else neg_inf_result
+
             result = factorialx([n], **kw)[0] if boxed else factorialx(n, **kw)
-            assert_really_equal(result, np.float64("nan"))
+            assert_really_equal(result, expected)
             # also tested in test_factorial{,2,k}_{array,scalar}_corner_cases
 
+    @pytest.mark.parametrize("extend", [0, 1.1, np.nan, "string"])
+    def test_factorialx_raises_extend(self, extend):
+        with pytest.raises(ValueError, match="argument `extend` must be.*"):
+            special.factorial(1, extend=extend)
+        with pytest.raises(ValueError, match="argument `extend` must be.*"):
+            special.factorial2(1, extend=extend)
+        with pytest.raises(ValueError, match="argument `extend` must be.*"):
+            special.factorialk(1, k=3, exact=True, extend=extend)
+
     @pytest.mark.parametrize("levels", range(1, 5))
     @pytest.mark.parametrize("exact", [True, False])
     def test_factorialx_array_shape(self, levels, exact):
@@ -2294,15 +2378,20 @@ def test_factorial_int_reference(self, n):
         assert_allclose(correct, special.factorial(n, exact=False), rtol=rtol)
         assert_allclose(correct, special.factorial([n], exact=False)[0], rtol=rtol)
 
+        # extend="complex" only works for exact=False
+        kw = {"exact": False, "extend": "complex"}
+        assert_allclose(correct, special.factorial(n, **kw), rtol=rtol)
+        assert_allclose(correct, special.factorial([n], **kw)[0], rtol=rtol)
+
     def test_factorial_float_reference(self):
         def _check(n, expected):
             rtol = 8e-14 if sys.platform == 'win32' else 1e-15
             assert_allclose(special.factorial(n), expected, rtol=rtol)
             assert_allclose(special.factorial([n])[0], expected, rtol=rtol)
             # using floats with `exact=True` raises an error for scalars and arrays
-            with pytest.raises(ValueError, match="Non-integer values.*"):
-                assert_allclose(special.factorial(n, exact=True), expected, rtol=rtol)
-            with pytest.raises(ValueError, match="factorial with `exact=Tr.*"):
+            with pytest.raises(ValueError, match="`exact=True` only supports.*"):
+                special.factorial(n, exact=True)
+            with pytest.raises(ValueError, match="`exact=True` only supports.*"):
                 special.factorial([n], exact=True)
 
         # Reference values from mpmath for gamma(n+1)
@@ -2317,8 +2406,25 @@ def _check(n, expected):
         # close to maximum for float64
         _check(170.6243, 1.79698185749571048960082e+308)
 
+    def test_factorial_complex_reference(self):
+        def _check(n, expected):
+            rtol = 3e-15 if sys.platform == 'win32' else 2e-15
+            kw = {"exact": False, "extend": "complex"}
+            assert_allclose(special.factorial(n, **kw), expected, rtol=rtol)
+            assert_allclose(special.factorial([n], **kw)[0], expected, rtol=rtol)
+
+        # Reference values from mpmath.gamma(n+1)
+        # negative & complex values
+        _check(-0.5, expected=1.7724538509055160276)
+        _check(-0.5 + 0j, expected=1.7724538509055160276 + 0j)
+        _check(2 + 2j, expected=-0.42263728631120216694 + 0.87181425569650686062j)
+        # close to poles
+        _check(-0.9999, expected=9999.422883232725532)
+        _check(-1 + 0.0001j, expected=-0.57721565582674219 - 9999.9999010944009697j)
+
     @pytest.mark.parametrize("dtype", [np.int64, np.float64,
                                        np.complex128, object])
+    @pytest.mark.parametrize("extend", ["zero", "complex"])
     @pytest.mark.parametrize("exact", [True, False])
     @pytest.mark.parametrize("dim", range(0, 5))
     # test empty & non-empty arrays, with nans and mixed
@@ -2327,7 +2433,7 @@ def _check(n, expected):
         [[], [1], [1.1], [np.nan], [np.nan + np.nan * 1j], [np.nan, 1]],
         ids=["[]", "[1]", "[1.1]", "[NaN]", "[NaN+i*NaN]", "[NaN, 1]"],
     )
-    def test_factorial_array_corner_cases(self, content, dim, exact, dtype):
+    def test_factorial_array_corner_cases(self, content, dim, exact, extend, dtype):
         if dtype is object and SCIPY_ARRAY_API:
             pytest.skip("object arrays unsupported in array API mode")
         # get dtype without calling array constructor (that might fail or mutate)
@@ -2336,17 +2442,23 @@ def test_factorial_array_corner_cases(self, content, dim, exact, dtype):
         if dtype == np.float64 and any(_is_subdtype(type(x), "c") for x in content):
             pytest.skip("impossible combination")
 
-        kw = {"exact": exact}
+        kw = {"exact": exact, "extend": extend}
         # np.array(x, ndim=0) will not be 0-dim. unless x is too
         content = content if (dim > 0 or len(content) != 1) else content[0]
         n = np.array(content, ndmin=dim, dtype=dtype)
 
         result = None
-        if not _is_subdtype(n.dtype, ["i", "f"]):
+        if extend == "complex" and exact:
+            with pytest.raises(ValueError, match="Incompatible options:.*"):
+                special.factorial(n, **kw)
+        elif not _is_subdtype(n.dtype, ["i", "f", "c"]):
             with pytest.raises(ValueError, match="Unsupported data type.*"):
                 special.factorial(n, **kw)
+        elif _is_subdtype(n.dtype, "c") and extend != "complex":
+            with pytest.raises(ValueError, match="In order to use non-integer.*"):
+                special.factorial(n, **kw)
         elif exact and not _is_subdtype(n.dtype, "i"):
-            with pytest.raises(ValueError, match="factorial with `exact=.*"):
+            with pytest.raises(ValueError, match="`exact=True` only supports.*"):
                 special.factorial(n, **kw)
         else:
             result = special.factorial(n, **kw)
@@ -2359,22 +2471,33 @@ def test_factorial_array_corner_cases(self, content, dim, exact, dtype):
             # result is empty if and only if n is empty, and has the same dimension
             # as n; dtype stays the same, except when not empty and not exact:
             if n.size:
-                dtype = native_int if exact else np.float64
+                cx = (extend == "complex") and _is_subdtype(n.dtype, "c")
+                dtype = np.complex128 if cx else (native_int if exact else np.float64)
             expected = np.array(ref, ndmin=dim, dtype=dtype)
-            assert_really_equal(result, expected)
+            assert_really_equal(result, expected, rtol=1e-15)
 
+    @pytest.mark.parametrize("extend", ["zero", "complex"])
     @pytest.mark.parametrize("exact", [True, False])
     @pytest.mark.parametrize("n", [1, 1.1, 2 + 2j, np.nan, np.nan + np.nan*1j, None],
                              ids=["1", "1.1", "2+2j", "NaN", "NaN+i*NaN", "None"])
-    def test_factorial_scalar_corner_cases(self, n, exact):
-        kw = {"exact": exact}
-        if not _is_subdtype(type(n), ["i", "f", type(None)]):
+    def test_factorial_scalar_corner_cases(self, n, exact, extend):
+        kw = {"exact": exact, "extend": extend}
+        if extend == "complex" and exact:
+            with pytest.raises(ValueError, match="Incompatible options:.*"):
+                special.factorial(n, **kw)
+        elif not _is_subdtype(type(n), ["i", "f", "c", type(None)]):
             with pytest.raises(ValueError, match="Unsupported data type.*"):
                 special.factorial(n, **kw)
+        elif _is_subdtype(type(n), "c") and extend != "complex":
+            with pytest.raises(ValueError, match="In order to use non-integer.*"):
+                special.factorial(n, **kw)
         elif n is None or np.isnan(n):
-            assert_really_equal(special.factorial(n, **kw), np.float64("nan"))
+            # account for dtype and whether extend="complex"
+            complexify = (extend == "complex") and _is_subdtype(type(n), "c")
+            expected = np.complex128("nan+nanj") if complexify else np.float64("nan")
+            assert_really_equal(special.factorial(n, **kw), expected)
         elif exact and _is_subdtype(type(n), "f"):
-            with pytest.raises(ValueError, match="Non-integer values.*"):
+            with pytest.raises(ValueError, match="`exact=True` only supports.*"):
                 special.factorial(n, **kw)
         else:
             assert_equal(special.factorial(n, **kw), special.gamma(n + 1))
@@ -2409,8 +2532,37 @@ def test_factorial2_int_reference(self, n):
         assert_allclose(correct, special.factorial2(n, exact=False), rtol=rtol)
         assert_allclose(correct, special.factorial2([n], exact=False)[0], rtol=rtol)
 
+        # extend="complex" only works for exact=False
+        kw = {"exact": False, "extend": "complex"}
+        # approximation only matches exactly for `n == 1 (mod k)`, see docstring
+        if n % 2 == 1:
+            assert_allclose(correct, special.factorial2(n, **kw), rtol=rtol)
+            assert_allclose(correct, special.factorial2([n], **kw)[0], rtol=rtol)
+
+    def test_factorial2_complex_reference(self):
+        # this tests for both floats and complex
+        def _check(n, expected):
+            rtol = 5e-15
+            kw = {"exact": False, "extend": "complex"}
+            assert_allclose(special.factorial2(n, **kw), expected, rtol=rtol)
+            assert_allclose(special.factorial2([n], **kw)[0], expected, rtol=rtol)
+
+        # Reference values from mpmath for:
+        # mpmath.power(2, n/2) * mpmath.gamma(n/2 + 1) * mpmath.sqrt(2 / mpmath.pi)
+        _check(3, expected=3)
+        _check(4, expected=special.factorial2(4) * math.sqrt(2 / math.pi))
+        _check(20, expected=special.factorial2(20) * math.sqrt(2 / math.pi))
+        # negative & complex values
+        _check(-0.5, expected=0.82217895866245855122)
+        _check(-0.5 + 0j, expected=0.82217895866245855122 + 0j)
+        _check(3 + 3j, expected=-1.0742236630142471526 + 1.4421398439387262897j)
+        # close to poles
+        _check(-1.9999, expected=7978.8918745523440682)
+        _check(-2 + 0.0001j, expected=0.0462499835314308444 - 7978.84559148876374493j)
+
     @pytest.mark.parametrize("dtype", [np.int64, np.float64,
                                        np.complex128, object])
+    @pytest.mark.parametrize("extend", ["zero", "complex"])
     @pytest.mark.parametrize("exact", [True, False])
     @pytest.mark.parametrize("dim", range(0, 5))
     # test empty & non-empty arrays, with nans and mixed
@@ -2419,22 +2571,28 @@ def test_factorial2_int_reference(self, n):
         [[], [1], [1.1], [np.nan], [np.nan + np.nan * 1j], [np.nan, 1]],
         ids=["[]", "[1]", "[1.1]", "[NaN]", "[NaN+i*NaN]", "[NaN, 1]"],
     )
-    def test_factorial2_array_corner_cases(self, content, dim, exact, dtype):
+    def test_factorial2_array_corner_cases(self, content, dim, exact, extend, dtype):
         # get dtype without calling array constructor (that might fail or mutate)
         if dtype == np.int64 and any(np.isnan(x) or (x != int(x)) for x in content):
             pytest.skip("impossible combination")
         if dtype == np.float64 and any(_is_subdtype(type(x), "c") for x in content):
             pytest.skip("impossible combination")
 
-        kw = {"exact": exact}
+        kw = {"exact": exact, "extend": extend}
         # np.array(x, ndim=0) will not be 0-dim. unless x is too
         content = content if (dim > 0 or len(content) != 1) else content[0]
         n = np.array(content, ndmin=dim, dtype=dtype)
 
         result = None
-        if not _is_subdtype(n.dtype, "i"):
+        if extend == "complex" and exact:
+            with pytest.raises(ValueError, match="Incompatible options:.*"):
+                special.factorial2(n, **kw)
+        elif not _is_subdtype(n.dtype, ["i", "f", "c"]):
             with pytest.raises(ValueError, match="Unsupported data type.*"):
                 special.factorial2(n, **kw)
+        elif _is_subdtype(n.dtype, ["f", "c"]) and extend != "complex":
+            with pytest.raises(ValueError, match="In order to use non-integer.*"):
+                special.factorial2(n, **kw)
         else:
             result = special.factorial2(n, **kw)
 
@@ -2446,20 +2604,34 @@ def test_factorial2_array_corner_cases(self, content, dim, exact, dtype):
             # result is empty if and only if n is empty, and has the same dimension
             # as n; dtype stays the same, except when not empty and not exact:
             if n.size:
-                dtype = native_int if exact else np.float64
+                cx = (extend == "complex") and _is_subdtype(n.dtype, "c")
+                dtype = np.complex128 if cx else (native_int if exact else np.float64)
             expected = np.array(ref, ndmin=dim, dtype=dtype)
-            assert_really_equal(result, expected, rtol=1e-15)
+            assert_really_equal(result, expected, rtol=2e-15)
 
+    @pytest.mark.parametrize("extend", ["zero", "complex"])
     @pytest.mark.parametrize("exact", [True, False])
     @pytest.mark.parametrize("n", [1, 1.1, 2 + 2j, np.nan, np.nan + np.nan*1j, None],
                              ids=["1", "1.1", "2+2j", "NaN", "NaN+i*NaN", "None"])
-    def test_factorial2_scalar_corner_cases(self, n, exact):
-        kw = {"exact": exact}
-        if not _is_subdtype(type(n), "i"):
+    def test_factorial2_scalar_corner_cases(self, n, exact, extend):
+        kw = {"exact": exact, "extend": extend}
+        if extend == "complex" and exact:
+            with pytest.raises(ValueError, match="Incompatible options:.*"):
+                special.factorial2(n, **kw)
+        elif not _is_subdtype(type(n), ["i", "f", "c", type(None)]):
             with pytest.raises(ValueError, match="Unsupported data type.*"):
                 special.factorial2(n, **kw)
+        elif _is_subdtype(type(n), ["f", "c"]) and extend != "complex":
+            with pytest.raises(ValueError, match="In order to use non-integer.*"):
+                special.factorial2(n, **kw)
+        elif n is None or np.isnan(n):
+            # account for dtype and whether extend="complex"
+            complexify = (extend == "complex") and _is_subdtype(type(n), "c")
+            expected = np.complex128("nan+nanj") if complexify else np.float64("nan")
+            assert_really_equal(special.factorial2(n, **kw), expected)
         else:
-            assert_equal(special.factorial2(n, **kw), 1)
+            expected = self.factorialk_ref(n, k=2, **kw)
+            assert_really_equal(special.factorial2(n, **kw), expected, rtol=1e-15)
 
     @pytest.mark.parametrize("k", range(1, 5))
     # note that n=170 is the last integer such that factorial(n) fits float64;
@@ -2495,8 +2667,47 @@ def test_factorialk_int_reference(self, n, k):
         assert_allclose(correct, special.factorialk(n, k, exact=False), rtol=rtol)
         assert_allclose(correct, special.factorialk([n], k, exact=False)[0], rtol=rtol)
 
+        # extend="complex" only works for exact=False
+        kw = {"k": k, "exact": False, "extend": "complex"}
+        # approximation only matches exactly for `n == 1 (mod k)`, see docstring
+        if n % k == 1:
+            rtol = 2e-14
+            assert_allclose(correct, special.factorialk(n, **kw), rtol=rtol)
+            assert_allclose(correct, special.factorialk([n], **kw)[0], rtol=rtol)
+
+    def test_factorialk_complex_reference(self):
+        # this tests for both floats and complex
+        def _check(n, k, exp):
+            rtol = 1e-14
+            kw = {"k": k, "exact": False, "extend": "complex"}
+            assert_allclose(special.factorialk(n, **kw), exp, rtol=rtol)
+            assert_allclose(special.factorialk([n], **kw)[0], exp, rtol=rtol)
+
+        # Reference values from mpmath for:
+        # mpmath.power(k, (n-1)/k) * mpmath.gamma(n/k + 1) / mpmath.gamma(1/k + 1)
+        _check(n=4, k=3, exp=special.factorialk(4, k=3, exact=True))
+        _check(n=5, k=3, exp=7.29011132947227083)
+        _check(n=6.5, k=3, exp=19.6805080113566010)
+        # non-integer k
+        _check(n=3, k=2.5, exp=2.58465740293218541)
+        _check(n=11, k=2.5, exp=1963.5)  # ==11*8.5*6*3.5; c.f. n == 1 (mod k)
+        _check(n=-3 + 3j + 1, k=-3 + 3j, exp=-2 + 3j)
+        # complex values
+        _check(n=4 + 4j, k=4, exp=-0.67855904082768043854 + 2.1993925819930311497j)
+        _check(n=4, k=4 - 4j, exp=1.9775338957222718742 + 0.92607172675423901371j)
+        _check(n=4 + 4j, k=4 - 4j, exp=0.1868492880824934475 + 0.87660580316894290247j)
+        # negative values
+        _check(n=-0.5, k=3, exp=0.72981013240713739354)
+        _check(n=-0.5 + 0j, k=3, exp=0.72981013240713739354 + 0j)
+        _check(n=2.9, k=-0.7, exp=0.45396591474966867296 + 0.56925525174685228866j)
+        _check(n=-0.6, k=-0.7, exp=-0.07190820089634757334 - 0.090170031876701730081j)
+        # close to poles
+        _check(n=-2.9999, k=3, exp=7764.7170695908828364)
+        _check(n=-3 + 0.0001j, k=3, exp=0.1349475632879599864 - 7764.5821055158365027j)
+
     @pytest.mark.parametrize("dtype", [np.int64, np.float64,
                                        np.complex128, object])
+    @pytest.mark.parametrize("extend", ["zero", "complex"])
     @pytest.mark.parametrize("exact", [True, False])
     @pytest.mark.parametrize("dim", range(0, 5))
     # test empty & non-empty arrays, with nans and mixed
@@ -2505,22 +2716,28 @@ def test_factorialk_int_reference(self, n, k):
         [[], [1], [1.1], [np.nan], [np.nan + np.nan * 1j], [np.nan, 1]],
         ids=["[]", "[1]", "[1.1]", "[NaN]", "[NaN+i*NaN]", "[NaN, 1]"],
     )
-    def test_factorialk_array_corner_cases(self, content, dim, exact, dtype):
+    def test_factorialk_array_corner_cases(self, content, dim, exact, extend, dtype):
         # get dtype without calling array constructor (that might fail or mutate)
         if dtype == np.int64 and any(np.isnan(x) or (x != int(x)) for x in content):
             pytest.skip("impossible combination")
         if dtype == np.float64 and any(_is_subdtype(type(x), "c") for x in content):
             pytest.skip("impossible combination")
 
-        kw = {"k": 3, "exact": exact}
+        kw = {"k": 3, "exact": exact, "extend": extend}
         # np.array(x, ndim=0) will not be 0-dim. unless x is too
         content = content if (dim > 0 or len(content) != 1) else content[0]
         n = np.array(content, ndmin=dim, dtype=dtype)
 
         result = None
-        if not _is_subdtype(n.dtype, "i"):
+        if extend == "complex" and exact:
+            with pytest.raises(ValueError, match="Incompatible options:.*"):
+                special.factorialk(n, **kw)
+        elif not _is_subdtype(n.dtype, ["i", "f", "c"]):
             with pytest.raises(ValueError, match="Unsupported data type.*"):
                 special.factorialk(n, **kw)
+        elif _is_subdtype(n.dtype, ["f", "c"]) and extend != "complex":
+            with pytest.raises(ValueError, match="In order to use non-integer.*"):
+                special.factorialk(n, **kw)
         else:
             result = special.factorialk(n, **kw)
 
@@ -2532,22 +2749,35 @@ def test_factorialk_array_corner_cases(self, content, dim, exact, dtype):
             # result is empty if and only if n is empty, and has the same dimension
             # as n; dtype stays the same, except when not empty and not exact:
             if n.size:
-                dtype = native_int if exact else np.float64
+                cx = (extend == "complex") and _is_subdtype(n.dtype, "c")
+                dtype = np.complex128 if cx else (native_int if exact else np.float64)
             expected = np.array(ref, ndmin=dim, dtype=dtype)
-            assert_really_equal(result, expected, rtol=1e-15)
+            assert_really_equal(result, expected, rtol=2e-15)
 
+    @pytest.mark.parametrize("extend", ["zero", "complex"])
     @pytest.mark.parametrize("exact", [True, False])
     @pytest.mark.parametrize("k", range(1, 5))
     @pytest.mark.parametrize("n", [1, 1.1, 2 + 2j, np.nan, np.nan + np.nan*1j, None],
                              ids=["1", "1.1", "2+2j", "NaN", "NaN+i*NaN", "None"])
-    def test_factorialk_scalar_corner_cases(self, n, k, exact):
-        kw = {"k": k, "exact": exact}
-        if not _is_subdtype(type(n), "i"):
+    def test_factorialk_scalar_corner_cases(self, n, k, exact, extend):
+        kw = {"k": k, "exact": exact, "extend": extend}
+        if extend == "complex" and exact:
+            with pytest.raises(ValueError, match="Incompatible options:.*"):
+                special.factorialk(n, **kw)
+        elif not _is_subdtype(type(n), ["i", "f", "c", type(None)]):
             with pytest.raises(ValueError, match="Unsupported data type.*"):
                 special.factorialk(n, **kw)
+        elif _is_subdtype(type(n), ["f", "c"]) and extend != "complex":
+            with pytest.raises(ValueError, match="In order to use non-integer.*"):
+                special.factorialk(n, **kw)
+        elif n is None or np.isnan(n):
+            # account for dtype and whether extend="complex"
+            complexify = (extend == "complex") and _is_subdtype(type(n), "c")
+            expected = np.complex128("nan+nanj") if complexify else np.float64("nan")
+            assert_really_equal(special.factorialk(n, **kw), expected)
         else:
-            # factorialk(1, k) == 1 for all k
-            assert_equal(special.factorialk(n, **kw), 1)
+            expected = self.factorialk_ref(n, **kw)
+            assert_really_equal(special.factorialk(n, **kw), expected, rtol=1e-15)
 
     @pytest.mark.parametrize("k", range(1, 5))
     def test_factorialk_deprecation_exact(self, k):
@@ -2557,29 +2787,42 @@ def test_factorialk_deprecation_exact(self, k):
             special.factorialk(1, k=k)
 
     @pytest.mark.parametrize("boxed", [True, False])
-    @pytest.mark.parametrize("exact", [True, False])
-    @pytest.mark.parametrize("k", [0, 1.1, np.nan])
-    def test_factorialk_raises_k_complex(self, k, exact, boxed):
+    @pytest.mark.parametrize("exact,extend",
+                             [(True, "zero"), (False, "zero"), (False, "complex")])
+    @pytest.mark.parametrize("k", [-1, -1.0, 0, 0.0, 0 + 1j, 1.1, np.nan])
+    def test_factorialk_raises_k_complex(self, k, exact, extend, boxed):
         n = [1] if boxed else 1
-        kw = {"k": k, "exact": exact}
-        with pytest.raises(ValueError, match="k must be a positive integer*"):
+        kw = {"k": k, "exact": exact, "extend": extend}
+        if extend == "zero":
+            msg = "In order to use non-integer.*"
+            if _is_subdtype(type(k), "i") and (k < 1):
+                msg = "For `extend='zero'`.*"
+            with pytest.raises(ValueError, match=msg):
+                special.factorialk(n, **kw)
+        elif k == 0:
+            with pytest.raises(ValueError, match="Parameter k cannot be zero!"):
+                special.factorialk(n, **kw)
+        else:
+            # no error
             special.factorialk(n, **kw)
 
     @pytest.mark.parametrize("boxed", [True, False])
-    @pytest.mark.parametrize("exact", [True, False])
+    @pytest.mark.parametrize("exact,extend",
+                             [(True, "zero"), (False, "zero"), (False, "complex")])
     # neither integer, float nor complex
     @pytest.mark.parametrize("k", ["string", np.datetime64("nat")],
                              ids=["string", "NaT"])
-    def test_factorialk_raises_k_other(self, k, exact, boxed):
+    def test_factorialk_raises_k_other(self, k, exact, extend, boxed):
         n = [1] if boxed else 1
-        kw = {"k": k, "exact": exact}
-        with pytest.raises(ValueError, match="k must be a positive integer*"):
+        kw = {"k": k, "exact": exact, "extend": extend}
+        with pytest.raises(ValueError, match="Unsupported data type.*"):
             special.factorialk(n, **kw)
 
-    @pytest.mark.parametrize("exact", [True, False])
+    @pytest.mark.parametrize("exact,extend",
+                             [(True, "zero"), (False, "zero"), (False, "complex")])
     @pytest.mark.parametrize("k", range(1, 12))
-    def test_factorialk_dtype(self, k, exact):
-        kw = {"k": k, "exact": exact}
+    def test_factorialk_dtype(self, k, exact, extend):
+        kw = {"k": k, "exact": exact, "extend": extend}
         if exact and k in _FACTORIALK_LIMITS_64BITS.keys():
             n = np.array([_FACTORIALK_LIMITS_32BITS[k]])
             assert_equal(special.factorialk(n, **kw).dtype, np_long)
@@ -2594,7 +2837,7 @@ def test_factorialk_dtype(self, k, exact):
         else:
             n = np.array([_FACTORIALK_LIMITS_64BITS.get(k, 1)])
             # for exact=True and k >= 10, we always return object;
-            # for exact=False it's always float
+            # for exact=False it's always float (unless input is complex)
             dtype = object if exact else np.float64
             assert_equal(special.factorialk(n, **kw).dtype, dtype)
 
@@ -2602,7 +2845,7 @@ def test_factorial_mixed_nan_inputs(self):
         x = np.array([np.nan, 1, 2, 3, np.nan])
         expected = np.array([np.nan, 1, 2, 6, np.nan])
         assert_equal(special.factorial(x, exact=False), expected)
-        with pytest.raises(ValueError, match="factorial with `exact=True.*"):
+        with pytest.raises(ValueError, match="`exact=True` only supports.*"):
             special.factorial(x, exact=True)