This reads rather like a self-study question so I'll give an outline.
I'll assume that the intent was actually to consider $M^\top M + a^\top M$.
Consider the $i$-th component:
$M_i^2+a_iM_i = (M_i+a_i/2)^2 - a_i^2/4$
Now $\sum_i (M_i+a_i/2)^2$ will be a scaled non-central chi-squared
(To deal with the scale, let $Z=M/\sigma$ and let $a^*=a/\sigma$ and consider the above kind of manipulation in $Z$ and $a^*$ instead of $M$ and $a$; then the desired result is $\sigma^2$ times that one.)
Consequently, $M^\top M + a^\top M = \sum_i (M_i+a_i/2)^2 - \sum_i a_i^2/4$
will simply have a shifted and scaled non-central chi-squared distribution where the scale parameter, the noncentrality parameter and the shift parameter can be immediately written down by inspection.
[The scale and shift can be dealt with easily and R has the usual built-in d-,p-,q-, and r- functions for the non-central chi-square, so this leaves nothing more to do, aside the trivial details.]