Factorial should be done using evaluate=False, but then numeric checking is nigh impossible because we sample in the half-open interval [0 - 1) which is not integer.

Currenly factorials above 20 are also too large integers to be cast into a datatype NumPy can deal with and so cause errors if evaluated by SymPy directly.