We know that if three different points are not collinear then there is a unique circle that passes through them.
When three points are collinear we get a unique straight line passing through them and a straight line can be
considered as part of a circle of infinite radius. Given several thousand points and three special points
A, B and C in a two dimensional Cartesian coordinate system your job is to find out how
many points are contained in the circle passing through the three special points. Please note that your
solution must be very efficient.
In the figure above the three special points are shown as small filled circles and normal points are shown
as small filled squares. In Figure 1 we can see a circle passing through three special points A, B
and C. There is also a normal point D. There are total five points strictly within the circle and
three points are on the boundary. So the circle contains eight points in total. When the point D is
turned into special point C we get another circle, which is shown, in Figure 2. This larger circle
contains total nine points.
The input file contains at most ten sets of input. The description of each set is given below.
First line of each set contains six integers Ax, Ay,
Bx, By, N and Q. Here
(Ax, Ay) is the coordinate of A and
(Bx, By)
is the coordinate of B, N (
0 < N
50000) is the
total number of points (excluding A and B) and Q (
1Q1000) is the number of queries.
Each of the next N lines contains two integers xi, yi, where
(xi, yi) is the Cartesian
coordinate of the i-th point (
1iN,
-10000xi, yi10000). Each of the
next Q lines contains one integer SC (
1SCN), which indicates that
the SC-th point is to be considered as point C for the current query. For example, in the first
sample input the point with serial no 4 is (5, 5) as there are three points before
it ((-1, -1), (-1, -1) and (3,3)) in the input sequence.
Please note that for a single set of input, point A and B are fixed. Only the position of
point C is being altered for each query.
The input is terminated by a case where N = Q = 0. This case should not be processed. You can see from the
first sample input that more than one point can have the same coordinate. Such points should be considered as
different points.
For each set of input produce Q + 1 lines of output. The description of the output
for each set is given below:
The first line contains the serial number of the output as shown in the output for sample
input. Each of the next Q lines contains an integer D which indicates how many points will be within
the circle passing through A, B and C. The points on the boundary of the circle are also
considered inside the circle. So you should
also consider A, B and C within the circle. If the radius of the circle is greater
than 100000 then simply print the line `Impossible' instead of the integer D.
0 10 10 0 6 2
-1 -1
-1 -1
3 3
5 5
-1 -1
12 16
4
1
0 10 10 0 8 2
-1 -1
2 2
4 4
6 6
8 8
12 12
14 14
16 16
3
5
0 10 10 0 0 0
Case 1:
Impossible
7
Case 2:
8
7
將固定兩點的 AB 拉成一條直線，
然後討論每個點 P 對於過 AB 的圓，
過 AB 的圓與其半徑 r，會產生兩個圓，
假使圓心 Ci 與 P 同側(在 AB 的同左或者同右側)，盡可能讓 r 越小越好，
如果得知 Ci 與 P 同側，可以用二分搜尋找到所有比 query_r 還小的所有個數。
相反地，異側 r 則越大越好，可以找到所有比 query_r 還大的所有個數。
這樣對於每個 query 可以達到 O(logn)
在建立 Ci 與 P 討論的 r 時，可以利用二分搜尋來完成。
#include <stdio.h>
#include <math.h>
#include <algorithm>
#include <vector>
#define eps 1e-6
using namespace std;
struct Pt {
double x, y;
};
double dist(Pt &a, Pt &b) {
return pow(a.x-b.x,2)+pow(a.y-b.y,2);
}
Pt circle(Pt &a, Pt &b, Pt &c) {
static Pt ab, ac, o;
static double a1, b1, c1, a2, b2, c2, D, D1, D2;
ab.x = (a.x+b.x)/2, ab.y = (a.y+b.y)/2;
ac.x = (a.x+c.x)/2, ac.y = (a.y+c.y)/2;
a1 = a.x-b.x, b1 = a.y-b.y, c1 = a1*ab.x + b1*ab.y;
a2 = a.x-c.x, b2 = a.y-c.y, c2 = a2*ac.x + b2*ac.y;
D = a1*b2-a2*b1;
D1 = c1*b2-c2*b1, D2 = a1*c2-a2*c1;
o.x = D1/D, o.y = D2/D;
return o;
}
Pt A, B, C, mAB;
Pt p[50000];
int main() {
int cases = 0;
int n, q;
int i, j, k;
while(scanf("%lf %lf", &A.x, &A.y) == 2) {
scanf("%lf %lf", &B.x, &B.y);
scanf("%d %d", &n, &q);
if(n == 0 && q == 0) break;
for(i = 0; i < n; i++)
scanf("%lf %lf", &p[i].x, &p[i].y);
double a, b, c; // ax + by = c
a = A.y-B.y, b = -A.x+B.x;
if(a < 0) a = -a, b = -b;
c = a*A.x + b*A.y;
mAB.x = (A.x+B.x)/2, mAB.y = (A.y+B.y)/2;
double na = a, nb = b;
na /= sqrt(a*a+b*b);
nb /= sqrt(a*a+b*b);
if(na < 0) na = -na, nb = -nb;
vector<double> rSame, rDif, lSame, lDif;
int base = 0;
for(i = 0; i < n; i++) {
double argv = a*p[i].x+b*p[i].y-c;
if(fabs((A.x-B.x)*(p[i].y-A.y)-(A.y-B.y)*(p[i].x-A.x)) < eps) {
if(min(A.x, B.x) <= p[i].x && p[i].x <= max(A.x, B.x) &&
min(A.y, B.y) <= p[i].y && p[i].y <= max(A.y, B.y))
base++;
continue;
}
if(argv > 0) { // right side
double l = 0, r = 150000, m, cr2;
Pt C;
int cnt = 0;
while(l <= r && cnt < 50) {
m = (l+r)/2;
C.x = mAB.x + na*m, C.y = mAB.y + nb*m;
cr2 = dist(C, A);
if(dist(C, p[i]) <= cr2)
r = m;
else
l = m;
cnt++;
}
rSame.push_back(cr2);
// circle center in left side
l = 0, r = 150000, cnt = 0;
while(l <= r && cnt < 50) {
m = (l+r)/2;
C.x = mAB.x - na*m, C.y = mAB.y - nb*m;
cr2 = dist(C, A);
if(dist(C, p[i]) <= cr2)
l = m;
else
r = m;
cnt++;
}
rDif.push_back(cr2);
} else {// left side
double l = 0, r = 150000, m, cr2;
Pt C;
int cnt = 0;
while(l <= r && cnt < 50) {
m = (l+r)/2;
C.x = mAB.x + na*m, C.y = mAB.y + nb*m;
cr2 = dist(C, A);
if(dist(C, p[i]) <= cr2)
l = m;
else
r = m;
cnt++;
}
lDif.push_back(cr2);
// circle center in left side
l = 0, r = 150000, cnt = 0;
while(l <= r && cnt < 50) {
m = (l+r)/2;
C.x = mAB.x - na*m, C.y = mAB.y - nb*m;
cr2 = dist(C, A);
if(dist(C, p[i]) <= cr2)
r = m;
else
l = m;
cnt++;
}
lSame.push_back(cr2);
}
}
sort(rSame.begin(), rSame.end());
sort(rDif.begin(), rDif.end());
sort(lSame.begin(), lSame.end());
sort(lDif.begin(), lDif.end());
printf("Case %d:\n", ++cases);
while(q--) {
int Idx;
scanf("%d", &Idx);
C = p[Idx-1];
if(fabs((A.x-B.x)*(C.y-A.y)-(A.y-B.y)*(C.x-A.x)) < eps) {
puts("Impossible");
continue;
}
C = circle(A, B, C);
double argv = a*C.x + b*C.y - c, cr2;
cr2 = dist(C, A);
if(cr2 > (100000-eps)*(100000.0-eps)) {
puts("Impossible");
continue;
}
int v1, v2;
#define eps2 1e-3
if(argv > 0) { // center in right side
v1 = upper_bound(rSame.begin(), rSame.end(), cr2+eps2) - rSame.begin();
v2 = upper_bound(lDif.begin(), lDif.end(), cr2-eps2) - lDif.begin();
v2 = lDif.size()-v2;
} else { // center in left side
v1 = upper_bound(lSame.begin(), lSame.end(), cr2+eps2) - lSame.begin();
v2 = upper_bound(rDif.begin(), rDif.end(), cr2-eps2) - rDif.begin();
v2 = rDif.size()-v2;
}
printf("%d\n", v1+v2+2+base);
}
}
return 0;
}
文章定位：
