Recursive formula for number of $n\times m\times p$ cubes of $0$'s and $1$'s with a special property.

One approach is classic inclusion-exclusion.

The argument is much simpler for $L_2(n,m,p)$. I'll use $1 \leq i \leq n$, $1 \leq j \leq m$, $1 \leq k \leq p$ throughout.

• Let $P_i$ be the set of such matrices with $B_{i,j,k} = 0$ for all $j,k$.
• Let $Q_j$ be the set of such matrices with $B_{i,j,k} = 0$ for all $i,k$.
• Let $R_k$ be the set of such matrices with $B_{i,j,k} = 0$ for all $i,j$.

$P_i$ forces $mp$ entries to be $0$, so $|P_i| = 2^{nmp - mp}$. Similarly an $a$-fold intersection of $P_i$'s does not force $(n-a)mp$ entries to be $0$, so the result has size $2^{(n-a)mp}$. Likewise, an $a$-fold intersection of $P_i$'s intersected with a $b$-fold intersection of $Q_j$'s and a $c$-fold intersection of $R_k$'s does not force $(n-a)(m-b)(p-c)$ entries to be $0$, so has size $2^{(n-a)(m-b)(p-c)}$.

By the principle of inclusion-exclusion, the number of such matrices in none of these sets is

$$L_2(n,m,p) = \sum_{a=0}^n \sum_{b=0}^m \sum_{c=0}^p (-1)^{a+b+c} \binom{n}{a} \binom{m}{b} \binom{p}{c} 2^{(n-a)(m-b)(p-c)}. \qquad(*)$$

Analogous reasoning gives $$L(n,m) = \sum_{a=0}^n \sum_{b=0}^m (-1)^{a+b} \binom{n}{a} \binom{m}{b} 2^{(n-a)(m-b)}. \qquad(**)$$

You can go from (**) to your formula by repeated applications of the binomial theorem: $$\begin{align*} \sum_{m=0}^M \binom{M}{m} L(N,m) &= \sum_{m=0}^M \binom{M}{m} \sum_{a=0}^N \sum_{b=0}^m (-1)^{a+b} \binom{N}{a} \binom{m}{b} 2^{(N-a)(m-b)} \\ &= \sum_{a=0}^N (-1)^a \binom{N}{a} \sum_{m=0}^M \binom{M}{m} (2^{N-a}-1)^m \\ &= \sum_{a=0}^N (-1)^a \binom{N}{a} (2^{N-a})^M \\ &= (2^M-1)^N. \end{align*}$$

Similarly from (*) we get $$\sum_{m=0}^M \sum_{p=0}^P \binom{M}{m} \binom{P}{p} L_2(N,m,p) = (2^{MP}-1)^N$$ since the left-hand side is $$\begin{align*} &\sum_{a=0}^N (-1)^a \binom{N}{a} \sum_{m=0}^M \binom{M}{m} \sum_{b=0}^m (-1)^b \binom{m}{b} \sum_{p=0}^P \binom{P}{c} \sum_{c=0}^p (-1)^c \binom{p}{c} 2^{(N-a)(m-b)(p-c)} \\ &= \sum_{a=0}^N (-1)^a \binom{N}{a} \sum_{m=0}^M \binom{M}{m} \sum_{b=0}^m (-1)^b \binom{m}{b} \sum_{p=0}^P \binom{P}{c} (2^{(N-a)(m-b)}-1)^p \\ &= \sum_{a=0}^N (-1)^a \binom{N}{a} \sum_{m=0}^M \binom{M}{m} \sum_{b=0}^m (-1)^b \binom{m}{b} (2^{P(N-m)})^{m-b} \\ &= \sum_{a=0}^N (-1)^a \binom{N}{a} \sum_{m=0}^M \binom{M}{m} (2^{P(N-a)}-1)^m \\ &= \sum_{a=0}^N (-1)^a \binom{N}{a} (2^{MP})^{N-a} \\ &= (2^{MP}-1)^N. \end{align*}$$

(Obviously this will all generalize to hypercubes.)

You can get a formula for $L_1(m, p, n)$ along these lines, but the multi-fold intersections of the analogues of the $P$'s, $Q$'s, and $R$'s are significantly more annoying, since lines can either be skew or intersect. You'd need to index the intersections by subsets of $(i, j)$'s, $(j, k)$'s, and $(k, i)$'s, tracking how many times these collections themselves intersect. One would expect a recursive formula to exist, but it may not be transparent how to get one from this approach. Perhaps I'm overly pessimistic, but I see no reason to expect there to be a "nice" answer for $L_1(n, m, p)$.