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.$$