Below, I adapt this answer of usul to this problem to derive a different "explicit" formula for $g(m,n)$. Unfortunately, this is not a complete answer; I do not (yet?) see how to simplify the weighted sum over paths below to the conjectured result.

First note that the "collinear grid fillings" are in bijection with shellings of the complete bipartite graph $K_{m,n}$ viewed as a 1D simplicial complex. In this context a shelling is equivalent to an ordering of the edges $e_1,\dots,e_{mn}$ of $K_{m,n}$ so that the graph $G_j$ induced by the edge set $E_j=\{e_i|i\leq j\}$ is connected for all $1\leq j\leq mn$. We will write $V_j$ for the vertex set of $G_j$.

Let $p(m,n)=\frac{g(m,n)}{(mn)!}$ denote the probability that a given ordering of the edges of $K_{m,n}$ is a shelling.

The two vertex partitions of $K_{m,n}$ will be denoted A and B, respectively, so that $|A|=m,|B|=n$. For any valid shelling $(e_j)_{j=1}^{mn}$ of $K_{m,n}$ note that the set of ordered pairs:

$$(\text{number of A vertices in }V_j, \text{number of B vertices in }V_j)$$

as $j$ runs from 1 to $mn$, yields a subset of $\mathbb{Z}^2$ which form a walk on the standard grid graph from $(1,1)$ to $(m,n)$ taking steps only in the (1,0) and (0,1) directions. From now on, the term "walk" will refer to such "North or East" walks on the grid graph.

Let $\mathcal{R}=\{(k,l)\in\mathbb{Z}^2|1\leq k\leq m,1\leq l \leq n\}$ be the set of possible such A, B vertex counts reachable by a shelling of $K_{m,n}$. Let $d$ denote an auxiliary "dead" state.

Then $\mathcal{R}\cup\{d\}$ form the states of a Markov chain which coarse-grain the dynamics of edge addition. In addition to the (1,0) and (0,1) transitions which occur if new vertices are attached one at a time (thus preserving the shelling property), there are also transitions to $d$, which occur if an edge is added which disconnects $G_j$. $d$ and $(m,n)$ are the absorbing states.

We will compute the transition probabilities between states in this Markov chain which are induced by a random edge ordering; then $p(m,n)$ will be the probability that a walk from (1,1) survives to reach (m,n).

Suppose we consider a random edge ordering which is at the state $(k,l)\neq(m,n)$. As we add edges, the state must eventually change for some added edge $e_j$ and there are three possibilities:

- $e_j$ connects $V_{j-1}$ to A, ($l(m-k)$ choices of $e_j$)
- $e_j$ connects $V_{j-1}$ to B, ($k(n-l)$ choices)
- $e_j$ is not adjacent to any vertices in $V_{j-1}$ ($(m-k)(n-l)$ choices)

The total number of ways that we can leave $(k,l)$ is $mn-kl$ (equal to the number of edges of $K_{m,n}$ with at least one endpoint outside the $k$ vertices of A and $l$ vertices of B already reached). Thus the transition probabilities are:

- $\tilde{p}((k,l),(k+1,l))=\frac{l(m-k)}{mn-kl}$
- $\tilde{p}((k,l),(k,l+1))=\frac{k(n-l)}{mn-kl}$
- $\tilde{p}((k,l),d)=\frac{(m-k)(n-l)}{mn-kl}$

Here is an image showing the transition probabilities in the case $m=5,n=3$:

The probability that a random edge ordering induces a given walk $W$ from $(1,1)$ to $(m,n)$ is equal to the product of the transition probabilities $\tilde{p}(a,b)$ for all edges $(W_j,W_{j+1})$ of $W$. So for instance, the walk in the image above which goes north from $(1,1)$ to $(1,3)$ and then east to $(5,3)$ has probability $\frac{1}{7}\times\frac{1}{13}=\frac{1}{1001}$.

Finally, $p(m,n)$ is the sum of these probabilities over all possible walks. This can be written as:

$$p(m,n)=\sum_{W\text{ walk from }(1,1)\text{ to }(m,n)}\prod_{j=1}^{m+n-2}\tilde{p}(W_{j},W_{j+1}).$$

For the example above, $p(5,3)=\frac{1}{1001}+\frac{16}{7007}+\frac{32}{7007}+\frac{8}{1911}+\frac{16}{1911}+\frac{1}{91}+\frac{128}{21021}+\frac{256}{21021}+\frac{16}{1001}+\frac{16}{1001}+\frac{128}{21021}+\frac{256}{21021}+\frac{16}{1001}+\frac{16}{1001}+\frac{1}{91}=1/7.$

I have not yet managed to figure out how to simplify this sum in general; the conjecture gives $p(m,n)=\frac{m+n}{\binom{m+n}{m}}$. Curiously, $\binom{m+n}{n}$ is the number of walks from **(0,0)** to $(m,n)$.

Here is some SageMath code which implements this "weighted path sum" for $p(m,n)$ to compute $g(m,n)=(mn)!p(m,n)$:

```
def usulformula(m,n):
def weightA(k,l):
return l*(m-k)/(m*n-k*l)
def weightB(k,l):
return k*(n-l)/(m*n-k*l)
def recursewalk(k,l,w):
if k<m:
if l<n:
return (recursewalk(k+1,l,weightA(k,l)*w)
+recursewalk(k,l+1,weightB(k,l)*w))
else:
return recursewalk(k+1,l,weightA(k,l)*w)
else:
if l<n:
return recursewalk(k,l+1,weightB(k,l)*w)
else:
return w
return recursewalk(1,1,1)*factorial(m*n)
```

While this code is recursive, it only has two recursions as compared to the three recursions used in the `grid`

function in the question.

The values of $g(m,n)$ from the above code agree with those from the conjectured formula for all values of $m,n$ that I've tried.