AOJ 1132 Circle and Points

問題

Circle and Points | Aizu Online Judge

解法

ある円が答えだったとします。この円を少し動かしても, この円が最適であることに変わりはありません。この円を適当に動かして, 初めてある点とぶつかったとき, 円をこの点からはみ出さないように動かすと, もう一つの点とぶつかります。

ということで, 円としては, 2 点を通るような円のみを考えれば良いです。

2 点 pnt[i], pnt[j] を通る円の中心はどうなるのかを考えます。この円の中心は, pnt[i], pnt[j] からの距離がどちらも 1 になっているようなものなので, pnt[i], pnt[j] を結ぶ線分の垂直二等分線上にあります。後は適当にホイしましょう。

typedef long double Real;

Real eps = 1e-9;

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);
    }
};
const int MAXN = 333;
P pnt[MAXN];
int N;

int calc(P center) {
	int ret = 0;
	for (int i = 0; i < N; i++) {
		P vec = pnt[i]-center;
		if (vec.dot(vec) < 1+eps) {
			ret++;
		}
	}
	return ret;
}

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 >> pnt[i].x >> pnt[i].y;
    	int ans = 1;
    	for (int i = 0; i < N; i++) for (int j = i+1; j < N; j++) {
    		if (pnt[i].dist(pnt[j]) > 2) continue;
    		// diff を取る
    		P vec = pnt[j]-pnt[i];
    		// vec に垂直なベクトルを取る
    		P n = P(vec.y, -vec.x);
    		n.normalize();
    		vec = vec*0.5;
    		Real len = sqrt(1-vec.dot(vec));
    		P cand1 = pnt[i] + vec + n*len;
    		P cand2 = pnt[i] + vec - n*len;
    		ans = max(ans, calc(cand1));
    		ans = max(ans, calc(cand2));
    	}
    	cout << ans << endl;
    }
    return 0;
}