Repository: scipy/scipy
Pull request: ENH: stats.mannwhitneyu: replace exact p-value calculation
#### Reference issue
gh-4933
Closes gh-17559

#### What does this implement/fix?
https://github.com/scipy/scipy/pull/4933#issuecomment-1898082691 suggested a different algorithm for the `mannwhitneyu` exact p-value calculation that greatly improves performance (e.g. factor of 675 for $m \approx n \approx 100$) and reduces memory consumption from $O(m^2n^2)$ to $O(mn)$. This PR implements the suggestion.

#### Additional information
@AtsushiSakai let's wait for confirmation from @toobaz on use of the code before starting a detailed review. In the meantime, I thought you'd be interested in this:

```python3
import numpy as np
from scipy import stats
from scipy.stats._mannwhitneyu import _mwu_state
rng = np.random.default_rng(1305672164505329978)
x = rng.random(100)
y = rng.random(101)
```

Before:
```python3
res = stats.mannwhitneyu(x, y, method='exact')  # 16.9 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)
res  # MannwhitneyuResult(statistic=4792.0, pvalue=0.5331021234219595)
_mwu_state._fmnks.shape  #  (101, 102, 4793)
```

After:
```python3
res = stats.mannwhitneyu(x, y, method='exact')  # 25 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)
res  # MannwhitneyuResult(statistic=4792.0, pvalue=0.5331021234219553)
_mwu_state.configurations.shape  #  (4793,)
```



Problem statement:
ENH: _mannwhitneyu.py computation of exact MWU statistics may be simplified
### Describe your issue.

This is not really a bug but an invitation for improvement and simplification. Imho the _MWU class in _mannwhitneyu.py may be abolished.
In pmf_recursive() or pmf_iterative() a 3-dimensional array is used, which stores fmnk s for all combinations of m, n, k i.e. stores k-distributions for all values of m and n.
This is not necessary as formula (5) of DiBucchianico:Combinatorics_computer_algebra_and_wilcoxon_mann_whitney_test
allows for recursive updating from all f(m)nk to f(m+1)nk by adding a k-shifted f(m)(n-1)k to f(m)nk.
An example: f2,0,k = (1,0,..) f2,1,k = (1,1,1,0,..) f2,2,k=(1,1,2,1,1,0,..),..
Now f3,1,k=(1,1,1,1,0,..) = f2,1,k + rightshift(f3,0,k|3) and f3,2,k=(1,1,2,2,2,1,1,0,..) = f2,2,k + rightshift( f3,1,k |3).
This is faster and requires only a 2-dimensional array. This way fmnk may be computed for higher m and n values.
Code for this is below, free for your adaptation. 

### Reproducing Code Example

```python
import numpy as np
from scipy import special

# Computes f_mnk distribution needed for exact Mann-Whitney U statistics using formula (5) of
# DiBucchianico (1996): "Combinatorics, computer algebra and Wilcoxon-Mann-Whitney test".
#
# We use unsigned quads with max value of 2**64-1 = 18446744073709551615
#    as long as a safety test to check if math_comb(m+n,n) > 2**64-1 is fulfilled
#
def fmnk(m, n, k=None):  # if k is None, return full dist: (fmn0, fmn1,.., fmn(m*n)), else only (fmn0, fmn1,.., fmnk)
    dtyp = np.longdouble if special.binom(n+m, n) > 2**64 else np.uint64  # overflow protection: sum of all densities fits -> every entry fits
    n, m = sorted((m, n))  # make n the smaller of both (n requires memory), we iterate through m
    if k is None:          # v-- original: l = m*n+1   # return array with full fmn distribution for k  = 0,..,m*n
        l =  (m*n+1)//2 if m*n % 2 else (m*n)//2 + 1   # we compute only half of all distributions (including center), for symmetry reasons
    else:
        l = k+1            # compute and return only k-densities for k = 0,..,k
    fm = np.zeros((n+1, l), dtype=dtyp)
    fm[:,0] = 1
    for j in range(1, m+1):       # update all fjik until fmik in place: j=1,..., m
        for i in range(1, n+1):   # fj0k remains (1,0,..0); fjik will be updated for i >= 1
            fm[i,j:]+=fm[i-1,:-j] # add right shifted (by m value which is j) previous (n-1-th) k-distributrion
    if k is not None:             # this is formula (5) of DiBucchianico:Combinatorics_computer_algebra_and_wilcoxon_mann_whitney_test
        return fm[n,:]     # return f_mn which is k-density for k = 0,..,k
    r = np.zeros((m*n+1), dtype=dtyp)  # k was None and we had only computed half of the k distribution
    r[:l] = fm[n,:l]
    print(f'')
    r[l:] = fm[n,l-1::-1] if m*n % 2 else fm[n,l-2::-1]  # so we return a new array and fill it symmetrically                 
    return r          # return full f_mn which is k-density for k = 0,..,m*n

print(fmnk(3,3))
```


### Error message

```shell
[1 1 2 3 3 3 3 2 1 1]
```


### SciPy/NumPy/Python version information

1.9.3 1.23.4 sys.version_info(major=3, minor=10, micro=8, releaselevel='final', serial=0)

Hints from issue discussion:

Commit message: MAINT: stats.mannwhitneyu: fix CI test failure

Patch to review:
diff --git a/scipy/stats/_mannwhitneyu.py b/scipy/stats/_mannwhitneyu.py
index a796edf075e7..77c7b27f8bc5 100644
--- a/scipy/stats/_mannwhitneyu.py
+++ b/scipy/stats/_mannwhitneyu.py
@@ -19,14 +19,29 @@ def _broadcast_concatenate(x, y, axis):
 
 class _MWU:
     '''Distribution of MWU statistic under the null hypothesis'''
-    # Possible improvement: if m and n are small enough, use integer arithmetic
 
