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.
We start with a small example $2n+1=7$. Then we prove the formulas in the additive and multiplicative case for general odd positive numbers $2n+1$. We will then compare the tables and the proofs to see if there are some commonalities.
Additive case: $n=7$:
\begin{align*}
\begin{array}{rrrrrrr|r}
7&+12&+15&+16&+15&+12&+7\\
\hline
=7&+7&+7&+7&+7&+7&+7&=7\cdot 7\\
&+5&+5&+5&+5&+5&&+5\cdot 5\\
&&+3&+3&+3&&&+3\cdot 3\\
&&&+1&&&&+1\cdot1
\end{array}
\end{align*}
The table above shows how the summands in the first row can be split into equal parts of odd numbers so that the sum of these equal parts is a square number given in the rightmost column.
In general we obtain
\begin{align*}
\color{blue}{\sum_{j=-n}^n}\color{blue}{\sum_{k=0}^{n-|j|}(2n+1-2k)}
&=\sum_{{-n\le j\le n}\atop{0\le k\le n-|j|}}(2n+1-2k)\tag{1.1}\\
&=\sum_{{0\le k\le n}\atop{0\le |j|\le n-k}}(2n+1-2k)\tag{1.2}\\
&=\sum_{k=0}^n\sum_{j=-n+k}^{n-k}(2n+1-2k)\tag{1.3}\\
&=\sum_{k=0}^n(2n+1-2k)^2\tag{1.4}\\
&\,\,\color{blue}{=\sum_{k=0}^n(2k+1)^2}\tag{1.5}
\end{align*}
Comments:
In (1.1) we write the index region a bit more conveniently to better seethe transformations which follow.
In (1.2) we change the order of summation starting with $k$.
In (1.3) we use two summation symbols again.
In (1.4) we see the summands do not depend on $j$.
In (1.5) we change the order of summation $k\to n-k$.
Multiplicative case $2n+1=7$:
Now we take a look at the multplicative case and start again with a small table.
\begin{align*}
\begin{array}{ccccccc}
7&\cdot 12&\cdot 15&\cdot 16&\cdot 15&\cdot 12&\cdot 7\\
\hline
=(7)&\cdot(7+5)&\cdot(7+5+3)&\cdot(7+5+3+1)&\cdot(7+5+3)&\cdot(7+5)&\cdot(7)\\
=(7)&\cdot(6+6)&\cdot(5+5+5)&\cdot(4+4+4+4)&\cdot(5+5+5)&\cdot(6+6)&\cdot(7)\\
=(7)&\cdot(2\cdot 6)&\cdot(3\cdot 5)&\cdot(4\cdot4)&(3\cdot 5)&\cdot(2\cdot 6)&\cdot(7)\\
\end{array}
\end{align*}
The rearrangements in the table above are significantly different compared with the additive case. Here are the arrangements within each factor of the first row only.
In general we take again odd positive numbers $2n+1$ and consider
\begin{align*}
\color{blue}{\prod_{j=-n}^n\sum_{k=0}^{n-|j|}(2n+1-2k)}
\end{align*}
We observe in case $2n+1=7$ the greatest factor $16=7+5+3+1$ has four summands. It is convenient to split the calculation into $n=2m$ and $n=2m+1$ and perform the calculation for even and odd number of summands of the greatest factor separately. Here we consider the even case $n=2m$ only. The odd case can be calculated similarly.
We have
\begin{align*}
\color{blue}{\prod_{j=-2m}^{2m}\sum_{k=0}^{2m-|j|}(4m+1-2k)}\tag{2.1}
\end{align*}
Since the calculation of the inner sum of (2.1) is a bit different depending on the parity of $j$, we consider two cases.
Case $|j|=2q$: We obtain
\begin{align*}
\color{blue}{\sum_{k=0}^{2m-2q}}&\color{blue}{(4m+1-2k)}\\
&=\sum_{k=0}^{m-q-1}((4m+1-2k)+(4m+1-2(2m-2q-k))\\
&\qquad\qquad+4m+1-2(m-q)\tag{2.2}\\
&=\sum_{k=0}^{m-q-1}(4m+4q+2)+(2m+2q+1)\\
&=2\sum_{k=0}^{m-q-1}(2m+2q+1)+(2m+2q+1)\\
&=(2m-2q+1)(2m+2q+1)\tag{2.3}\\
&=(2m-|j|+1)(2m+|j|+1)\\
&\,\,\color{blue}{=(2m-j+1)(2m+j+1)}\tag{2.4}
\end{align*}
Comments:
Case $|j|=2q+1$: This case can be calculated analogously resulting in
\begin{align*}
\color{blue}{\sum_{k=0}^{2m-2q-1}(4m+1-2k)}&=(2m-2q)(2m+2q+2)\\
&\,\,\color{blue}{=(2m-j+1)(2m+j+1)}\tag{2.5}
\end{align*}
From (2.1),(2.4) and (2.5) we derive for $n=2m$:
\begin{align*}
\color{blue}{\prod_{j=-2m}^{2m}}&\color{blue}{\sum_{k=0}^{2m-|j|}(4m+1-2k)}\\
&=\prod_{j=-2m}^{2m}(2m-j+1)(2m+j+1)\\
&=\prod_{j=-2m}^{2m}(2m+j+1)^2\\
&=\prod_{j=0}^{4m}(j+1)^2\\
&\,\,\color{blue}{=((4m+1)!)^2}
\end{align*}
Conclusion: Both, the different rearrangements of the tables for the special case, as well as the different derivations of the general formulas rather indicate that there is no hidden relationship between the formulas in the additive and in the multiplicative case. Here it seems we have independent formulas showing us some nice properties of odd positive numbers.
Best Answer
Then your plot is wrong, or better yet, you just misread it or you mistyped what you wanted to say.
In general, you have
$$N! = 1\cdot 2\cdot 3\cdot 4\cdots\cdot N > 1\cdot 2\cdot 3\cdot 3\cdot 3\cdots 3 = 2\cdot 3^{N-2}$$
so $2^N$ will be smaller. Indeed, the same line of reasoning can be used to prove that $N!$ is bigger than $a^N$ for any $a\in\mathbb R$.
More precisely, the idea of the proof above shows the statement:
This also means that, for any $a\in\mathbb R^n$, we have
$$\lim_{n\to\infty}\frac{a^n}{n!} = 0$$ because, for large values of $n$, we have (if we take $b=a+1$):
$$\frac{a^n}{n!} < \frac{a^n}{c\cdot (a+1)^n} = \frac1c\cdot \left(\frac{a}{a+1}\right)^n,$$
and $$\lim_{n\to\infty} \left(\frac{a}{a+1}\right)^n = 0.$$