Number of collinear ways to fill a grid A way to fill a finite grid (one box after the other) is called collinear if every newly filled box (the first excepted) is vertically or horizontally collinear with a previously filled box. See the following example:    

Let $g(m,n)$ be the number of collinear ways to fill a $m$-by-$n$ grid. Note that $g(m,n) = g(n,m)$.  
Question: What is an explicit formula for $g(m,n)$?
Conjecture (user44191): $g(m,n)=m!n!(mn)!/(m+n-1)!$.   
Definition (user44191): Let $g(m,n,i)$ be the number of collinear ways to fill $i$ boxes in a $m$-by-$n$ grid such that every row and every column contain at least one filled box.   
Remark: $g(m,n) = g(m,n,mn)$. 
Proposition (user44191): Here is a recursive formula for $g(m,n,i)$:   


*

*$g(1,1,1) = 1$.  

*If $m=0$ or $n=0$ or $ i< \min(m,n)$, then $g(m,n,i) = 0$.   

*$g(m,n,i+1)=(mn-i) g(m,n,i) + mn g(m-1,n,i) + mn g(m,n-1,i).$  


Proof: The two first points are obvious. We consider the number of collinear ways to fill $i+1$ boxes in a $m$-by-$n$ grid such that every row and every column contain at least one filled box.
There are three cases, corresponding to the three components of the recursive formula:  


*

*The last filled box is not the only filled box in its row and not the only filled box in its column.  

*The last filled box is the only filled box in its row.  

*The last filled box is the only filled box in its column.  


By the collinear assumption, 2. does not overlap 3.  $\square$  
One way to answer the question is to prove the conjecture using the above recursive formula.
We checked the conjecture for $1\le m \le n \le 5$, using the recursive formula (see below).
Remark: This question admits an extension to higher dimensional grids.
Remark: This question was inspired by that one. 

Sage program
# %attach SAGE/grid.sage

from sage.all import *

import copy

def grid(m,n,j):
    if [m,n,j]==[1,1,1]:
        return 1
    elif j < min(m,n) or m==0 or n==0:
        return 0
    else:
        i=j-1
        return (m*n-i)*grid(m,n,i) + m*n*grid(m-1,n,i) + m*n*grid(m,n-1,i)

def IsFormulaCorrect(m,n):
    return grid(m,n,m*n)==factorial(m)*factorial(n)*factorial(m*n)/factorial(m+n-1)

def CheckFormula(M,N):
    for m in range(1,M+1):
        for n in range(M,N+1):
            if not IsFormulaCorrect(m,n):
                return False
    return True

Computation 
sage: CheckFormula(5,5)
True

 A: Let me give an alternative proof with only finitely many variables used. It also gives some formulae for the following numbers: "when the filling becomes not collinear for the first time, there are $a$ columns and $b$ rows already used in the filling".