-    def __init__(self):
-        '''Minimal initializer'''
-        self._fmnks = -np.ones((1, 1, 1))
-        self._recursive = None
+    def __init__(self, n1, n2):
+        self._reset(n1, n2)
 
-    def pmf(self, k, m, n):
+    def set_shapes(self, n1, n2):
+        n1, n2 = min(n1, n2), max(n1, n2)
+        if (n1, n2) == (self.n1, self.n2):
+            return
+
+        self.n1 = n1
+        self.n2 = n2
+        self.s_array = np.zeros(0, dtype=int)
+        self.configurations = np.zeros(0, dtype=np.uint64)
+
+    def reset(self):
+        self._reset(self.n1, self.n2)
+
+    def _reset(self, n1, n2):
+        self.n1 = None
+        self.n2 = None
+        self.set_shapes(n1, n2)
+
+    def pmf(self, k):
 
         # In practice, `pmf` is never called with k > m*n/2.
         # If it were, we'd exploit symmetry here:
@@ -35,143 +50,100 @@ def pmf(self, k, m, n):
         # i = k2 < k
         # k[i] = k2[i]
 
-        if (self._recursive is None and m <= 500 and n <= 500
-                or self._recursive):
-            return self.pmf_recursive(k, m, n)
-        else:
-            return self.pmf_iterative(k, m, n)
-
-    def pmf_recursive(self, k, m, n):
-        '''Probability mass function, recursive version'''
-        self._resize_fmnks(m, n, np.max(k))
-        # could loop over just the unique elements, but probably not worth
-        # the time to find them
-        for i in np.ravel(k):
-            self._f(m, n, i)
-        return self._fmnks[m, n, k] / special.binom(m + n, m)
-
-    def pmf_iterative(self, k, m, n):
-        '''Probability mass function, iterative version'''
-        fmnks = {}
-        for i in np.ravel(k):
-            fmnks = _mwu_f_iterative(m, n, i, fmnks)
-        return (np.array([fmnks[(m, n, ki)] for ki in k])
-                / special.binom(m + n, m))
-
-    def cdf(self, k, m, n):
+        pmfs = self.build_u_freqs_array(np.max(k))
+        return pmfs[k]
+
+    def cdf(self, k):
         '''Cumulative distribution function'''
 
         # In practice, `cdf` is never called with k > m*n/2.
         # If it were, we'd exploit symmetry here rather than in `sf`
-        pmfs = self.pmf(np.arange(0, np.max(k) + 1), m, n)
+        pmfs = self.build_u_freqs_array(np.max(k))
         cdfs = np.cumsum(pmfs)
         return cdfs[k]
 
-    def sf(self, k, m, n):
+    def sf(self, k):
         '''Survival function'''
         # Note that both CDF and SF include the PMF at k. The p-value is
         # calculated from the SF and should include the mass at k, so this
         # is desirable
 
-        # Use the fact that the distribution is symmetric; i.e.
-        # _f(m, n, m*n-k) = _f(m, n, k), and sum from the left
-        kc = np.asarray(m*n - k)  # complement of k
+        # Use the fact that the distribution is symmetric and sum from the left
+        kc = np.asarray(self.n1*self.n2 - k)  # complement of k
         i = k < kc
         if np.any(i):
             kc[i] = k[i]
-            cdfs = np.asarray(self.cdf(kc, m, n))
-            cdfs[i] = 1. - cdfs[i] + self.pmf(kc[i], m, n)
+            cdfs = np.asarray(self.cdf(kc))
+            cdfs[i] = 1. - cdfs[i] + self.pmf(kc[i])
         else:
-            cdfs = np.asarray(self.cdf(kc, m, n))
+            cdfs = np.asarray(self.cdf(kc))
         return cdfs[()]
 
