The moving average filter (sometimes known colloquially as a boxcar filter) has a rectangular impulse response:
$$
h[n] = \frac{1}{N}\sum_{k=0}^{N-1} \delta[n-k]
$$
Or, stated differently:
$$
h[n] = \begin{cases}
\frac{1}{N}, && 0 \le n < N \\
0, && \text{otherwise}
\end{cases}
$$
Remembering that a discrete-time system's frequency response is equal to the discrete-time Fourier transform of its impulse response, we can calculate it as follows:
$$
\begin{align}
H(\omega) &= \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n} \\
&= \frac{1}{N}\sum_{n=0}^{N-1} e^{-j\omega n}
\end{align}
$$
To simplify this, we can use the known formula for the sum of the first $N$ terms of a geometric series:
$$
\sum_{n=0}^{N-1} e^{-j\omega n} = \frac{1-e^{-j \omega N}}{1 - e^{-j\omega}}
$$
What we're most interested in for your case is the magnitude response of the filter, $|H(\omega)|$. Using a couple simple manipulations, we can get that in an easier-to-comprehend form:
$$
\begin{align}
H(\omega) &= \frac{1}{N}\sum_{n=0}^{N-1} e^{-j\omega n} \\
&= \frac{1}{N} \frac{1-e^{-j \omega N}}{1 - e^{-j\omega}} \\
&= \frac{1}{N} \frac{e^{-j \omega N/2}}{e^{-j \omega/2}} \frac{e^{j\omega N/2} - e^{-j\omega N/2}}{e^{j\omega /2} - e^{-j\omega /2}}
\end{align}
$$
This may not look any easier to understand. However, due to Euler's identity, recall that:
$$
\sin(\omega) = \frac{e^{j\omega} - e^{-j\omega}}{j2}
$$
Therefore, we can write the above as:
$$
\begin{align}
H(\omega) &= \frac{1}{N} \frac{e^{-j \omega N/2}}{e^{-j \omega/2}} \frac{j2 \sin\left(\frac{\omega N}{2}\right)}{j2 \sin\left(\frac{\omega}{2}\right)} \\
&= \frac{1}{N} \frac{e^{-j \omega N/2}}{e^{-j \omega/2}} \frac{\sin\left(\frac{\omega N}{2}\right)}{\sin\left(\frac{\omega}{2}\right)}
\end{align}
$$
As I stated before, what you're really concerned about is the magnitude of the frequency response. So, we can take the magnitude of the above to simplify it further:
$$
|H(\omega)| = \frac{1}{N} \left|\frac{\sin\left(\frac{\omega N}{2}\right)}{\sin\left(\frac{\omega}{2}\right)}\right|
$$
Note: We are able to drop the exponential terms out because they don't affect the magnitude of the result; $|e^{j\omega}| = 1$ for all values of $\omega$. Since $|xy| = |x||y|$ for any two finite complex numbers $x$ and $y$, we can conclude that the presence of the exponential terms don't affect the overall magnitude response (instead, they affect the system's phase response).
The resulting function inside the magnitude brackets is a form of a Dirichlet kernel. It is sometimes called a periodic sinc function, because it resembles the sinc function somewhat in appearance, but is periodic instead.
Anyway, since the definition of cutoff frequency is somewhat underspecified (-3 dB point? -6 dB point? first sidelobe null?), you can use the above equation to solve for whatever you need. Specifically, you can do the following:
Set $|H(\omega)|$ to the value corresponding to the filter response that you want at the cutoff frequency.
Set $\omega$ equal to the cutoff frequency. To map a continuous-time frequency to the discrete-time domain, remember that $\omega = 2\pi \frac{f}{f_s}$, where $f_s$ is your sample rate.
Find the value of $N$ that gives you the best agreement between the left and right hand sides of the equation. That should be the length of your moving average.
