The following ideas are somewhat similar to that in Suvrit's answer.
Idea 1. Instead of first computing the square root of $A$, just solve the generalized eigenvalue problem $Ax = \lambda B x$ directly. This is especially tractable since according to the OP both $A$ and $B$ are symmetric and positive definite. You get a matrix $U$ such that $U^T B U = I$, $U^T A U = \mathrm{diag} (\lambda_i)$, and
\begin{align}
U^T (A + mB) U &= \mathrm{diag} (\lambda_i + m) , \\
U^{-1} (A + mB)^{-1} (U^{-1})^T &= \mathrm{diag} (\frac{1}{\lambda_i + m}) .
\end{align}
The formula that you want is then
$$
x^T (A + mB)^{-1} y = x^T U \mathrm{diag} (\frac{1}{\lambda_i + m}) U^{T} y .
$$
The cost of solving the generalized eigenvalue problem is still $O(n^3)$ and the cost of each of the above multiplications once $U$ is known is $O(n^2)$. Thus, the overall cost is $O(n^3 + Mn^2)$.
Idea 2. Now, an eigenvalue problem requires iteration. If you'd prefer an algorithm with a definite number of steps, you could instead simultaneously tridiagonalize $A$ and $B$, that is, find a matrix $U$ such that $U^T A U = T_A$ and $U^T B U = T_B$, where both $T_A$ and $T_B$ are tridiagonal. Tridiagonalization is in any case a common preconditioning step in eigenvalue problem algorithms. Similar to above, the formula that you want is then
$$
x^T (A + mB)^{-1} y = x^T U (T_A + m T_B)^{-1} U^{T} y .
$$
Since $T_A + m T_B$ is tridiagonal, you can compute its inverse in linear time, though slower by a constant factor than a diagonal matrix.
The cost of a tridiagonalization is once again $O(n^3)$ and the cost of each of the above inversions + multiplications is still $O(n^2)$ (dense matrix-vector multiplication dominates tridiagonal inversion). Thus, the overall cost is again $O(n^3 + Mn^2)$, though without requiring iteration.