A New Ground Heating Device UVALive - 7675 （圆的面积交）

  1 #include <bits/stdc++.h>
  2 using namespace std;
  3 const double eps = 1e-8;
  4 const int maxn = 210;
  5 const double pi = acos(-1.0);
  6 double area[maxn];
  7 
  8 #define sqr(x) (x)*(x)
  9 int dcmp(double x) {
 10     if (x < -eps) return -1; else return x > eps;
 11 }
 12 struct cp {
 13     double x, y, r, angle;
 14     int d;
 15     cp(){}
 16     cp(double xx, double yy, double ang = 0, int t = 0) {
 17         x = xx;  y = yy;  angle = ang;  d = t;
 18     }
 19     void get() {
 20         scanf("%lf%lf%lf", &x, &y, &r);
 21         d = 1;
 22     }
 23 }cir[maxn], tp[maxn * 2];
 24 double dis(cp a, cp b) {
 25     return sqrt(sqr(a.x - b.x) + sqr(a.y - b.y));
 26 }
 27 double cross(cp p0, cp p1, cp p2) {
 28     return (p1.x - p0.x) * (p2.y - p0.y) - (p1.y - p0.y) * (p2.x - p0.x);
 29 }
 30 int CirCrossCir(cp p1, double r1, cp p2, double r2, cp &cp1, cp &cp2) {
 31     double mx = p2.x - p1.x, sx = p2.x + p1.x, mx2 = mx * mx;
 32     double my = p2.y - p1.y, sy = p2.y + p1.y, my2 = my * my;
 33     double sq = mx2 + my2, d = -(sq - sqr(r1 - r2)) * (sq - sqr(r1 + r2));
 34     if (d + eps < 0) return 0; if (d < eps) d = 0; else d = sqrt(d);
 35     double x = mx * ((r1 + r2) * (r1 - r2) + mx * sx) + sx * my2;
 36     double y = my * ((r1 + r2) * (r1 - r2) + my * sy) + sy * mx2;
 37     double dx = mx * d, dy = my * d; sq *= 2;
 38     cp1.x = (x - dy) / sq; cp1.y = (y + dx) / sq;
 39     cp2.x = (x + dy) / sq; cp2.y = (y - dx) / sq;
 40     if (d > eps) return 2; else return 1;
 41 }
 42 bool circmp(const cp& u, const cp& v) {
 43     return dcmp(u.r - v.r) < 0;
 44 }
 45 bool cmp(const cp& u, const cp& v) {
 46     if (dcmp(u.angle - v.angle)) return u.angle < v.angle;
 47     return u.d > v.d;
 48 }
 49 double calc(cp cir, cp cp1, cp cp2) {
 50     double ans = (cp2.angle - cp1.angle) * sqr(cir.r) 
 51         - cross(cir, cp1, cp2) + cross(cp(0, 0), cp1, cp2);
 52     return ans / 2;
 53 }
 54 void CirUnion(cp cir[], int n) {
 55     cp cp1, cp2;
 56     sort(cir, cir + n, circmp);
 57     for (int i = 0; i < n; ++i)
 58         for (int j = i + 1; j < n; ++j)
 59             if (dcmp(dis(cir[i], cir[j]) + cir[i].r - cir[j].r) <= 0)
 60                 cir[i].d++;
 61     for (int i = 0; i < n; ++i) {
 62         int tn = 0, cnt = 0;
 63         for (int j = 0; j < n; ++j) {
 64             if (i == j) continue;
 65             if (CirCrossCir(cir[i], cir[i].r, cir[j], cir[j].r,
 66                 cp2, cp1) < 2) continue;
 67             cp1.angle = atan2(cp1.y - cir[i].y, cp1.x - cir[i].x);
 68             cp2.angle = atan2(cp2.y - cir[i].y, cp2.x - cir[i].x);
 69             cp1.d = 1;    tp[tn++] = cp1;
 70             cp2.d = -1;   tp[tn++] = cp2;
 71             if (dcmp(cp1.angle - cp2.angle) > 0) cnt++;
 72         }
 73         tp[tn++] = cp(cir[i].x - cir[i].r, cir[i].y, pi, -cnt);
 74         tp[tn++] = cp(cir[i].x - cir[i].r, cir[i].y, -pi, cnt);
 75         sort(tp, tp + tn, cmp);
 76         int p, s = cir[i].d + tp[0].d;
 77         for (int j = 1; j < tn; ++j) {
 78             p = s;  s += tp[j].d;
 79             area[p] += calc(cir[i], tp[j - 1], tp[j]);
 80         }
 81     }
 82 }
 83 
 84 int X[maxn], Y[maxn];
 85 double Z[maxn];
 86 
 87 int n, k, w, s;
 88 int ck(double x){
 89     memset(cir, 0, sizeof(cir));
 90     for(int i = 0; i < n; i++){
 91         cir[i].x = X[i];
 92         cir[i].y = Y[i];
 93         cir[i].r = w * 1.0 / Z[i] / sqrt(sqr(X[i]) + sqr(Y[i]) + sqr(x));
 94         cir[i].d = 1;
 95     }
 96     memset(area, 0, sizeof(area));
 97     CirUnion(cir, n);
 98     if(dcmp(area[k] - s) >= 0) return 1;
 99     return 0;
100 }
101 int main(){
102     int t;
103     //freopen("in.txt", "r", stdin);
104     scanf("%d", &t);
105     while(t--){
106         scanf("%d %d %d %d", &n, &w, &k, &s);
107         for(int i = 0; i < n; i++){
108             scanf("%d %d %lf", &X[i], &Y[i], &Z[i]);
109         }
110         double L = 0, R = 500;
111         if(ck(500)) {
112             puts("Oops!");
113         }else if(!ck(0)) {
114             puts("No solution!");
115         }else {
116             while(R - L > 1e-6) {
117                 double m = (R + L) / 2;
118                 if(ck(m)) L = m;
119                 else R = m;
120             }
121             printf("%.4lf\n", R);
122         }
123         int ty;
124         while(scanf("%d", &ty) != EOF) cout<<area[ty]<<endl;
125     }
126     return 0;
127 }