Consider a random enumeration $e_1,e_2,\ldots,e_{mn}$ of the edges of the graph $K_{n,m}$ with parts $N=\{v_1,\ldots,v_n\}$ and $M=\{u_1,\ldots,u_m\}$. We look it as a process: the edges $e_1,e_2,\ldots$ are added consecutively. Denote by $G(k)$ the graph formed by $e_1,\ldots,e_k$. For positive integers $a\leqslant n, b\leqslant m$ denote by $f(a,b)$ the probability of the following event $E(a,b)$:
(i) there exists $k$ such that the graphs $G(1),\ldots,G(k)$ are connected and the vertices of $G(k)$ are $v_1,\ldots,v_a;u_1,\ldots,u_b$; and
(ii) the vertices $v_1,v_2,\ldots,v_a$ appear in the sequence of graphs $G(1),G(2),\ldots$ in this order (that is, the first appearance of $v_i$ is before the first appearance of $v_{i+1}$, for each $i=1,\ldots,a-1$); the same for $u_1,\ldots,u_b$.
Clearly the probability of event (i) itself (without assuming (ii)) equals $a!b! f(a,b)$. So we should prove $f(n,m)=\frac1{(n+m-1)!}$. The key point is the following
Recursion. $f(1,1)=1$, also define $f(a,0)=f(0,b)=0$ for $n\geqslant a\geqslant 2$, $m\geqslant b\geqslant 2$. Then for $n\geqslant a\geqslant 1,m\geqslant b\geqslant 1$ and $(a,b)\ne (1,1)$ we have
$$
f(a,b)=\frac{a}{nm-a(b-1)}f(a,b-1)+\frac{b}{nm-(a-1)b}f(a-1,b). \,\quad(1)
$$
Proof of recursion. If the vertex $v_a$ appears before than $u_b$, then the event $E(a,b-1)$ happened. Consider the first moment after it happened and the new vertex (not $v_1,\ldots,v_a;u_1,\ldots,u_{b-1})$ appeared. It is the vertex $u_b$ with probability $\frac{a}{nm-a(b-1)}$ (there are $a$ admissible edges and $nm-a(b-1)$ edges using the new vertex in total). That's the first summand, the second analogously corresponds to the case when $u_b$ appears before $v_a$.
Now we solve recurrence (1). I do not know how it should be done in a regular way, so I guessed the answer and proved that it works.
Take a variable $t$ and call a double sequence of rational functions $\varphi(a,b)=\varphi(a,b)(t)$ defined for all non-negative integers $a,b$ appropriate if they satisfy a recurrence
$$(t-ab)\varphi (a,b)=a\varphi (a,b-1)+b\varphi (a-1,b)\quad \text{for}\,\,a,b\geqslant 1.$$
So admissible double sequences form a linear space.
Observation. For any $i=1,2,3,\ldots$ the function $$\psi_{i}(a,b):=(-1)^{b-i}\binom{a+b}{a+i}\frac1{(\frac{t}i+b)^{\underline{a+b+1}}}$$
is appropriate. Here we use Knuth's notation $x^{\underline{k}}=x(x-1)\ldots (x-k+1)$. The proof is straightforward.
Proposition. Consider the double sequence
$$
h(a,b):=\sum_{i=1}^\infty \left(1+\frac{t}{i^2}\right)\psi_i(a,b)
$$
(the summation is of course finite, only upto $i=b$, since for $i>b$ we have $\binom{a+b}{a+i}=0$). Then $h$ is an appropriate double sequence, satisfying $h(1,1)=\frac1{t(t-1)},h(a,0)=h(0,b)=0$ for $a,b\geqslant 2$.
Proof of the proposition. $h$ is appropriate by Observation and linearity of the set of appropriate functions. Relations $h(1,1)=\frac1{t(t-1)}$ and $h(a,0)=0$ are clear. For computing $h(0,b)$ we may look at $h(0,b)$ as a rational function of negative degree and check that all residues are equal to 0. The residue at $t=0$ equals $$\sum_{i=1}^b(-1)^{b-i}\binom{b}i\cdot i=b\sum_{i=1}^b(-1)^{b-i}\binom{b-1}{i-1}=0\quad \text{for}\quad b\geqslant 2.$$
The residues at $t=-ij$, $1\leqslant i,j\leqslant b, i\ne j$, for functions
$(1+\frac{t}{i^2})\psi_i(0,b)$ and $(1+\frac{t}{j^2})\psi_j(0,b)$ cancel out.
Main formula. $$f(a,b)=(t-ab)h(a,b)\quad \text{where} \quad t=mn$$
for $1\leqslant a\leqslant n, 1\leqslant b\leqslant m$. (For $a=n,b=m$ the function $h(n,m)$ has a pole at $t=mn$, but the product $(t-mn)h(n,m)$ is well-defined as ${\rm res}_{t=mn} h(n,m)$.)
Proof. Both parts satisfy the same initial conditions $f(1,1)=\frac1{mn}=\frac1{t}=(t-1)h(1,1)$, $f(a,0)=f(0,b)=0$ for $a\geqslant 2$ and the same recurrence (1).
Now for $a=n,b=m$ we get
$$f(n,m)={\rm res}_{t=mn}\sum_{i=1}^m\left(1+\frac{t}{i^2}\right)\psi_i(a,b)
={\rm res}_{t=mn} (1+\frac{n}m) \psi_m(n,m)=\frac{1}{(n+m-1)!}$$
A: Here is some further progress on the problem, without a full solution, based on ideas of Alex Postnikov. Alex's idea is to apply the "multivariate" (or "content-weighted") hook-length formula for Standard Young Tableaux. Everything (everything correct, that is) in this post is due to him.
Richard Stanley reduced the conjectured formula $g(m,n) = m!n!(mn)!/(m+n-1)!$ to showing that
$$ \sum_{\alpha} \frac{b^\alpha_1b^\alpha_2 \cdots b^\alpha_{m+n-2}}{b^\alpha_{m+n-2} (b^\alpha_{m+n-2}+b^\alpha_{m+n-3})\cdots (b^\alpha_{m+n-2}+b^\alpha_{m+n-3}+\cdots+b^\alpha_1) } = \frac{mn}{(m+n-1)!}$$
where the sum is over all $\alpha=(\alpha_1,\alpha_2,\ldots,\alpha_{m+n-2})$ consisting of $m-1$ 0's and $n-1$ 1's and where
$$ b^\alpha_i := \#\{1\leq j < i\colon \alpha_i\neq\alpha_j\} +1.$$
I don't know how to evaluate the above sum, but Alex's technique can be 
used to give a product formula for a very similar weighted sum over the same objects. Let us define
$$ c^\alpha_i := \#\{1\leq j \leq i\colon \alpha_i=\alpha_j\}.$$
Observe that $b^\alpha_i =(i+1)-c^\alpha_i$.
The claim is that
$$ \sum_{\alpha} \frac{1}{c^\alpha_{m+n-2} (c^\alpha_{m+n-2}+c^\alpha_{m+n-3})\cdots (c^\alpha_{m+n-2}+c^\alpha_{m+n-3}+\cdots+c^\alpha_1) } = \frac{2^{m-1}2^{n-1}}{(2m-2)! (2n-2)!}$$
where the sum is over the same $\alpha=(\alpha_1,\alpha_2,\ldots,\alpha_{m+n-2})$.
To see this, we biject such an $\alpha$ to a Standard Young Tableau of hook shape $(m,1^{n-1})$ in the obvious way: $i+1$ is in the arm of the hook iff $\alpha_i=0$. Then $c^\alpha_i$ is the absolute value of the content of the box containing $i+1$. And now we can directly apply the multivariate hook-length formula; see e.g. Theorem 18 of these notes: RSK via local transformations. In the notation of those notes we set $x_i := |i|$ for all $i \in \mathbb{Z}$.
As mentioned, I don't see immediately how this can be used to say anything about the relevant sum involving $b^{\alpha}_i$'s, but maybe someone else will figure out the approrpriate transformation.
A: 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.
A: By the method I used to solve
  Counting "connected" edge orderings (shellings) of the complete graph,
I can show the following:
  $$ g(m,n) = (mn-1)!\,m!\,n!\sum \frac{b_1 b_2\cdots b_{m+n-2}}
           {b_{m+n-2}(b_{m+n-3}+b_{m+n-2})\cdots(b_1+b_2+\cdots+b_{m+n-2})}, $$
where the sum is over all sequences $a_1 a_2\cdots a_{m+n-2}$ of $m-1$
$0$'s and $n-1$ $1$'s, and where
  $$ b_i=\#\{1\leq j<i\,\colon\, a_i\neq a_j\}+1. $$
Does anyone see why this is equal to $m!\,n!\,(mn)!/(m+n-1)!$?
A: We managed to obtain a solution via Stanley's comment above and some manipulations of binomial coefficients. See https://arxiv.org/abs/1809.10263
