but why the modulus?
The short answer is: because the modulus tells us whether one number divides another.
Some elaboration on that cannot hurt, of course. So let us elaborate. The notation
$$a \equiv b \pmod{c},$$
where $a,b,c$ are integers, by definition means that $c$ divides $a-b$, i.e. there is an integer $k$ such that $a-b = k\cdot c$. Hence $c$ divides $a$ if and only if $a \equiv 0 \pmod{c}$, and to see whether $i$ divides $n!$, we must check whether $n! \equiv 0 \pmod{i}$. The thing that makes this feasible for large $n$ (and $i > n$, it's trivially true for $i \leqslant n$) is that we do not need to compute $n!$ to see what remainder that would leave when divided by $i$. The property that allows this is that for all integers $m$ we have the implication
$$\bigl(a \equiv b \pmod{c}\bigr) \implies \bigl( ma \equiv mb \pmod{c}\bigr).\tag{$\ast$}$$
This is easy to see: since $a\equiv b \pmod{c}$ means $a-b = k\cdot c$ for some integer $k$, it follows that $ma - mb = m(a-b) = mk\cdot c$ for that same $k$, and $mk$ is of course an integer. So $c$ divides $ma - mb$, i.e. $ma \equiv mb \pmod{c}$.
And your algorithm is just applying $(\ast)$ at each step in the calculation of the factorial. Inductively, if $k! \equiv r_k \pmod{i}$, then by $(\ast)$ we have
$$(k+1)! = (k+1)\cdot k! \equiv (k+1)\cdot r_k \pmod{i}.$$
And then $r_{k+1}$ is the remainder of $(k+1)\cdot r_k$ upon division by $i$.
And clearly, if $r_k = 0$, then $r_{k+1} = 0$ too, so if we ever have a remainder of $0$ before reaching $n$, then we also have $r_n = 0$, i.e. $i$ divides $n!$.
Reducing modulo $i$ at each step solves the space problem, since none of the intermediate results is larger than $i^2$. But it's not the most efficient algorithm, since it typically requires $O(n)$ multiplications and divisions (reductions modulo $i$). It is more efficient to find the prime factorisation of $i$ (in the range up to $2^{32}$ that is still efficiently done by trial division), and then check whether each prime dividing $i$ occurs with a high enough exponent in the prime factorisation of $n!$. The exponent of a prime $p$ in the factorisation of $n!$ is quickly found using de Polignac's formula.
Best Answer
If you're extending factorials to real values, you'll probably want to look into the Gamma Function.
$\Gamma(n+1) = n!$ when $n$ is a natural number, but $\Gamma$ can be defined for all complex numbers except for the negative integers. In Python you can use the function
lgammafrom themathmodule in the standard library.lgammacomputes the log of the absolute value of the Gamma function.Something like this could work to compute the function you're interested in
though whether or not the Gamma function is actually what you want may be use case dependent.
Edit: I neglected to mention. If you plan to use Newton's method, you'll also need to be able to calculate the derivative of
lgamma. In Python you can use scipy.special.digamma from the 3rd party scipy package.