ボ〜ル (AOJ 2276)

http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=2276

問題

略

解法

まず地点を[0,PI]の区間に直す。次に明らかにいらない区間を全て削除する。後はDPをすれば良い。普通にDPするとO(N^2K)かかるが、その区間を選ぶよりも明らかに良い区間がある場合をスキップするとO(NK)になる。

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

using namespace std;
typedef long long ll;
typedef unsigned int uint;
typedef unsigned long long ull;
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))

#include <complex>
#include <vector>

typedef complex<double> Point;

static const double INF = 1e+1;

struct Line : public vector<Point> {
  Line() {;}
  Line(Point a, Point b) { push_back(a); push_back(b); }
};

struct Circle {
  Point p;
  double r;
  Circle() {;}
  Circle(Point p, double r) : p(p), r(r) {;}
};

vector<Point> tangentCP(const Circle &c, const Point &p) {
  vector<Point> ret;
  Point vect = c.p - p;
  double d = abs(vect);
  double l = sqrt(norm(vect) - c.r * c.r);
  if (isnan(l)) { return ret; }
  Point v1 = vect * Point(l / d,  c.r / d); 
  Point v2 = vect * Point(l / d, -c.r / d); 
  ret.push_back(v1);
  if (l > EPS) {
    ret.push_back(v2);
  }
  return ret;
}


struct Interval {
  double start;
  double end;
  Interval() {;}
  Interval(double s, double e) : start(s), end(e) {;}
  bool operator<(const Interval &rhs) const {
    return start < rhs.start;
  }
};

int n, k;
vector<Interval> interval;
map<double, double> memo[1550];

double calc(double start, int index, int rest) {
  if (start < interval[index].start) { start = interval[index].start; }
  if (memo[rest].count(start)) { return memo[rest][start]; }
  if (start >= PI || index == n || rest == 0) { return 0.0; }
  if (interval[index + 1].start <= start) { return calc(start, index + 1, rest); }
  double ret = calc(interval[index].end, index + 1, rest - 1) + interval[index].end - start;
  if (rest < n - index) {
    ret = max(ret, calc(start, index + 1, rest));
  }
  return memo[rest][start] = ret;
}

int main() {
  scanf("%d %d", &n, &k);
  Point center = Point(0, 0);
  REP(i, n) {
    double x, y, r;
    scanf("%lf %lf %lf", &x, &y, &r);
    Point p(x, y);
    Circle c = Circle(p, r);
    vector<Point> tangent = tangentCP(c, center);
    double start = arg(tangent[1]);
    double end = arg(tangent[0]);
    if (start < 0.0 && end < 0.0) { start = 0.0; end = 0.0; }
    if (start <= 0.0) { start = 0.0; }
    if (end < 0.0) { end = PI; }
    interval.push_back(Interval(start, end));
  }
  REP(i, interval.size()) {
    REP(j, interval.size()) {
      if (i == j) { continue; }
      if (interval[j].start <= interval[i].start && interval[i].end <= interval[j].end) {
        interval.erase(interval.begin() + i);
        i--;
        break;
      }
    }
  }
  n = interval.size();
  sort(interval.begin(), interval.end());
  interval.push_back(Interval(PI, PI));
  double ans = calc(0.0, 0, k) / PI;
  printf("%.8f\n", ans);
}