Amazing Graze (AOJ 1023)

http://rose.u-aizu.ac.jp/onlinejudge/ProblemSet/description.jsp?id=1023&lang=jp

問題

略

解法

平面走査を行う。
各点のイベントとしてx-2*rの位置にinイベント、x+2*rの位置にoutイベントを挿入してx座標でソートしておく。
イベントが起こるとy座標でソートされた点の集合への挿入、削除が行わる。inイベントではその点から[y-2*r,y+2*r]の範囲に存在する全ての点について半径が4r以内に存在するかチェックして戦闘力の更新を行う。
範囲の取得には二部探索が使用できO(logn)で行える。一方、範囲に含まれる点の数は面積の関係でO(1)個しか存在しない。
全体でO(nlogn)で実行可能。
コードは続きから。

#include <stdio.h>
#include <string.h>
#include <algorithm>
#include <iostream>
#include <math.h>
#include <assert.h>
#include <vector>
#include <set>

using namespace std;
typedef long long ll;
static const double EPS = 1e-9;
static const double PI = acos(-1.0);

#define REP(i, n) for (int i = 0; i < (int)(n); i++)
#define FOR(i, s, n) for (int i = (s); i < (int)(n); i++)
#define FOREQ(i, s, n) for (int i = (s); i <= (int)(n); i++)
#define FORIT(it, c) for (__typeof((c).begin())it = (c).begin(); it != (c).end(); it++)
#define MEMSET(v, h) memset((v), h, sizeof(v))

struct Point {
  int x;
  int y;
  char type;
  Point() {;}
  Point(double x, double y, int type) : x(x), y(y), type(type) {;}
  bool operator<(const Point &rhs) const {
    if (y != rhs.y) { return y < rhs.y; }
    return x < rhs.x;
  }
};

struct Event {
  int x;
  Point p;
  char in;
  Event() {;}
  Event(int x, Point p, int in) : x(x), p(p), in(in) {;}
  bool operator<(const Event &rhs) const {
    if (x != rhs.x) { return x < rhs.x; }
    return in > rhs.in;
  }
};

inline int square(int x) { return x * x; }

inline int dist2(const Point &a, const Point &b) {
  return square(a.x - b.x) + square(a.y - b.y);
}

int an, bn, n;
int r;
Event event[500000];

int PlaneSweep() {
  int ans = 0;
  int r4 = 4 * r;
  int r2 = 16 * r * r;
  set<Point> open[2];
  REP(iter, n) {
    Event &e = event[iter];
    Point &p = e.p;
    if (e.in) {
      set<Point>::iterator lower = open[(int)p.type ^ 1].lower_bound(Point(0, p.y - r4, 0));
      set<Point>::iterator upper = open[(int)p.type ^ 1].upper_bound(Point(1000000000, p.y + r4, 0));
      for (set<Point>::iterator it = lower; it != upper; it++) {
        if (dist2(*it, p) <= r2) { ans++; }
      }
      open[(int)p.type].insert(p);
    } else {
      open[(int)p.type].erase(open[(int)p.type].find(p));
    }
  }
  return ans;
}

int main() {
  while (scanf("%d %d %d", &an, &bn, &r), an|bn) {
    n = (an + bn) * 2;
    REP(i, an) {
      int x, y;
      scanf("%d %d", &x, &y);
      event[2 * i + 0] = Event(x - 2 * r, Point(x, y, 0), 1);
      event[2 * i + 1] = Event(x + 2 * r, Point(x, y, 0), 0);
    }
    REP(i, bn) {
      int x, y;
      scanf("%d %d", &x, &y);
      event[2 * an + 2 * i + 0] = Event(x - 2 * r, Point(x, y, 1), 1);
      event[2 * an + 2 * i + 1] = Event(x + 2 * r, Point(x, y, 1), 0);
    }
    sort(event, event + n);
    int ans = PlaneSweep();
    printf("%d\n", ans);
  }
}