BUG: fix incorrect values in factorial for 0 with uint dtype (scipy#22168)

Author: h-vetinari

scipy/special/_basic.py

@@ -2972,8 +2972,9 @@ def _factorialx_approx_core(n, k, extend):
         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)
+            # don't use `(n-1)/k` in np.power; underflows if 0 is of a uintX type
+            result = np.power(k, n / k, dtype=p_dtype) * _gamma1p(n / k)
+            result *= rgamma(1 / k + 1) / np.power(k, 1 / k, dtype=p_dtype)
         if isinstance(n, np.ndarray):
             # ensure we keep array-ness for 0-dim inputs; already n/k above loses it
             result = np.array(result)

scipy/special/tests/test_basic.py

@@ -2109,11 +2109,12 @@ def assert_really_equal(x, y, rtol=None):
     Sharper assertion function that is stricter about matching types, not just values
 
     This is useful/necessary in some cases:
-      * handled by xp_assert_* functions
       * dtypes for arrays that have the same _values_ (e.g. element 1.0 vs 1)
       * distinguishing complex from real NaN
+      * result types for scalars
 
     We still want to be able to allow a relative tolerance for the values though.
+    The main logic comparison logic is handled by the xp_assert_* functions.
     """
     def assert_func(x, y):
         xp_assert_equal(x, y) if rtol is None else xp_assert_close(x, y, rtol=rtol)
@@ -2350,6 +2351,24 @@ def _nest_me(x, k=1):
         assert_func(special.factorialk(n, 3, exact=exact),
                     np.array(exp_nucleus[3], ndmin=level))
 
+    @pytest.mark.parametrize("dtype", [np.uint8, np.uint16, np.uint32, np.uint64])
+    @pytest.mark.parametrize("exact,extend",
+                             [(True, "zero"), (False, "zero"), (False, "complex")])
+    def test_factorialx_uint(self, exact, extend, dtype):
+        # ensure that uint types work correctly as inputs
+        kw = {"exact": exact, "extend": extend}
+        assert_func = assert_array_equal if exact else assert_allclose
+        def _check(n):
+            n_ref = n.astype(np.int64) if isinstance(n, np.ndarray) else np.int64(n)
+            assert_func(special.factorial(n, **kw), special.factorial(n_ref, **kw))
+            assert_func(special.factorial2(n, **kw), special.factorial2(n_ref, **kw))
+            assert_func(special.factorialk(n, k=3, **kw),
+                        special.factorialk(n_ref, k=3, **kw))
+        _check(dtype(0))
+        _check(dtype(1))
+        _check(np.array(0, dtype=dtype))
+        _check(np.array([0, 1], dtype=dtype))
+
     # note that n=170 is the last integer such that factorial(n) fits float64
     @pytest.mark.parametrize('n', range(30, 180, 10))
     def test_factorial_accuracy(self, n):