The docstring for special.factorial says that the function accepts ints, and int arrays.
However, it also seems to accept float arrays (including those that contain np.nan).
Is there a rationale for it accepting floats/float arrays? One can't calculate factorials of non-integers.
When `exact=False`, the default, the result is computed with a gamma function, which generalizes the factorial function to non-integers. The answer that you get is a perfectly reasonable, meaningful, non-arbitrary result.
When `exact=True`, you will see a `ValueError`, because we defer computation over to `math.factorial()` which does check for integralness because it's doing an exact computation that requires integralness.
I guess that we could add in a check for integralness in the `exact=False` case, but to what end? Would we catch anyone's real errors? There's actually a possible use case here for representing integral data that have missing data represented as NaNs (i.e. a float array where all the non-NaN entries are integral-valued despite being floating point numbers).