code: def get_hypergeometric_stats(N, indices):
      assert isinstance(N, (int, np.integer))
      assert isinstance(indices, np.ndarray) and \
              np.issubdtype(indices.dtype, np.uint16)
      K = indices.size
      pvals = np.empty(N+1, dtype=np.float64)
      folds = np.empty(N+1, dtype=np.float64)
      pvals[0] = 1.0
      folds[0] = 1.0
      n = 0
      k = 0
      p = 1.0
      while n < N:
          if k < K and indices[k] == n:
              p *= (float((n+1) * (K-k)) / \
                    float((N-n) * (k+1)))
              k += 1
          else:
              p *= (float((n+1) * (N-K-n+k)) /
                    float((N-n) * (n-k+1)))
          n += 1
          pvals[n] = get_hgp(p, k, N, K, n)
          folds[n] = k / (K*(n/float(N)))
      return pvals, folds