Generate algorithmically all grid points inside a hypercube
Here's another way of looking at the problem: You have a lattice generated by the columns of $M$. Use the Lenstra–Lenstra–Lovász (LLL) algorithm to obtain a reduced basis of this lattice. If you replace $M$ by a new matrix formed by the output of LLL, then the columns of $M$ will still generate the same lattice, but the basis vectors will be closer to being orthogonal to one another, and the entries of $M^{-1}$ should have smaller magnitude. 

From there, it would also help to bound each component of $v$ separately: i.e., you can bound the $i$th component $|v_i|$ by $\sum_{j=1}^d |(M^{-1})_{ij}|$. (By the way, the bound $\|v\|_\infty \leq \|M^{-1}\|$ is not correct; we need to use the sum of the elements on each row, not the maximum.)

For values of $d$ up to about 30, the LLL algorithm will finish practically instantly. Asymptotically, it takes $O(d^6)$, so it will slow down for very large $d$, but at the same time the number of points we need to check grows exponentially in $d$, so the run-time of the LLL is not really the bottleneck. On the other hand, the savings in the number of points needing to be checked can be enormous. I wrote some GAP code to generate a random regular (stochastic) matrix $M$ and compare the bounds on the components of $v$ that we obtain using the original basis, compared to the LLL-reduced basis (By the way, we do not need to assume that the matrix is regular; I made this restriction only because this was the case in your application):


d:=8;
M:=IdentityMat(d);
for i in [1..d] do
  for j in [1..d] do
    M[i][j]:=Random([-10^8..10^8]);
  od;
  M[i]:=M[i]/Sum(M[i]);
od;
L:=LLLReducedBasis(M).basis;
MM:=M^-1*1.0;
LL:=L^-1*1.0;
for i in [1..d] do
  for j in [1..d] do
    MM[i][j] := MM[i][j]*SignFloat(MM[i][j]);
    LL[i][j] := LL[i][j]*SignFloat(LL[i][j]);
  od;
od;

Print("Bounds for original basis: ");
ones:=[1..d]*0+1;
v:=MM*ones;
for i in [1..d] do
  v[i]:=Int(Floor(v[i]));
  Print(v[i]);
  Print(" ");
od;
Print("\n(");
Print(Product(v*2+1));
Print(" points to check)\n");

Print("Bounds for LLL basis: ");
v:=LL*ones;
for i in [1..d] do
  v[i]:=Int(Floor(v[i]));
  Print(v[i]);
  Print(" ");
od;
Print("\n(");
Print(Product(v*2+1));
Print(" points to check)\n");


The following output (based on the default random seed, with $d=8$) is not atypical:


Bounds for original basis: 9 23 24 4 23 16 23 4 
(258370076349 points to check)
Bounds for LLL basis: 3 3 2 2 3 4 2 3 
(2701125 points to check)


Edit: This problem is a special case of the general problem of enumerating lattice points in convex polytopes, which it turns out is a well-studied problem, and there are more efficient algorithms than the one described above. See this this paper for a survey.