docs: def ks_unif_durbin_matrix(samples, statistic):
      """
      Calculates the probability that the statistic is less than the given value,
      using a fairly accurate implementation of the Durbin's matrix formula.
  
      Not an exact transliteration of the Marsaglia code, but using the same
      ideas. Assumes samples > 0. See: doi:10.18637/jss.v008.i18.
      """
      # Construct the Durbin matrix.
      h, k = modf(samples * statistic)
      k = int(k)
      h = 1 - h
      m = 2 * k + 1
      A = tri(m, k=1)
      hs = h ** arange(1, m + 1)
      A[:, 0] -= hs
      A[-1] -= hs[::-1]
      if h > .5:
          A[-1, 0] += (2 * h - 1) ** m
      A /= fromfunction(lambda i, j: gamma(fmax(1, i - j + 2)), (m, m))
      # Calculate A ** n, expressed as P * 2 ** eP to avoid overflows.
      P = identity(m)
      s = samples
      eA, eP = 0, 0
      while s != 1:
          s, b = divmod(s, 2)
          if b == 1:
              P = dot(P, A)
              eP += eA
              if P[k, k] > factor:
                  P /= factor
                  eP += shift
          A = dot(A, A)
          eA *= 2
          if A[k, k] > factor:
              A /= factor
              eA += shift
      P = dot(P, A)
      eP += eA
      # Calculate n! / n ** n * P[k, k].
      x = P[k, k]
      for i in arange(1, samples + 1):
          x *= i / samples
          if x < factorr:
              x *= factor
              eP -= shift
      return x * 2 ** eP