I don't have an answer to your question, but here are some thoughts.
If your normal distribution is able to generate a $B$ whose spectral radius is $\ge1$, then any procedure to modify $B$ to make its spectral radius less than one would clearly result in non-conformance to the normal distribution. So, I assume that you want to modify $B$ such that it follows a conditional distribution of normal distribution.
This seems to me a difficult problem, but the difficulty goes beyond that. In fact, when $p$ and $M$ are large, any such procedure would necessarily mean that you are drawing samples from the tails of normal distribution, i.e. you are mostly using extreme events as samples.
Here is an explanation. Let $v=(v_1^T,v_2^T,\ldots,v_n^T)^T$ be an eigenvector of $B$. If $v$ corresponds to the zero eigenvalue, then the equation $Bv=0$ means
$$
\begin{pmatrix}
A_1 & A_2 & \ldots & A_{p-1} & A_p\\
I & I & \ldots & I & 0
\end{pmatrix}v=0
$$
The rank of the matrix on the LHS is almost surely $2M$, so its null space is $(p-2)M$-dimensional and $B$ almost surely has $2M$ nonzero eigenvalues.
Now, suppose $\lambda$ is a nonzero eigenvalue of $B$. From $Bv=\lambda v$, we get
\begin{align}
A_1v_1 + \sum_{i=2}^p A_pv_p &= \lambda v_1,\\
v_1 + v_2+v_3+\ldots+v_{p-1} &= \lambda v_2,\\
v_1 + v_2+v_3+\ldots+v_{p-1} &= \lambda v_3,\\
&\vdots\\
v_1 + v_2+v_3+\ldots+v_{p-1} &= \lambda v_p.
\end{align}
In other words, $v_2=v_3=\ldots=v_p=$ some vector $u$, $v_1 = [\lambda-(p-2)]u$ and
\begin{align}
A_1[\lambda-(p-2)]u + \sum_{i=2}^p A_pu = \lambda [\lambda-(p-2)]u,\\
\textrm{i.e. }\quad
C_\lambda u:=\left(\lambda^2 I - \lambda [(p-2) I + A_1] + \left[(p-2)A_1 - \sum_{i=2}^p A_p\right]\right)u = 0.
\end{align}
Let $q(\lambda) = \det\,C_\lambda$. Then $q(\lambda)=0$ is a degree-$2M$ polynomial equation in $\lambda$. The sum $S$ of its $2M$ roots is the coefficient of $\lambda^{2M-1}$ in $q(\lambda)$, i.e. $S = \mathrm{trace}\left((p-2) I + A_1\right) = M(p-2)+\mathrm{trace}\,A_1$. In order that these $2M$ roots all have moduli smaller than $1$, $|S|$ must be smaller than $2M$. Hence $\left|(p-2)+\frac{\mathrm{trace}\,A_1}{M}\right|<2$. When $M$ is large, however $\frac{\mathrm{trace}\,A_1}{M}$ will be distributed around the mean $\mu$ of the normal distribution, with variance of order $O(1/M)$. Since $\left|(p-2)+\mu\right|\ge2$ when $p$ is large, we conclude that, when both $M$ and $p$ are large, the condition $|S|<2$ is satisfied only when $\frac{\mathrm{trace}\,A_1}{M}$ is an extreme value.
Let's call your matrix $A$. The matrix $A$ is real and symmetric and so it is similar to diagonal matrix with entries $\lambda_1, \dots, \lambda_n$. Also, $A$ has rank two so the diagonal matrix must also have rank two so let's say $\lambda_3 = \dots = \lambda_n = 0$. Now you need two equations to find $\lambda_1,\lambda_2$. A possible set of such equations is:
- $\operatorname{tr}(A) = \lambda_1 + \lambda_2 + \dots + \lambda_n = \lambda_1 + \lambda_2 = n$.
- $\operatorname{tr}(A^2) = \lambda_1^2 + \lambda_2^2 + \dots + \lambda_n^2 = \lambda_1^2 + \lambda_2^2.$
One readily computes $\operatorname{tr}(A^2)$ to be
$$ \operatorname{tr}(A^2) = 1^2 + 2^2 + \dots + (n-1)^2 + 1^2 + 2^2 + \dots + n^2 = \frac{n(2n^2+1)}{3}. $$
By plugging $\lambda_2 = n - \lambda_1$ into $\lambda_1^2 + \lambda_2^2 = \operatorname{tr}(A^2)$ you get a quadratic equation for $\lambda_1$ whose solutions are
$$ \frac{n}{2} \pm \sqrt{\frac{8n^3 - 6n^2 + 4n}{12}} $$
and those are the eigenvalues.
Best Answer
The definition of the characteristic polynomial is the same, with the obvious caveat that $$p_A(\lambda)\,\equiv\,det(\lambda\boldsymbol{I}\,-\,A)=\,\Pi_i(\lambda\,-\,\lambda_i)\,\mod\,2^{m+1}$$ where $2^{m+1}>\lambda_i\in\mathbb{Z}$ are the eigenvalues.
For you particular matrix, you may try to find a recurrence that will allow to get $A^{2^m-1}$ or diagonalize it first, exponentiate, add and change back to the original base.