factorial(N), for scalar N, is the product of all the integers from 1 to N, i.e. prod(1:n). When N is an N-dimensional array, factorial(N) is the factorial for each element of N.
Since double precision numbers only have about 15 digits, the answer is only accurate for n <= 21. For larger n, the answer will have the right magnitude, and is accurate for the first 15 digits.