ARC068 E. Snuke Line - phyllo’s algorithm note

問題

直線上の鉄道で、0～Mまでの番号が付いたM+1個の駅がある。
N種類の名産品が与えられ、名産品iは駅lから駅rまでの区間で売られている。
最初0番の駅にいるとして、d駅ごとに停車する電車の場合、購入可能な名産品の種類数を知りたい。
dが1～Mの場合それぞれについて、種類数を答えよ。

制約

1 <= N <= 3*10^5
1 <= M <= 10^5
1 <= l_i <= r_i <= M

解説

AtCoder Regular Contest 068/ Beginner Contest 053 解説放送 - YouTube
解説放送の方の方法で解いてみた。

よく知られた事実として、以下のループはO(M log M)程度であることが使える。

for(int d=1; d<=M; d++){
  for(int p=1; d*(long long)p<=M; p++){
    //d*pについて、何かする
  }
}

単純な方法だと、「l_i <= d*p <= r_i」であるかをN個の名産品でチェックする方法がある。
しかし、これだとTLEしてしまうし、気を付けないと、d*pもd*(p+1)も範囲内であるような「l_i <= d*(p+1) <= r_i」のケースでダブルカウントしてしまうかもしれない。

範囲l_i～r_iに対して、d*pが1回しかカウントしないようにするうまい方法として、「d*(p-1)は範囲外で、d*pが範囲内の場合カウント」することで、範囲内の停車駅のうち、一番左の駅を見つけることができる。
X=d*(p-1)、Y=d*pとおいて考えると、「X < l_i、かつ、l_i <= Y <= r_i」を満たすような特産品を高速に求められれば良いことになる。

範囲は扱いにくいので、「l～r」という情報を2次元平面の点p(l,r)として考える。
上のXとYの条件はこの2次元平面上でどのような条件になるか考えると、分解して、

X < l
l <= Y
Y <= r

を満たす範囲の点の数を見つけることに相当する。
この3つの条件は、以下の*の部分のような平面上の長方形になっている。

M |-----+-------+
  |     |*******|
r |-----+***p***|
  |     |*******|
Y |-----+---+---|
  |     |   |   |
  +--------------
        X   l   Y

これを高速に処理することを考える。

今回、問い合わせクエリ(X,Y)と特産品の点(l,r)が先に求めておける。
そこで、問い合わせと特産品の点を合わせてl軸方向でソートしておき、r軸方向に対する1次元BITを使って、l軸の小さい方から点を置いたり、問い合わせに答えたりをすれば必要な(X,Y)の答えを求めておける。
あとは、X,Yの個数を使って該当する特産品数を求めて表示すればよい。

反省

(X,Y)の条件を満たす点の数を保存するのに、

map<pair<int,int>,int> memo;

としてしまったが、O(M log M)のループとはいえ、10^6個以上要素がある。
このままだと、memo[make_pair(X,Y)]などとアクセスした場合、木が深いため、時間が無視できない程度かかってしまう。
代わりに、

map<int,int> memo[100005];

などと木の深さが深くならないようにすることで、時間をかけずにアクセスできるようになる。

コード

#include <bits/stdc++.h>
using namespace std;
#define REP(i,a,n) for(int i=(a); i<(int)(n); i++)
#define rep(i,n) REP(i,0,n)
#define FOR(it,c) for(__typeof((c).begin()) it=(c).begin(); it!=(c).end(); ++it)
#define ALLOF(c) (c).begin(), (c).end()
typedef long long ll;
typedef unsigned long long ull;

#define MAX 100005

class BIT{
  static const int MAX_N = 1 << 18;
  int n, bit[MAX_N+1];
public:
  BIT(int n_){
    n = n_;
    for(int i=0; i<n+1; i++) bit[i] = 0;
  }
  int sum(int i){
    int s = 0;
    while(i>0){
      s+=bit[i];
      i-=i&(-i);
    }
    return s;
  }
  void add(int i, int x){
    while(i<=n){
      bit[i]+=x;
      i+=i&(-i);
    }
  }
};
 
struct ST {
  int type;
  int l, r;
};
bool operator<(const ST& a, const ST& b){
  if(a.l == b.l) return a.type < b.type;
  return a.l < b.l;
}
 
map<int,int> query[MAX];
 
int main(){
  ios::sync_with_stdio(false);
 
  int N, M;
  cin >> N >> M;
 
  vector<ST> pts;
  rep(i,N){
    int l, r;
    cin >> l >> r;
    pts.push_back((ST){0,l,r});
  }

  for(int d=2; d<=M; d++){
    for(int p=1; d*(ll)p<=M; p++){
      int X = d*(p-1);
      int Y = d*p;
 
      if(query[X].count(Y-1)==0){
        query[X][Y-1] = 1;
        pts.push_back((ST){1,X,Y-1});
      }
      if(query[X].count(MAX)==0){
        query[X][MAX] = 1;
        pts.push_back((ST){1,X,MAX});
      }
      if(query[Y].count(Y-1)==0){
        query[Y][Y-1] = 1;
        pts.push_back((ST){1,Y,Y-1});
      }
      if(query[Y].count(MAX)==0){
        query[Y][MAX] = 1;
        pts.push_back((ST){1,Y,MAX});
      }
    }
  }

  sort(ALLOF(pts));
 
  BIT bit(1<<18);
 
  rep(i,pts.size()){
    if(pts[i].type == 0){
      bit.add(pts[i].r+1, 1);
    }else{
      query[pts[i].l][pts[i].r] = bit.sum(pts[i].r+1);
    }
  }

  for(int d=1; d<=M; d++){
    if(d==1){
      cout << N << endl;
    }else{
      int ret = 0;
      for(int p=1; d*(ll)p<=M; p++){
        int X = d*(p-1);
        int Y = d*p;
 
        ret += query[Y][MAX] - query[Y][Y-1] - query[X][MAX] + query[X][Y-1];
      }
      cout << ret << endl;
    }
  }
 
  return 0;
}