読者です 読者をやめる 読者になる 読者になる

mayoko’s diary

プロコンとかいろいろ。

AOJ 2201 Immortal Jewels

解法

ある最適な直線があったとします。この直線を少し動かしても結果は変わりません。結果が変わる直線の動かし方は, その直線が ある円の内部を貫くか, ある円の中心からの距離が Ri + Mi となるときです。

なので, 考えるべき直線は, 二つの円の中心を考えたとき,

  • 円 i からは Ri 離れていて, 円 j からは Rj 離れているような直線
  • 円 i からは Ri 離れていて, 円 j からは Rj+Mj 離れているような直線
  • 円 i からは Ri+Mi 離れていて, 円 j からは Rj 離れているような直線
  • 円 i からは Ri+Mi 離れていて, 円 j からは Rj+Mj 離れているような直線

のみを考えれば良いです。

これを全部列挙すれば勝利ですが, 意外に直線を求めるのが面倒です。これを見ましょう。
shogo82148.github.io

国内予選でこれを出すのはなかなかハードですね…(考えの道筋が正しければ部分点をくれるような形式じゃないので)

typedef double Real;
Real eps = 1e-12;

Real add(Real a, Real b) {
    if (abs(a+b) < eps * (abs(a)+abs(b))) return 0;
    return a+b;
}

bool equal(Real a, Real b) {
    return add(a, -b) == 0;
}

struct P {
    Real x, y;
    P() {}
    P(Real x, Real y) : x(x), y(y) {}
    P operator+(P p) const {return P(add(x, p.x), add(y, p.y));}
    P operator-(P p) const {return P(add(x, -p.x), add(y, -p.y));}
    P operator*(Real d) const {return P(x*d, y*d);}
    Real dot(P p) const {return add(x*p.x, y*p.y);} // 内積
    Real det(P p) const {return add(x*p.y, -y*p.x);} // 外積
    Real dist(P p) const {return sqrt((x-p.x)*(x-p.x) + (y-p.y)*(y-p.y));} // 距離
    void normalize() {Real d = sqrt(x*x+y*y); x /= d; y /= d;} // 正規化
    bool operator<(const P& rhs) const {
        if (x != rhs.x) return x < rhs.x;
        return y < rhs.y;
    }
    bool operator==(const P& rhs) const {
        return equal(x, rhs.x) && equal(y, rhs.y);
    }
};

// 線分p1-p2上に点qがあるかを判定する
bool on_seg(P p1, P p2, P q) {
    return (p1-q).det(p2-q) == 0 && (p1-q).dot(p2-q) <= 0;
}

// 直線p1-p2と直線q1-q2が平行かどうかの判定
bool parallel(P p1, P p2, P q1, P q2) {
    P a = p2-p1;
    P b = q2-q1;
    return a.det(b) == 0;
}

// 直線p1-p2と直線q1-q2の交点
P intersection(P p1, P p2, P q1, P q2) {
    return p1+(p2-p1)*((q2-q1).det(q1-p1)/(q2-q1).det(p2-p1));
}
inline Real square(Real x) {return x*x;}
// 直線p1-p2と点qの距離
Real dist(P p1, P p2, P q) {
    q = q-p1;
    p2 = p2-p1;
    return sqrt((q.dot(q)*p2.dot(p2) - square(q.dot(p2))) / p2.dot(p2));
}

const int MAXN = 55;
int N;
P ps[MAXN];
Real R[MAXN], M[MAXN];

// ax+by+c = 0 の直線でカウント
int calc(Real a, Real b, Real c) {
	int ret = 0;
	for (int i = 0; i < N; i++) {
		double d = abs(a*ps[i].x + b*ps[i].y + c) / sqrt(a*a+b*b);
		if (d >= R[i]-eps && d <= R[i]+M[i]+eps) ret++;
	}
	return ret;
}

void getCoeff1(Real xp, Real yp, Real inSqrt, Real* xq, Real* yq, Real r, Real R) {
	*xq = r*(xp*(r+R)+yp*sqrt(inSqrt))/(xp*xp+yp*yp);
	*yq = r*(yp*(r+R)-xp*sqrt(inSqrt))/(xp*xp+yp*yp);
}
void getCoeff2(Real xp, Real yp, Real inSqrt, Real* xq, Real* yq, Real r, Real R) {
	*xq = r*(xp*(r+R)-yp*sqrt(inSqrt))/(xp*xp+yp*yp);
	*yq = r*(yp*(r+R)+xp*sqrt(inSqrt))/(xp*xp+yp*yp);
}

int ans;

void hoge(int i, int j, Real Rj) {
	Real xp = ps[j].x-ps[i].x;
	Real yp = ps[j].y-ps[i].y;
	Real tmp = (xp*xp+yp*yp)-(R[i]+Rj)*(R[i]+Rj);
	if (tmp >= 0) {
		double xq, yq;
		getCoeff1(xp, yp, tmp, &xq, &yq, R[i], Rj);
		ans = max(ans, calc(xq, yq, -xq*ps[i].x-yq*ps[i].y-R[i]*R[i]));
		getCoeff2(xp, yp, tmp, &xq, &yq, R[i], Rj);
		ans = max(ans, calc(xq, yq, -xq*ps[i].x-yq*ps[i].y-R[i]*R[i]));
	}
	tmp = (xp*xp+yp*yp)-(R[i]-Rj)*(R[i]-Rj);
	if (tmp >= 0) {
		double xq, yq;
		getCoeff1(xp, yp, tmp, &xq, &yq, R[i], -Rj);
		ans = max(ans, calc(xq, yq, -xq*ps[i].x-yq*ps[i].y-R[i]*R[i]));
		getCoeff2(xp, yp, tmp, &xq, &yq, R[i], -Rj);
		ans = max(ans, calc(xq, yq, -xq*ps[i].x-yq*ps[i].y-R[i]*R[i]));
	}
}

int main() {
	cin.tie(0);
	ios::sync_with_stdio(false);
	while (cin >> N) {
		if (N==0) break;
		for (int i = 0; i < N; i++)
			cin >> ps[i].x >> ps[i].y >> R[i] >> M[i];
		ans = 1;
		for (int i = 0; i < N; i++) for (int j = 0; j < N; j++) if (i!=j) {
			hoge(i, j, R[j]);
			hoge(i, j, R[j]+M[j]);
		}
		cout << ans << endl;
	}
	return 0;
}