全ての点が同じ座標にある場合、答えは \(N\) です。
そうでない場合、以下の補題が成り立ちます。

以下の条件を満たす最適解が存在する。

• 直線 \(l\) が以下の条件を満たすような整数 \(1 \le i < j \le N\) (ただし、点 \(i,j\) の座標は異なるとする) が存在する。
  • \(l\) と点 \(i,j\) を結んだ直線とが平行である。

上記の条件を満たさない最適解をひとつ選び出した時、直線を適切に回転・平行移動させることで上記の条件を満たすように変形できます( \(l\) から距離 \(D\) となる境界線のどちらかに点が \(2\) つ乗るまで適切に回転・平行移動を行うイメージです)。

傾きを \(O(N^2)\) 通り全探索することを考えます (縦軸に平行となるケースは \(x,y\) 座標を入れ替えて \(2\) 度解くことにして対処します)。
このとき、各傾き \(s\) に対して以下のように解けばよいです。

• 各点について、傾き \(s\) に沿って \(y\) 軸上 (\(x=0\)) に来るまで移動させる。
  • 傾きを決め打った時、このように移動させても領域内に包含されるかどうかは変わりません。
• 全ての点が \(y\) 軸上にあるため、点を \(y\) 座標の昇順にソートしてどの点からどの点まで含むことができるかを尺取り法で求める。

実装例の時間計算量は \(O(N^3 \log N)\) です。
また、実装例では全ての数値計算を有理数で済ませています。 long long では桁数が足りませんが、 __int128 を用いると十分な桁数を確保できます。

実装例 (C++):

#include<bits/stdc++.h>

using namespace std;
using ll=long long;
using i128=__int128;

struct frac{
  i128 p;
  i128 q;
  friend frac operator+(frac a, frac b) {
    frac res={a.p*b.q+a.q*b.p, a.q*b.q};
    if(res.q<0){res.p*=-1; res.q*=-1;}
    return res;
  }
  friend frac operator-(frac a, frac b) {
    frac res={a.p*b.q-a.q*b.p, a.q*b.q};
    if(res.q<0){res.p*=-1; res.q*=-1;}
    return res;
  }
  friend frac operator*(frac a, frac b) {
    frac res={a.p*b.p, a.q*b.q};
    if(res.q<0){res.p*=-1; res.q*=-1;}
    return res;
  }
  friend frac operator/(frac a, frac b) {
    frac res={a.p*b.q, a.q*b.p};
    if(res.q<0){res.p*=-1; res.q*=-1;}
    return res;
  }
  bool operator ==(const frac &r){
    return (((this->p)*r.q) == ((this->q)*r.p));
  }
  bool operator <(const frac &r){
    return (((this->p)*r.q) < ((this->q)*r.p));
  }
  bool operator <=(const frac &r){
    return (((this->p)*r.q) <= ((this->q)*r.p));
  }
};

struct pnt{
  frac x;
  frac y;
};

bool check(frac tg,frac D,pnt vec){
  tg=tg*tg;
  D=D*D;
  frac lef=vec.x*vec.x;
  frac rig=vec.y*vec.y;
  D=D/lef;
  D=D*(lef+rig);
  return (tg<=D);
}

ll solve(ll N,frac D,vector<pnt> X){
  ll res=1;
  {
    // dummy
    pnt del={{1ll,1ll},{1ll,1ll}};
    vector<frac> fv;
    for(ll k=0;k<N;k++){
      frac gap=(X[k].x/del.x)*del.y;
      fv.push_back(X[k].y-gap);
    }
    sort(fv.begin(),fv.end());
    ll r=0;
    for(ll l=0;l<N;l++){
      while(r<N && check(fv[r]-fv[l],D,del)){
        r++;
      }
      res=max(res,r-l);
    }
  }
  for(ll i=0;i<N;i++){
    for(ll j=i+1;j<N;j++){
      if(X[i].x==X[j].x){continue;}
      pnt del={X[j].x-X[i].x,X[j].y-X[i].y};
      vector<frac> fv;
      for(ll k=0;k<N;k++){
        frac gap=(X[k].x/del.x)*del.y;
        fv.push_back(X[k].y-gap);
      }
      sort(fv.begin(),fv.end());
      ll r=0;
      for(ll l=0;l<N;l++){
        while(r<N && check(fv[r]-fv[l],D,del)){
          r++;
        }
        res=max(res,r-l);
      }
    }
  }
  return res;
}

int main(){
  ll N,D;
  cin >> N >> D;
  vector<pnt> X(N);
  for(auto &nx : X){
    ll rx,ry;
    cin >> rx >> ry;
    nx.x.p=rx;
    nx.y.p=ry;
    nx.x.q=1; nx.y.q=1;
  }
  ll res=solve(N,frac{2*D,1},X);
  for(auto &nx : X){
    swap(nx.x,nx.y);
  }
  res=max(res,solve(N,frac{2*D,1},X));
  cout << res << "\n";
  return 0;
}