So indeed it exists an efficient algorithm for computing the solution to this non-convex problem.

Let $\delta=2\pi/M$ and let $\angle$ denote the phase of a complex number (in radians). Upon defining $r_n=e^{i\phi_n}a_n$, the optimization problem can be recast as
\begin{equation}
 \max_{r_n} \left|\sum_{n=1}^N r_n\right|,
\end{equation}
subject to $|r_n|=|a_n|$ and $\angle r_n=\angle a_n+\delta \, k_n$ for some integer $k_n$ and for all $n$. If this problem is solved, then the solution to the original problem is $\phi_n=\delta \, k_n$ for all $n$.

It is not too hard to prove (left as an exercise to the reader) that any globoal maximizer satisfies
\begin{equation}
 \left|\angle r_n- \angle r_p\right|_{\text{wa}}\leq\delta
\end{equation}
for any $n,p$. Here, $|\angle r_n- \angle r_p|_{\text{wa}}$ is the wrap-around phase difference, so as an example, for $\angle r_n =\frac{3\pi}{2}$ and $\angle r_p =0$, the phase difference is $\frac{\pi}{2}$ and not $\frac{3\pi}{2}$.

The above conclusion can be expressed as
\begin{equation}
 \psi\leq \angle r_n \leq\psi+\delta,
\end{equation}
for an unknown $\psi$ and for all $n$ such that $a_n \neq 0$. Furthermore, it exists only one feasible value of $r_n$ satisfying the above condition (except if $\angle a_n = \psi$, in which case $\angle a_n = \psi+\delta$ is also feasible). Therefore, the problem boils down to a 1-dimensional search over $\psi$.

With this in mind, the global maximizer can be found with the following algorithm:

 1. let $\Omega_n$ be a permutation such that $0\leq \operatorname{mod}(\angle a_{\Omega_1},\delta)\leq\cdots\leq\operatorname{mod}(\angle a_{\Omega_N},\delta)$
 2. set $v_n$ such that $|v_n|=|a_{\Omega_n}|$ and $\angle v_n=\operatorname{mod}(\angle a_{\Omega_n},\delta)$ for $n=1,\ldots,N$
 3. initialize $S=0$
 4. **for** $p=1,\ldots,N$ **do**
 5. $\hspace{3ex}$ **if** $|\sum_{n=1}^N v_n|>S$ **then**
 6. $\hspace{6ex} S=|\sum_{n=1}^N v_n|$
 7. $\hspace{6ex} \phi_n = \angle v_n-\angle a_{\Omega_n}$
 8. $\hspace{3ex}$ **end if**
 9. $\hspace{3ex} \angle v_p \leftarrow  \angle v_p + \delta$
 10. **end for**