Repository: scipy/scipy
Issue category: Feature Request

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



Functions to edit:
scipy/special/_basic.py:_factorialx_array_approx
scipy/special/_basic.py:_factorialx_approx_core
scipy/special/_basic.py:factorial
scipy/special/_basic.py:factorial2
scipy/special/_basic.py:factorialk
Functions added: scipy/special/_basic.py:_gamma1p

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)