Stanley's Enumerative Combinatorics (http://www-math.mit.edu/~rstan/ec/ec1.pdf) contains next fact: 1.1.3 Example. Let f(n) be the number of n × n matrices M of $0$’s and $1$’s such that every row and column of M has three 1’s. For example, f(0) = 1, f(1) = f(2) = 0, f(3) = 1. The most explicit formula known at present for f(n) is $$f(n)=6^{-n}{(n!)}^2\sum\frac{(-1)^{\beta}(\beta+3\gamma)!2^\alpha 3^{\beta}}{\alpha!\beta!\gamma!^26^{\gamma}}$$ where the sum ranges over all (n + 2)(n + 1)/2 solutions to α + β + γ = n in nonnegative integers.
I need proof of this fact. (i.e. reference to book or articles that contains proof this fact).