-    def _resize_fmnks(self, m, n, k):
-        '''If necessary, expand the array that remembers PMF values'''
-        # could probably use `np.pad` but I'm not sure it would save code
-        shape_old = np.array(self._fmnks.shape)
-        shape_new = np.array((m+1, n+1, k+1))
-        if np.any(shape_new > shape_old):
-            shape = np.maximum(shape_old, shape_new)
-            fmnks = -np.ones(shape)             # create the new array
-            m0, n0, k0 = shape_old
-            fmnks[:m0, :n0, :k0] = self._fmnks  # copy remembered values
-            self._fmnks = fmnks
-
-    def _f(self, m, n, k):
-        '''Recursive implementation of function of [3] Theorem 2.5'''
-
-        # [3] Theorem 2.5 Line 1
-        if k < 0 or m < 0 or n < 0 or k > m*n:
-            return 0
-
-        # if already calculated, return the value
-        if self._fmnks[m, n, k] >= 0:
-            return self._fmnks[m, n, k]
-
-        if k == 0 and m >= 0 and n >= 0:  # [3] Theorem 2.5 Line 2
-            fmnk = 1
-        else:   # [3] Theorem 2.5 Line 3 / Equation 3
-            fmnk = self._f(m-1, n, k-n) + self._f(m, n-1, k)
-
-        self._fmnks[m, n, k] = fmnk  # remember result
-
-        return fmnk
-
-
-# Maintain state for faster repeat calls to mannwhitneyu w/ method='exact'
-_mwu_state = _MWU()
-
-
-def _mwu_f_iterative(m, n, k, fmnks):
-    '''Iterative implementation of function of [3] Theorem 2.5'''
-
-    def _base_case(m, n, k):
-        '''Base cases from recursive version'''
-
-        # if already calculated, return the value
-        if fmnks.get((m, n, k), -1) >= 0:
-            return fmnks[(m, n, k)]
-
-        # [3] Theorem 2.5 Line 1
-        elif k < 0 or m < 0 or n < 0 or k > m*n:
-            return 0
-
-        # [3] Theorem 2.5 Line 2
-        elif k == 0 and m >= 0 and n >= 0:
-            return 1
-
-        return None
-
-    stack = [(m, n, k)]
-    fmnk = None
-
-    while stack:
-        # Popping only if necessary would save a tiny bit of time, but NWI.
-        m, n, k = stack.pop()
-
-        # If we're at a base case, continue (stack unwinds)
-        fmnk = _base_case(m, n, k)
-        if fmnk is not None:
-            fmnks[(m, n, k)] = fmnk
-            continue
-
-        # If both terms are base cases, continue (stack unwinds)
-        f1 = _base_case(m-1, n, k-n)
-        f2 = _base_case(m, n-1, k)
-        if f1 is not None and f2 is not None:
-            # [3] Theorem 2.5 Line 3 / Equation 3
-            fmnk = f1 + f2
-            fmnks[(m, n, k)] = fmnk
-            continue
-
-        # recurse deeper
-        stack.append((m, n, k))
-        if f1 is None:
-            stack.append((m-1, n, k-n))
-        if f2 is None:
-            stack.append((m, n-1, k))
-
-    return fmnks
+    # build_sigma_array and build_u_freqs_array adapted from code
+    # by @toobaz with permission. Thanks to @andreasloe for the suggestion.
+    # See https://github.com/scipy/scipy/pull/4933#issuecomment-1898082691
+    def build_sigma_array(self, a):
+        n1, n2 = self.n1, self.n2
+        if a + 1 <= self.s_array.size:
+            return self.s_array[1:a+1]
+
+        s_array = np.zeros(a + 1, dtype=int)
+
+        for d in np.arange(1, n1 + 1):
+            # All multiples of d, except 0:
+            indices = np.arange(d, a + 1, d)
+            # \epsilon_d = 1:
+            s_array[indices] += d
+
+        for d in np.arange(n2 + 1, n2 + n1 + 1):
+            # All multiples of d, except 0:
+            indices = np.arange(d, a + 1, d)
+            # \epsilon_d = -1:
+            s_array[indices] -= d
+
+        # We don't need 0:
+        self.s_array = s_array
+        return s_array[1:]
+
+    def build_u_freqs_array(self, maxu):
+        """
+        Build all the array of frequences for u from 0 to maxu.
+        Assumptions:
+          n1 <= n2
+          maxu <= n1 * n2 / 2
+        """
+        n1, n2 = self.n1, self.n2
+        total = special.binom(n1 + n2, n1)
+
+        if maxu + 1 <= self.configurations.size:
+            return self.configurations[:maxu + 1] / total
+
+        s_array = self.build_sigma_array(maxu)
+
+        # Start working with ints, for maximum precision and efficiency:
+        configurations = np.zeros(maxu + 1, dtype=np.uint64)
+        # How many ways to have U=0? 1
+        configurations[0] = 1
+
+        with np.errstate(invalid='raise', over='raise'):
+            for u in np.arange(1, maxu + 1):
+                coeffs = s_array[u - 1::-1]
+                new_val = np.dot(configurations[:u], coeffs) / u
+                try:
+                    configurations[u] = new_val
+                except (FloatingPointError, OverflowError):
+                    # OK, we got into numbers too big for uint64.
+                    # So now we start working with floats.
+                    # By doing this since the beginning, we would have lost precision.
+                    # (And working on python long ints would be unbearably slow)
+                    configurations = configurations.astype(float)
+                    configurations[u] = new_val
+
+        self.configurations = configurations
+        return configurations / total
+
+
+_mwu_state = _MWU(0, 0)
 
 
 def _get_mwu_z(U, n1, n2, t, axis=0, continuity=True):
