This problem can be rephrased as as a problem that, given \(N\) plots on a plane, enumerating the pair of points of distance at most \(K\), so we call a pair of integer \((x_i,y_i)\) as a point.

To tell the conclusion first, we can answer this problem by this solution.

1. For each point, assign \((x_i,y_i)\) to “Packet \((\lfloor \frac{x_i}{K} \rfloor , \lfloor \frac{y_i}{K} \rfloor)\).”
2. For each point, if the point is assigned to Packet \((X,Y)\), then the point to be paired can only be contained in the Packet itself or adjacent Packets \((X+l,Y+m)\) (where \(l\) and \(m\) are integers between \(-1\) and \(1\), inclusive), so exhaustively search for all points in these packets whether it can be paired.

The complexity of this solution is \(O(N+M)\) (ignoring the complexities like \(\log\) arising from the use of data structures), where \(M\) is the number of points to be output. How can we estimate it?

Let \(B_{X,Y}\) be the number of points in packet \((X,Y)\), and \(f(x)\) be “the least number of pairs of points of distance less than \(K\) in a single packet when the packet contains \(x\) points.”
First, we will show that \(f(x)\) is \(\Theta(x^2)\). It is obvious that \(f(x)\) is \(O(x^2)\), since we consider pairing \(2\) points out of \(x\) points. The proof that it is \(\Omega(x^2)\) (i.e. the lower bound is \(x^2\) (multiplied by a constant)) is given below.

In “a subsquare whose side’s length is less than \(\frac{K}{\sqrt{2}}\)”, for every pair of points contained in the subsquare, the distance between them is less than \(K\). Also, a \(K\) by \(K\) square can be covered with four square of length of side \(\frac{K}{\sqrt{2}}\), so we cover it with them. Here, if we plot \(4n\) number of points in a \(K\) by \(K\), at least one of the four small squares contains at least \(n\) points in it. Since the subsquare alone contains \(\frac{n(n-1)}{2}\) pairs of points of distance less than \(K\), so \(f(x)\) is \(\Omega(x^2)\).

Since \(\sum_{X,Y} f(B_{X,Y}) \le M\), the upper bound of \(\sum_{X,Y} B_{X,Y}^2\) is \(N+M\), multiplied by a constant.

Next, the number of pairs of points to be enumerated in this solution is \(\sum_{X,Y} \sum_{l=-1}^{1} \sum_{m=-1}^{1} B_{X,Y}B_{X+l,Y+m}\). Here, since the upper bound of \(\sum_{X,Y} B_{X,Y}^2\) is \(N+M\) multiplied by a constant, and \(B_{X,Y}B_{X+l,Y+m} \le \frac{1}{2}(B_{X,Y}^2+B_{X+l,Y+m}^2)\) due to inequality of arithmetic and geometric means, this can be bounded by \(N+M\), multiplied by the constant factor, too.

Therefore, we have proven that the complexity of this solution is \(O(N+M)\).

Sample code (C++):

#include<bits/stdc++.h>
#define big 2000000000

using namespace std;
using pi=pair<int,int>;
using pl=pair<long long,long long>;

vector<long long> dlt={-big-1,-big,-big+1,-1,0,1,big-1,big,big+1};

int main(){
  int n;
  long long k;
  cin >> n >> k;
  vector<pl> pts(n);
  map<long long,vector<int>> mp;
  for(int i=0;i<n;i++){
    cin >> pts[i].first >> pts[i].second;
    long long bkid=(pts[i].first/k)*big+(pts[i].second/k);
    mp[bkid].push_back(i);
  }
  vector<pi> res;
  for(int i=0;i<n;i++){
    long long bkid=(pts[i].first/k)*big+(pts[i].second/k);
    for(auto &d : dlt){
      long long cur=bkid+d;
      for(auto &nx : mp[cur]){
        if(i>=nx){continue;}
        long long dist=(pts[i].first-pts[nx].first)*(pts[i].first-pts[nx].first);
        dist+=(pts[i].second-pts[nx].second)*(pts[i].second-pts[nx].second);
        if(dist<=k*k){res.push_back({i+1,nx+1});}
      }
    }
  }
  sort(res.begin(),res.end());
  cout << res.size() << '\n';
  for(auto &nx : res){cout << nx.first << ' ' << nx.second << '\n';}
  return 0;
}