Determinant of a matrix involving factorials Consider the matrix
$$M_n=\begin{bmatrix}1&\frac1{2!}&\frac1{3!}&\dotsb&\frac1{n!}\\
\frac1{2!}&\frac1{3!}&\frac1{4!}&\dotsb&\frac1{(n+1)!}\\
\frac1{3!}&\frac1{4!}&\frac1{5!}&\dotsb&\frac1{(n+2)!}\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
\frac1{n!}&\frac1{(n+1)!}&\frac1{(n+2)!}&\dotsb&\frac1{(2n-1)!}\end{bmatrix}$$
In other words, $M=[\frac1{(i+j-1)!}]_{1\le i,j\le n}$. Is there any formula for $\det M_n$?
Now, I know combinatorialists have lots of clever ways of evaluating determinants, so it doesn't seem so impossible to me that there might be a closed form solution to this. This is far from guaranteed, but I've been surprised before.
(By the way, I conjecture that $\det M_n\ne0$ for all $n$, but I don't immediately see a way to prove this.)
 A: This is not a complete answer but perhaps contributes another approach to help you discover / verify your formula for $\det M_n$.
Let us denote the value $p(a,n)$ to be the determinant of a similar looking $n\times n$ matrix starting with $1/a!$ instead on the top-left corner:
$$
p(a,n) = \det \begin{bmatrix} \frac{1}{a!} & \frac{1}{(a+1)!} & \cdots & \frac{1}{(a+n-1)!} \\ \frac{1}{(a+1)!} & \frac{1}{(a+2)!} & \cdots \\ \vdots  &  & \ddots \\ \frac{1}{(a+n-1)!} & \frac{1}{(a+n)!}  & \cdots & \frac{1}{(a+2n-2)!}\end{bmatrix}
$$
Then we have the recursive relation:
$$
p(a,n) p(a+2,n-2) = p(a+2,n-1)p(a,n-1) - p(a+1,n-1)^{2}  
$$
This comes from using the condensation relation for computing determinants:
Theorem. (Desnanot-Jacobi / Sylvester / Dodgson condensation)
Let $M$ be an $n\times n$ square matrix, and denote $M_{i}^{i'}$ to be the submatrix obtained from $M$ by removing the $i$-th row and $i'$-th column, and denote $M_{i,j}^{i',j'}$ to be the submatrix obtained from $M$ by removing the $i$-th and $j$-th row, and removing the $i'$-th and $j'$-th column. Then
$$
\det M \det M_{1,n}^{1,n} = \det M_{1}^{1}\det M_{n}^{n}-det M_{1}^{n}\det M_{n}^{1}
$$
(See: https://en.wikipedia.org/wiki/Dodgson_condensation, https://en.wikipedia.org/wiki/Sylvester%27s_determinant_identity)
Remark. This is often useful for finding $\det M$ when the entries of $M$ has a nice structure like in this case. I sometimes refer to this whimsically as "Alice in wonderland identity", as Charles Dodgson is more known as Lewis Carroll.
With a few initial conditions
$$
\begin{align*}
p(a,1) &= \det\left[ \frac{1}{a!}\right]= \frac{1}{a!}\\
p(a,2) &= \det\begin{bmatrix} \frac{1}{a!} & \frac{1}{(a+1)!}\\
\frac{1}{(a+1)!} & \frac{1}{(a+2)!}\end{bmatrix} = \frac{-1}{(a+1)!(a+2)!}
\end{align*}
$$
one could then compute a few more terms of $p(a,n)$ to guess a relation. By playing around a bit it does seem like we have $p(a,n) =\pm \frac{c(a,n)}{(a+n-1)!\cdots(a+2n-2)!}$ for some coefficient $c(a,n)$, though I did not work it out carefully. And finally to obtain your formula, we specialize $a=1$ to have $\det M_n = p(1,n)$.
Added edit.
By following this recursion, I got
$$
p(a,n) = \frac{c_n}{(a+n-1)!\cdots(a+2n-2)!},
$$
where $c_n = - \frac{c_{n-1}^2}{c_{n-2}}(n-1)$ with $c_1 = 1, c_2 = -1$.
This gives $c_n = (-1)^{\lfloor n/2\rfloor}(n-1)!\cdots 0!$ which gives
$$
p(a,n) = (-1)^{\lfloor n/2\rfloor}\frac{(n-1)!\cdots 0!}{(a+n-1)!\cdots(a+2n-2)!},
$$
which agrees with Misha's answer above.
A: Here is at least an easy proof that $\det M_n \neq 0$: $\det M_n$ is (up to sign) the Wronskian of the functions $\frac{x^n}{n!}, \frac{x^{n+1}}{(n+1)!}, \dots \frac{x^{2n-1}}{(2n-1)!}$, evaluated at $x = 1$. These functions are polynomials of distinct degree so they are clearly linearly independent; hence the Wronskian does not vanish identically. (This is not true in general but it is true for analytic functions.) On the other hand, by repeatedly factoring out all the factors of $x$ from both rows and columns, the Wronskian is some power of $x$ times a constant, namely its value at $x = 1$. So this constant, which is $\det M_n$ up to sign, also does not vanish.
A: Let $M_n'$ be the matrix obtained from $M_n$ by (1) multiplying the $i^{\text{th}}$ row by $(n+i-1)!$ and (2) dividing the $j^{\text{th}}$ column by $(n-j)!$. Then the $(i,j)$ entry of $M_n'$ is $\frac{(n+i-1)!}{(n-j)!(i+j-1)!} = \binom{n+i-1}{j+i-1}$. For example,
$$
   M_5' = \begin{bmatrix}
 5 & 10 & 10 & 5 & 1 \\
 15 & 20 & 15 & 6 & 1 \\
 35 & 35 & 21 & 7 & 1 \\
 70 & 56 & 28 & 8 & 1 \\
 126 & 84 & 36 & 9 & 1 \\
\end{bmatrix}.
$$
We see that each entry is the sum of the entry above it, plus the entry above and to its left. Formalizing this in matrix notation gives us the factorization $M_n' = A_n B_n$, where $a_{ij} = \binom{i-1}{j-1}$ and $b_{ij} = \binom{n}{i+j-1}$. For example,
$$
   M_5' = A_5 B_5 = \begin{bmatrix}
 1 & 0 & 0 & 0 & 0 \\
 1 & 1 & 0 & 0 & 0 \\
 1 & 2 & 1 & 0 & 0 \\
 1 & 3 & 3 & 1 & 0 \\
 1 & 4 & 6 & 4 & 1 \\
\end{bmatrix} \begin{bmatrix}
 5 & 10 & 10 & 5 & 1 \\
 10 & 10 & 5 & 1 & 0 \\
 10 & 5 & 1 & 0 & 0 \\
 5 & 1 & 0 & 0 & 0 \\
 1 & 0 & 0 & 0 & 0 \\
\end{bmatrix}.
$$
To prove that $M_n' = A_n B_n$ for all $n$, we need to check that the $(i,j)^{\text{th}}$ entry is equal to $\sum_{k=1}^n a_{ik} b_{kj}$, which means we must prove that $$\binom{n+i-1}{j+i-1} = \sum_{k=1}^n \binom{i-1}{k-1} \binom{n}{k+j-1}.$$ This is a special case of Vandermonde's identity.
It follows that $\det(M_n') = \det(A_n) \det(B_n)$. Here, $\det(A_n) = 1$, since $A_n$ is a lower triangular matrix. To put $B_n$ in the same form, we'd need to reflect it vertically, which means doing $\lfloor \frac n2\rfloor$ row swaps; therefore $\det(M_n') = \det(B_n) = (-1)^{\lfloor \frac n2\rfloor}$.
To go back to $\det(M_n)$, we need to undo the row operations we did when going from $M_n$ to $M_n'$. We divide the $i^{\text{th}}$ row by $(n+i-1)!$, which multiplies the determinant by $\prod_{i=1}^n \frac1{(n+i-1)!}$. We multiply the $j^{\text{th}}$ column by $(n-j)!$, which multiplies the determinant by $\prod_{j=1}^n (n-j)! = \prod_{i=1}^n (i-1)!$. Therefore $$\det(M_n) = (-1)^{\lfloor \frac n2\rfloor} \prod_{i=1}^n \frac{(i-1)!}{(n+i-1)!}.$$
Since each $\frac{(i-1)!}{(n+i-1)!}$ can be written as $\frac 1i \cdot \frac1{i+1} \cdots \frac1{(i+n-1)} = \prod_{j=1}^n \frac1{i+j-1}$, another handy way to write this expression is  $$\det(M_n) = (-1)^{\lfloor \frac n2\rfloor} \prod_{i=1}^n \prod_{j=1}^n \frac1{i+j-1},$$ which makes it clear that $\frac1{\det(M_n)}$ is an integer. In fact, the absolute value $|\frac1{\det(M_n)}|$ is sequence A107254 in the OEIS, which was what originally led me to guess this formula.
