You can use the formula for geometric sum if the modular multiplicative inverse of the denominator of the geometric sum exists. Refer to @joffan 's answer.
Otherwise, you can solve this recursively. Note that:
1 + a + a^2 + a^3 + \dots + a^{2n+1} = (1 + a)(1 + (a^2) + (a^2)^2 + (a^2)^3 + \dots + (a^2)^n)
If we denote 1 + a + a^2 + a^3 + \dots + a^n as S(a,\ n), then
S(a, 2n + 1) = (1 + a) \cdot S(a^2, n)
S(a, 2n) = 1 + a \cdot (1 + a) \cdot S(a^2, n-1)
with base cases:
S(a, 0) = 1
S(a, 1) = 1 + a
Since at each step we are reducing the number of terms to half, the complexity is \mathcal{O}(\log{n}).
You just have to store the answer modulo M at every step.
Here is the Python
[1].
[1]: https://ideone.com/gWBvJl