@@ -330,8 +302,8 @@ def mannwhitneyu(x, y, use_continuity=True, alternative="two-sided",
     consider `scipy.stats.wilcoxon`.
 
     `method` ``'exact'`` is recommended when there are no ties and when either
-    sample size is less than 8 [1]_. The implementation follows the recurrence
-    relation originally proposed in [1]_ as it is described in [3]_.
+    sample size is less than 8 [1]_. The implementation follows the algorithm
+    reported in [3]_.
     Note that the exact method is *not* corrected for ties, but
     `mannwhitneyu` will not raise errors or warnings if there are ties in the
     data. If there are ties and either samples is small (fewer than ~10
@@ -353,9 +325,9 @@ def mannwhitneyu(x, y, use_continuity=True, alternative="two-sided",
            Mathematical Statistics, Vol. 18, pp. 50-60, 1947.
     .. [2] Mann-Whitney U Test, Wikipedia,
            http://en.wikipedia.org/wiki/Mann-Whitney_U_test
-    .. [3] A. Di Bucchianico, "Combinatorics, computer algebra, and the
-           Wilcoxon-Mann-Whitney test", Journal of Statistical Planning and
-           Inference, Vol. 79, pp. 349-364, 1999.
+    .. [3] Andreas Löffler,
+           "Über eine Partition der nat. Zahlen und ihr Anwndung beim U-Test",
+           Wiss. Z. Univ. Halle, XXXII'83 pp. 87-89.
     .. [4] Rosie Shier, "Statistics: 2.3 The Mann-Whitney U Test", Mathematics
            Learning Support Centre, 2004.
     .. [5] Michael P. Fay and Michael A. Proschan. "Wilcoxon-Mann-Whitney
@@ -495,7 +467,8 @@ def mannwhitneyu(x, y, use_continuity=True, alternative="two-sided",
         method = _mwu_choose_method(n1, n2, np.any(t > 1))
 
     if method == "exact":
-        p = _mwu_state.sf(U.astype(int), min(n1, n2), max(n1, n2))
+        _mwu_state.set_shapes(n1, n2)
+        p = _mwu_state.sf(U.astype(int))
     elif method == "asymptotic":
         z = _get_mwu_z(U, n1, n2, t, continuity=use_continuity)
         p = stats.norm.sf(z)
diff --git a/scipy/stats/tests/test_hypotests.py b/scipy/stats/tests/test_hypotests.py
index c6f9adbcdc99..e229f89a3713 100644
--- a/scipy/stats/tests/test_hypotests.py
+++ b/scipy/stats/tests/test_hypotests.py
@@ -162,8 +162,6 @@ def test_callable_cdf(self):
 
 
 class TestMannWhitneyU:
-    def setup_method(self):
-        _mwu_state._recursive = True
 
     # All magic numbers are from R wilcox.test unless otherwise specified
     # https://rdrr.io/r/stats/wilcox.test.html
@@ -359,21 +357,23 @@ def test_exact_distribution(self):
             for m, p in table.items():
                 # check p-value against table
                 u = np.arange(0, len(p))
-                assert_allclose(_mwu_state.cdf(k=u, m=m, n=n), p, atol=1e-3)
+                _mwu_state.set_shapes(m, n)
+                assert_allclose(_mwu_state.cdf(k=u), p, atol=1e-3)
 
                 # check identity CDF + SF - PMF = 1
                 # ( In this implementation, SF(U) includes PMF(U) )
                 u2 = np.arange(0, m*n+1)
-                assert_allclose(_mwu_state.cdf(k=u2, m=m, n=n)
-                                + _mwu_state.sf(k=u2, m=m, n=n)
-                                - _mwu_state.pmf(k=u2, m=m, n=n), 1)
+                assert_allclose(_mwu_state.cdf(k=u2)
+                                + _mwu_state.sf(k=u2)
+                                - _mwu_state.pmf(k=u2), 1)
 
                 # check symmetry about mean of U, i.e. pmf(U) = pmf(m*n-U)
-                pmf = _mwu_state.pmf(k=u2, m=m, n=n)
+                pmf = _mwu_state.pmf(k=u2)
                 assert_allclose(pmf, pmf[::-1])
 
                 # check symmetry w.r.t. interchange of m, n
-                pmf2 = _mwu_state.pmf(k=u2, m=n, n=m)
+                _mwu_state.set_shapes(n, m)
+                pmf2 = _mwu_state.pmf(k=u2)
                 assert_allclose(pmf, pmf2)
 
     def test_asymptotic_behavior(self):
@@ -615,24 +615,25 @@ def test_gh19692_smaller_table(self):
         # sizes n1 and n2 changed. This was indicative of unnecessary cache
         # growth and redundant calculation. Check that this is resolved.
         rng = np.random.default_rng(7600451795963068007)
-        x = rng.random(size=5)
-        y = rng.random(size=11)
-        _mwu_state._fmnks = -np.ones((1, 1, 1))  # reset cache
-        stats.mannwhitneyu(x, y, method='exact')
-        shape = _mwu_state._fmnks.shape
-        assert shape[0] <= 6 and shape[1] <= 12  # one more than sizes
+        m, n = 5, 11
+        x = rng.random(size=m)
+        y = rng.random(size=n)
+        _mwu_state.reset()  # reset cache
+        res = stats.mannwhitneyu(x, y, method='exact')
+        shape = _mwu_state.configurations.shape
+        assert shape[-1] == min(res.statistic, m*n - res.statistic) + 1
         stats.mannwhitneyu(y, x, method='exact')
-        assert shape == _mwu_state._fmnks.shape  # unchanged when sizes are reversed
+        assert shape == _mwu_state.configurations.shape  # same when sizes are reversed
 
         # Also, we weren't exploiting the symmmetry of the null distribution
         # to its full potential. Ensure that the null distribution is not
         # evaluated explicitly for `k > m*n/2`.
-        _mwu_state._fmnks = -np.ones((1, 1, 1))  # reset cache
+        _mwu_state.reset()  # reset cache
         stats.mannwhitneyu(x, 0*y, method='exact', alternative='greater')
-        shape = _mwu_state._fmnks.shape
+        shape = _mwu_state.configurations.shape
         assert shape[-1] == 1  # k is smallest possible
         stats.mannwhitneyu(0*x, y, method='exact', alternative='greater')
-        assert shape == _mwu_state._fmnks.shape
+        assert shape == _mwu_state.configurations.shape
 
     @pytest.mark.parametrize('alternative', ['less', 'greater', 'two-sided'])
     def test_permutation_method(self, alternative):
@@ -646,43 +647,6 @@ def test_permutation_method(self, alternative):
         assert_allclose(res.statistic, res2.statistic, rtol=1e-15)
         assert_allclose(res.pvalue, res2.pvalue, rtol=1e-15)
 
-    def teardown_method(self):
-        _mwu_state._recursive = None
-
-
-class TestMannWhitneyU_iterative(TestMannWhitneyU):
-    def setup_method(self):
-        _mwu_state._recursive = False
-
-    def teardown_method(self):
-        _mwu_state._recursive = None
-
-
-@pytest.mark.xslow
-def test_mann_whitney_u_switch():
-    # Check that mannwhiteneyu switches between recursive and iterative
-    # implementations at n = 500
-
-    # ensure that recursion is not enforced
-    _mwu_state._recursive = None
-    _mwu_state._fmnks = -np.ones((1, 1, 1))
-
-    rng = np.random.default_rng(9546146887652)
-    x = rng.random(5)
-
-    # use iterative algorithm because n > 500
-    y = rng.random(501)
-    stats.mannwhitneyu(x, y, method='exact')
-    # iterative algorithm doesn't modify _mwu_state._fmnks
-    assert np.all(_mwu_state._fmnks == -1)
-
-    # use recursive algorithm because n <= 500
-    y = rng.random(500)
-    stats.mannwhitneyu(x, y, method='exact')
-
-    # recursive algorithm has modified _mwu_state._fmnks
-    assert not np.all(_mwu_state._fmnks == -1)
-
 
 class TestSomersD(_TestPythranFunc):
     def setup_method(self):


Code review comments:
diff_hunk: @@ -35,143 +50,100 @@ def pmf(self, k, m, n):
         # i = k2 < k
         # k[i] = k2[i]
 
-        if (self._recursive is None and m <= 500 and n <= 500
-                or self._recursive):
-            return self.pmf_recursive(k, m, n)
-        else:
-            return self.pmf_iterative(k, m, n)
-
-    def pmf_recursive(self, k, m, n):
-        '''Probability mass function, recursive version'''
-        self._resize_fmnks(m, n, np.max(k))
-        # could loop over just the unique elements, but probably not worth
-        # the time to find them
-        for i in np.ravel(k):
-            self._f(m, n, i)
-        return self._fmnks[m, n, k] / special.binom(m + n, m)
-
-    def pmf_iterative(self, k, m, n):
-        '''Probability mass function, iterative version'''
-        fmnks = {}
-        for i in np.ravel(k):
-            fmnks = _mwu_f_iterative(m, n, i, fmnks)
-        return (np.array([fmnks[(m, n, ki)] for ki in k])
-                / special.binom(m + n, m))
-
-    def cdf(self, k, m, n):
+        pmfs = self.build_u_freqs_array(np.max(k))
+        return pmfs[k]
+
+    def cdf(self, k):
         '''Cumulative distribution function'''
 
         # In practice, `cdf` is never called with k > m*n/2.
         # If it were, we'd exploit symmetry here rather than in `sf`
-        pmfs = self.pmf(np.arange(0, np.max(k) + 1), m, n)
+        pmfs = self.build_u_freqs_array(np.max(k))
         cdfs = np.cumsum(pmfs)
         return cdfs[k]
 
-    def sf(self, k, m, n):
+    def sf(self, k):
         '''Survival function'''
         # Note that both CDF and SF include the PMF at k. The p-value is
         # calculated from the SF and should include the mass at k, so this
         # is desirable
 
-        # Use the fact that the distribution is symmetric; i.e.
-        # _f(m, n, m*n-k) = _f(m, n, k), and sum from the left
-        kc = np.asarray(m*n - k)  # complement of k
+        # Use the fact that the distribution is symmetric and sum from the left
+        kc = np.asarray(self.n1*self.n2 - k)  # complement of k
         i = k < kc
         if np.any(i):
             kc[i] = k[i]
-            cdfs = np.asarray(self.cdf(kc, m, n))
-            cdfs[i] = 1. - cdfs[i] + self.pmf(kc[i], m, n)
+            cdfs = np.asarray(self.cdf(kc))
+            cdfs[i] = 1. - cdfs[i] + self.pmf(kc[i])
         else:
-            cdfs = np.asarray(self.cdf(kc, m, n))
+            cdfs = np.asarray(self.cdf(kc))
         return cdfs[()]
 
-    def _resize_fmnks(self, m, n, k):
-        '''If necessary, expand the array that remembers PMF values'''
-        # could probably use `np.pad` but I'm not sure it would save code
-        shape_old = np.array(self._fmnks.shape)
-        shape_new = np.array((m+1, n+1, k+1))
-        if np.any(shape_new > shape_old):
-            shape = np.maximum(shape_old, shape_new)
-            fmnks = -np.ones(shape)             # create the new array
-            m0, n0, k0 = shape_old
-            fmnks[:m0, :n0, :k0] = self._fmnks  # copy remembered values
-            self._fmnks = fmnks
-
-    def _f(self, m, n, k):
-        '''Recursive implementation of function of [3] Theorem 2.5'''
-
-        # [3] Theorem 2.5 Line 1
-        if k < 0 or m < 0 or n < 0 or k > m*n:
-            return 0
-
-        # if already calculated, return the value
-        if self._fmnks[m, n, k] >= 0:
-            return self._fmnks[m, n, k]
-
-        if k == 0 and m >= 0 and n >= 0:  # [3] Theorem 2.5 Line 2
-            fmnk = 1
-        else:   # [3] Theorem 2.5 Line 3 / Equation 3
-            fmnk = self._f(m-1, n, k-n) + self._f(m, n-1, k)
-
-        self._fmnks[m, n, k] = fmnk  # remember result
-
-        return fmnk
-
-
-# Maintain state for faster repeat calls to mannwhitneyu w/ method='exact'
-_mwu_state = _MWU()
-
-
-def _mwu_f_iterative(m, n, k, fmnks):
-    '''Iterative implementation of function of [3] Theorem 2.5'''
-
-    def _base_case(m, n, k):
-        '''Base cases from recursive version'''
-
-        # if already calculated, return the value
-        if fmnks.get((m, n, k), -1) >= 0:
-            return fmnks[(m, n, k)]
-
-        # [3] Theorem 2.5 Line 1
-        elif k < 0 or m < 0 or n < 0 or k > m*n:
-            return 0
-
-        # [3] Theorem 2.5 Line 2
-        elif k == 0 and m >= 0 and n >= 0:
-            return 1
-
-        return None
-
-    stack = [(m, n, k)]
-    fmnk = None
-
-    while stack:
-        # Popping only if necessary would save a tiny bit of time, but NWI.
-        m, n, k = stack.pop()
-
-        # If we're at a base case, continue (stack unwinds)
-        fmnk = _base_case(m, n, k)
-        if fmnk is not None:
-            fmnks[(m, n, k)] = fmnk
-            continue
-
-        # If both terms are base cases, continue (stack unwinds)
-        f1 = _base_case(m-1, n, k-n)
-        f2 = _base_case(m, n-1, k)
-        if f1 is not None and f2 is not None:
-            # [3] Theorem 2.5 Line 3 / Equation 3
-            fmnk = f1 + f2
-            fmnks[(m, n, k)] = fmnk
-            continue
-
-        # recurse deeper
-        stack.append((m, n, k))
-        if f1 is None:
-            stack.append((m-1, n, k-n))
-        if f2 is None:
-            stack.append((m, n-1, k))
-
-    return fmnks
+    # build_sigma_array and build_u_freqs_array adapted from code
+    # by @toobaz with permission. Thanks to @andreasloe for the suggestion.
+    # See https://github.com/scipy/scipy/pull/4933#issuecomment-1898082691
+    def build_sigma_array(self, a):
+        n1, n2 = self.n1, self.n2
+        if a + 1 <= self.s_array.size:
+            return self.s_array[1:a+1]
+
+        s_array = np.zeros(a + 1, dtype=int)
+
+        for d in np.arange(1, n1 + 1):
+            # All multiples of d, except 0:
+            indices = np.arange(d, a + 1, d)
+            # \epsilon_d = 1:
+            s_array[indices] += d
+
+        for d in np.arange(n2 + 1, n2 + n1 + 1):
+            # All multiples of d, except 0:
+            indices = np.arange(d, a + 1, d)
+            # \epsilon_d = -1:
+            s_array[indices] -= d
+
+        # We don't need 0:
+        self.s_array = s_array
+        return s_array[1:]
+
+    def build_u_freqs_array(self, maxu):
+        """
+        Build all the array of frequences for u from 0 to maxu.
+        Assumptions:
+          n1 <= n2
+          maxu <= n1 * n2 / 2
+        """
+        n1, n2 = self.n1, self.n2
+        total = special.binom(n1 + n2, n1)
+
+        if maxu + 1 <= self.configurations.size:
+            return self.configurations[:maxu + 1] / total
+
+        s_array = self.build_sigma_array(maxu)
+
+        # Start working with ints, for maximum precision and efficiency:
+        configurations = np.zeros(maxu + 1, dtype=np.uint64)
+        # How many ways to have U=0? 1
+        configurations[0] = 1
+
+        with np.errstate(invalid='raise', over='raise'):
path: scipy/stats/_mannwhitneyu.py
text: @author:
```suggestion
        with np.errstate(all='raise'):
```

diff_hunk: @@ -35,143 +50,100 @@ def pmf(self, k, m, n):
         # i = k2 < k
         # k[i] = k2[i]
 
-        if (self._recursive is None and m <= 500 and n <= 500
-                or self._recursive):
-            return self.pmf_recursive(k, m, n)
-        else:
-            return self.pmf_iterative(k, m, n)
-
-    def pmf_recursive(self, k, m, n):
-        '''Probability mass function, recursive version'''
-        self._resize_fmnks(m, n, np.max(k))
-        # could loop over just the unique elements, but probably not worth
-        # the time to find them
-        for i in np.ravel(k):
-            self._f(m, n, i)
-        return self._fmnks[m, n, k] / special.binom(m + n, m)
-
-    def pmf_iterative(self, k, m, n):
-        '''Probability mass function, iterative version'''
-        fmnks = {}
-        for i in np.ravel(k):
-            fmnks = _mwu_f_iterative(m, n, i, fmnks)
-        return (np.array([fmnks[(m, n, ki)] for ki in k])
-                / special.binom(m + n, m))
-
-    def cdf(self, k, m, n):
+        pmfs = self.build_u_freqs_array(np.max(k))
+        return pmfs[k]
+
+    def cdf(self, k):
         '''Cumulative distribution function'''
 
         # In practice, `cdf` is never called with k > m*n/2.
         # If it were, we'd exploit symmetry here rather than in `sf`
-        pmfs = self.pmf(np.arange(0, np.max(k) + 1), m, n)
+        pmfs = self.build_u_freqs_array(np.max(k))
         cdfs = np.cumsum(pmfs)
         return cdfs[k]
 
-    def sf(self, k, m, n):
+    def sf(self, k):
         '''Survival function'''
         # Note that both CDF and SF include the PMF at k. The p-value is
         # calculated from the SF and should include the mass at k, so this
         # is desirable
 
-        # Use the fact that the distribution is symmetric; i.e.
-        # _f(m, n, m*n-k) = _f(m, n, k), and sum from the left
-        kc = np.asarray(m*n - k)  # complement of k
+        # Use the fact that the distribution is symmetric and sum from the left
+        kc = np.asarray(self.n1*self.n2 - k)  # complement of k
         i = k < kc
         if np.any(i):
             kc[i] = k[i]
-            cdfs = np.asarray(self.cdf(kc, m, n))
-            cdfs[i] = 1. - cdfs[i] + self.pmf(kc[i], m, n)
+            cdfs = np.asarray(self.cdf(kc))
+            cdfs[i] = 1. - cdfs[i] + self.pmf(kc[i])
         else:
-            cdfs = np.asarray(self.cdf(kc, m, n))
+            cdfs = np.asarray(self.cdf(kc))
         return cdfs[()]
 
-    def _resize_fmnks(self, m, n, k):
-        '''If necessary, expand the array that remembers PMF values'''
-        # could probably use `np.pad` but I'm not sure it would save code
-        shape_old = np.array(self._fmnks.shape)
-        shape_new = np.array((m+1, n+1, k+1))
-        if np.any(shape_new > shape_old):
-            shape = np.maximum(shape_old, shape_new)
-            fmnks = -np.ones(shape)             # create the new array
-            m0, n0, k0 = shape_old
-            fmnks[:m0, :n0, :k0] = self._fmnks  # copy remembered values
-            self._fmnks = fmnks
-
-    def _f(self, m, n, k):
-        '''Recursive implementation of function of [3] Theorem 2.5'''
-
-        # [3] Theorem 2.5 Line 1
-        if k < 0 or m < 0 or n < 0 or k > m*n:
-            return 0
-
-        # if already calculated, return the value
-        if self._fmnks[m, n, k] >= 0:
-            return self._fmnks[m, n, k]
-
-        if k == 0 and m >= 0 and n >= 0:  # [3] Theorem 2.5 Line 2
-            fmnk = 1
-        else:   # [3] Theorem 2.5 Line 3 / Equation 3
-            fmnk = self._f(m-1, n, k-n) + self._f(m, n-1, k)
-
-        self._fmnks[m, n, k] = fmnk  # remember result
-
-        return fmnk
-
-
-# Maintain state for faster repeat calls to mannwhitneyu w/ method='exact'
-_mwu_state = _MWU()
-
-
-def _mwu_f_iterative(m, n, k, fmnks):
-    '''Iterative implementation of function of [3] Theorem 2.5'''
-
-    def _base_case(m, n, k):
-        '''Base cases from recursive version'''
-
-        # if already calculated, return the value
-        if fmnks.get((m, n, k), -1) >= 0:
-            return fmnks[(m, n, k)]
-
-        # [3] Theorem 2.5 Line 1
-        elif k < 0 or m < 0 or n < 0 or k > m*n:
-            return 0
-
-        # [3] Theorem 2.5 Line 2
-        elif k == 0 and m >= 0 and n >= 0:
-            return 1
-
-        return None
-
-    stack = [(m, n, k)]
-    fmnk = None
-
-    while stack:
-        # Popping only if necessary would save a tiny bit of time, but NWI.
-        m, n, k = stack.pop()
-
-        # If we're at a base case, continue (stack unwinds)
-        fmnk = _base_case(m, n, k)
-        if fmnk is not None:
-            fmnks[(m, n, k)] = fmnk
-            continue
-
-        # If both terms are base cases, continue (stack unwinds)
-        f1 = _base_case(m-1, n, k-n)
-        f2 = _base_case(m, n-1, k)
-        if f1 is not None and f2 is not None:
-            # [3] Theorem 2.5 Line 3 / Equation 3
-            fmnk = f1 + f2
-            fmnks[(m, n, k)] = fmnk
-            continue
-
-        # recurse deeper
-        stack.append((m, n, k))
-        if f1 is None:
-            stack.append((m-1, n, k-n))
-        if f2 is None:
-            stack.append((m, n-1, k))
-
-    return fmnks
+    # build_sigma_array and build_u_freqs_array adapted from code
+    # by @toobaz with permission. Thanks to @andreasloe for the suggestion.
+    # See https://github.com/scipy/scipy/pull/4933#issuecomment-1898082691
+    def build_sigma_array(self, a):
+        n1, n2 = self.n1, self.n2
+        if a + 1 <= self.s_array.size:
+            return self.s_array[1:a+1]
+
+        s_array = np.zeros(a + 1, dtype=int)
+
+        for d in np.arange(1, n1 + 1):
+            # All multiples of d, except 0:
+            indices = np.arange(d, a + 1, d)
+            # \epsilon_d = 1:
+            s_array[indices] += d
+
+        for d in np.arange(n2 + 1, n2 + n1 + 1):
+            # All multiples of d, except 0:
+            indices = np.arange(d, a + 1, d)
+            # \epsilon_d = -1:
+            s_array[indices] -= d
+
+        # We don't need 0:
+        self.s_array = s_array
+        return s_array[1:]
+
+    def build_u_freqs_array(self, maxu):
+        """
+        Build all the array of frequences for u from 0 to maxu.
+        Assumptions:
+          n1 <= n2
+          maxu <= n1 * n2 / 2
+        """
+        n1, n2 = self.n1, self.n2
+        total = special.binom(n1 + n2, n1)
+
+        if maxu + 1 <= self.configurations.size:
+            return self.configurations[:maxu + 1] / total
+
+        s_array = self.build_sigma_array(maxu)
+
+        # Start working with ints, for maximum precision and efficiency:
+        configurations = np.zeros(maxu + 1, dtype=np.uint64)
+        # How many ways to have U=0? 1
+        configurations[0] = 1
+
+        with np.errstate(invalid='raise', over='raise'):
+            for u in np.arange(1, maxu + 1):
+                coeffs = s_array[u - 1::-1]
+                new_val = np.dot(configurations[:u], coeffs) / u
+                try:
+                    configurations[u] = new_val
+                except (FloatingPointError, OverflowError):
path: scipy/stats/_mannwhitneyu.py
text: @author:
```suggestion
                except Exception:
```

diff_hunk: @@ -35,143 +50,100 @@ def pmf(self, k, m, n):
         # i = k2 < k
         # k[i] = k2[i]
 
-        if (self._recursive is None and m <= 500 and n <= 500
-                or self._recursive):
-            return self.pmf_recursive(k, m, n)
-        else:
-            return self.pmf_iterative(k, m, n)
-
-    def pmf_recursive(self, k, m, n):
-        '''Probability mass function, recursive version'''
-        self._resize_fmnks(m, n, np.max(k))
-        # could loop over just the unique elements, but probably not worth
-        # the time to find them
-        for i in np.ravel(k):
-            self._f(m, n, i)
-        return self._fmnks[m, n, k] / special.binom(m + n, m)
-
-    def pmf_iterative(self, k, m, n):
-        '''Probability mass function, iterative version'''
-        fmnks = {}
-        for i in np.ravel(k):
-            fmnks = _mwu_f_iterative(m, n, i, fmnks)
-        return (np.array([fmnks[(m, n, ki)] for ki in k])
-                / special.binom(m + n, m))
-
-    def cdf(self, k, m, n):
+        pmfs = self.build_u_freqs_array(np.max(k))
+        return pmfs[k]
+
+    def cdf(self, k):
         '''Cumulative distribution function'''
 
         # In practice, `cdf` is never called with k > m*n/2.
         # If it were, we'd exploit symmetry here rather than in `sf`
-        pmfs = self.pmf(np.arange(0, np.max(k) + 1), m, n)
+        pmfs = self.build_u_freqs_array(np.max(k))
         cdfs = np.cumsum(pmfs)
         return cdfs[k]
 
-    def sf(self, k, m, n):
+    def sf(self, k):
         '''Survival function'''
         # Note that both CDF and SF include the PMF at k. The p-value is
         # calculated from the SF and should include the mass at k, so this
         # is desirable
 
-        # Use the fact that the distribution is symmetric; i.e.
-        # _f(m, n, m*n-k) = _f(m, n, k), and sum from the left
-        kc = np.asarray(m*n - k)  # complement of k
+        # Use the fact that the distribution is symmetric and sum from the left
+        kc = np.asarray(self.n1*self.n2 - k)  # complement of k
         i = k < kc
         if np.any(i):
             kc[i] = k[i]
-            cdfs = np.asarray(self.cdf(kc, m, n))
-            cdfs[i] = 1. - cdfs[i] + self.pmf(kc[i], m, n)
+            cdfs = np.asarray(self.cdf(kc))
+            cdfs[i] = 1. - cdfs[i] + self.pmf(kc[i])
         else:
-            cdfs = np.asarray(self.cdf(kc, m, n))
+            cdfs = np.asarray(self.cdf(kc))
         return cdfs[()]
 
-    def _resize_fmnks(self, m, n, k):
-        '''If necessary, expand the array that remembers PMF values'''
-        # could probably use `np.pad` but I'm not sure it would save code
-        shape_old = np.array(self._fmnks.shape)
-        shape_new = np.array((m+1, n+1, k+1))
-        if np.any(shape_new > shape_old):
-            shape = np.maximum(shape_old, shape_new)
-            fmnks = -np.ones(shape)             # create the new array
-            m0, n0, k0 = shape_old
-            fmnks[:m0, :n0, :k0] = self._fmnks  # copy remembered values
-            self._fmnks = fmnks
-
-    def _f(self, m, n, k):
-        '''Recursive implementation of function of [3] Theorem 2.5'''
-
-        # [3] Theorem 2.5 Line 1
-        if k < 0 or m < 0 or n < 0 or k > m*n:
-            return 0
-
-        # if already calculated, return the value
-        if self._fmnks[m, n, k] >= 0:
-            return self._fmnks[m, n, k]
-
-        if k == 0 and m >= 0 and n >= 0:  # [3] Theorem 2.5 Line 2
-            fmnk = 1
-        else:   # [3] Theorem 2.5 Line 3 / Equation 3
-            fmnk = self._f(m-1, n, k-n) + self._f(m, n-1, k)
-
-        self._fmnks[m, n, k] = fmnk  # remember result
-
-        return fmnk
-
-
-# Maintain state for faster repeat calls to mannwhitneyu w/ method='exact'
-_mwu_state = _MWU()
-
-
-def _mwu_f_iterative(m, n, k, fmnks):
-    '''Iterative implementation of function of [3] Theorem 2.5'''
-
-    def _base_case(m, n, k):
-        '''Base cases from recursive version'''
-
-        # if already calculated, return the value
-        if fmnks.get((m, n, k), -1) >= 0:
-            return fmnks[(m, n, k)]
-
-        # [3] Theorem 2.5 Line 1
-        elif k < 0 or m < 0 or n < 0 or k > m*n:
-            return 0
-
-        # [3] Theorem 2.5 Line 2
-        elif k == 0 and m >= 0 and n >= 0:
-            return 1
-
-        return None
-
-    stack = [(m, n, k)]
-    fmnk = None
-
-    while stack:
-        # Popping only if necessary would save a tiny bit of time, but NWI.
-        m, n, k = stack.pop()
-
-        # If we're at a base case, continue (stack unwinds)
-        fmnk = _base_case(m, n, k)
-        if fmnk is not None:
-            fmnks[(m, n, k)] = fmnk
-            continue
-
-        # If both terms are base cases, continue (stack unwinds)
-        f1 = _base_case(m-1, n, k-n)
-        f2 = _base_case(m, n-1, k)
-        if f1 is not None and f2 is not None:
-            # [3] Theorem 2.5 Line 3 / Equation 3
-            fmnk = f1 + f2
-            fmnks[(m, n, k)] = fmnk
-            continue
-
-        # recurse deeper
-        stack.append((m, n, k))
-        if f1 is None:
-            stack.append((m-1, n, k-n))
-        if f2 is None:
-            stack.append((m, n-1, k))
-
-    return fmnks
+    # build_sigma_array and build_u_freqs_array adapted from code
+    # by @toobaz with permission. Thanks to @andreasloe for the suggestion.
+    # See https://github.com/scipy/scipy/pull/4933#issuecomment-1898082691
+    def build_sigma_array(self, a):
+        n1, n2 = self.n1, self.n2
+        if a + 1 <= self.s_array.size:
+            return self.s_array[1:a+1]
+
+        s_array = np.zeros(a + 1, dtype=int)
+
+        for d in np.arange(1, n1 + 1):
+            # All multiples of d, except 0:
+            indices = np.arange(d, a + 1, d)
+            # \epsilon_d = 1:
+            s_array[indices] += d
+
+        for d in np.arange(n2 + 1, n2 + n1 + 1):
+            # All multiples of d, except 0:
+            indices = np.arange(d, a + 1, d)
+            # \epsilon_d = -1:
+            s_array[indices] -= d
+
+        # We don't need 0:
+        self.s_array = s_array
+        return s_array[1:]
+
+    def build_u_freqs_array(self, maxu):
+        """
+        Build all the array of frequences for u from 0 to maxu.
+        Assumptions:
+          n1 <= n2
+          maxu <= n1 * n2 / 2
+        """
+        n1, n2 = self.n1, self.n2
+        total = special.binom(n1 + n2, n1)
+
+        if maxu + 1 <= self.configurations.size:
+            return self.configurations[:maxu + 1] / total
+
+        s_array = self.build_sigma_array(maxu)
+
+        # Start working with ints, for maximum precision and efficiency:
+        configurations = np.zeros(maxu + 1, dtype=np.uint64)
+        # How many ways to have U=0? 1
+        configurations[0] = 1
+
+        with np.errstate(invalid='raise', over='raise'):
path: scipy/stats/_mannwhitneyu.py
text: @author:
@user1 On some CI platforms, I don't this this is causing errors to be raised when overflow occurs. I think that that's all that's causing the CI failures because there were no failures when I made the `configurations` array `dtype=np.float64` at the beginning.
Do you have any suggestions?