As you can imagine, evaluating every pair of points is very expansive ($O(N^2)$ for $N$ points).
If your set of point is sufficiently dense with respect to $M$, a simple solution is to use a grid. Build a grid of $M \times M$ cells, then for each cell, compute the list of points inside (doing one loop on points $O(N)$). Finally you just have to compare every point with every other point in the same cell or in one of the 8 surrounding it.
If you set of points is very dense you even can use a second grid of $\frac{M}{\sqrt{2}} \times \frac{M}{\sqrt{2}}$ cells so the points in the same cells have necessarly pair-wise distance lower than $M$. This may save some computation.
Note that in the case you have lot of empty cells (small density), it may generate memory issue. Imagine for instance 2 pairs of points very distant to each other (with respect to $M$). It may require a tremendeous grid, just to process 4 points... In this case, the basic "all-pair" comparison would be largely more efficient. k-d tree may be a solution to this problem. Or you can either identify the empty zones to process the points subset by subset.