HDU 3467 (求五个圆相交面积) Song of the Siren

  1 #include <cstdio>
  2 #include <cmath>
  3 #include <cstring>
  4 #include <algorithm>
  5 
  6 const int maxn = 10;
  7 const double eps = 1e-8;
  8 const double PI = acos(-1.0);
  9 
 10 int dcmp(double x)
 11 { return (x > eps) - (x < -eps); }
 12 
 13 struct Point
 14 {
 15     double x, y;
 16     Point(double x=0, double y=0):x(x), y(y) {}
 17     void read() { scanf("%lf%lf", &x, &y); }
 18 };
 19 typedef Point Vector;
 20 Point operator + (const Vector& a, const Vector& b)
 21 { return Point(a.x+b.x, a.y+b.y); }
 22 Point operator - (const Vector& a, const Vector& b)
 23 { return Point(a.x-b.x, a.y-b.y); }
 24 Vector operator * (const Vector& a, double p)
 25 { return Point(a.x*p, a.y*p); }
 26 Vector operator / (const Vector& a, double p)
 27 { return Point(a.x/p, a.y/p); }
 28 bool operator == (const Point& a, const Point& b)
 29 { return dcmp(a.x-b.x) == 0 && dcmp(a.y-b.y) == 0; }
 30 
 31 double Dot(const Vector& a, const Vector& b)
 32 { return a.x*b.x + a.y*b.y; }
 33 double Cross(const Vector& a, const Vector& b)
 34 { return a.x*b.y - a.y*b.x; }
 35 double Length(const Vector& a)
 36 { return sqrt(Dot(a, a)); }
 37 Vector unit(const Vector& a)
 38 { return a / Length(a); }
 39 Vector Normal(const Vector& a)
 40 {
 41     double l = Length(a);
 42     return Vector(-a.y/l, a.x/l);
 43 }
 44 double Angle(const Vector& a)
 45 { return atan2(a.y, a.x); }
 46 
 47 Point Rotate(const Point& p, double angle, const Point& o = Point(0, 0))
 48 {
 49     Vector t = p - o;
 50     t = Vector(t.x*cos(angle)-t.y*sin(angle), t.x*sin(angle)+t.y*cos(angle));
 51     return t + o;
 52 }
 53 
 54 struct Region
 55 {
 56     double st, ed;
 57     Region(double s=0, double e=0):st(s), ed(e) {}
 58 };
 59 
 60 struct Circle
 61 {
 62     Point c;
 63     double r;
 64     Circle() {}
 65     Circle(Point c, double r):c(c), r(r) {}
 66     
 67     void read() { c.read(); scanf("%lf", &r); }
 68     
 69     double area() const { return PI * r * r; }
 70     
 71     bool contain(const Circle& rhs) const
 72     { return dcmp(Length(c-rhs.c) + rhs.r - r) <= 0; }
 73     
 74     bool contain(const Point& p) const
 75     { return dcmp(Length(c-p) - r) <= 0; }
 76     
 77     bool intersect(const Circle& rhs) const
 78     { return dcmp(Length(c-rhs.c) - r - rhs.r) < 0; }
 79     
 80     bool tangency(const Circle& rhs) const
 81     { return dcmp(Length(c-rhs.c) - r - rhs.r) == 0; }
 82     
 83     Point get_point(double ang) const
 84     { return Point(c.x + r * cos(ang), c.y + r * sin(ang)); }
 85 };
 86 
 87 void IntersectionPoint(const Circle& c1, const Circle& c2, Point& p1, Point& p2)
 88 {
 89     double d = Length(c1.c - c2.c);
 90     double l = (c1.r*c1.r + d*d - c2.r*c2.r) / (2 * d);
 91     double h = sqrt(c1.r*c1.r - l*l);
 92     Point mid = c1.c + unit(c2.c-c1.c) * l;
 93     Vector t = Normal(c2.c - c1.c) * h;
 94     p1 = mid + t;
 95     p2 = mid - t;
 96 }
 97 
 98 double IntersectionArea(const Circle& c1, const Circle& c2)
 99 {
100     double area = 0.0;
101     const Circle& M = c1.r > c2.r ? c1 : c2;
102     const Circle& N = c1.r > c2.r ? c2 : c1;
103     double d = Length(c1.c-c2.c);
104     
105     if(d < M.r + N.r && d > M.r - N.r)
106     {
107         double Alpha = 2.0 * acos((M.r*M.r + d*d - N.r*N.r) / (2 * M.r * d));
108         double Beta  = 2.0 * acos((N.r*N.r + d*d - M.r*M.r) / (2 * N.r * d));
109         area = ( M.r*M.r*(Alpha - sin(Alpha)) + N.r*N.r*(Beta - sin(Beta)) ) / 2.0;
110     }
111     else if(d <= M.r - N.r) area = N.area();
112     
113     return area;
114 }
115 
116 struct Region_vector
117 {
118     int n;
119     Region v[5];
120     void clear() { n = 0; }
121     void add(const Region& r) { v[n++] = r; }
122 } *last, *cur;
123 
124 Circle cir[maxn];
125 bool del[maxn];
126 double r;
127 int n = 5;
128 
129 bool IsOnlyOnePoint()
130 {
131     bool flag = false;
132     Point t;
133     for(int i = 0; i < n; ++i)
134     {
135         for(int j = i + 1; j < n; ++j)
136         {
137             if(cir[i].tangency(cir[j]))
138             {
139                 t = (cir[i].c + cir[j].c) / 2;
140                 flag = true;
141                 break;
142             }
143         }
144     }
145     
146     if(!flag) return false;
147     for(int i = 0; i < n; ++i)
148         if(!cir[i].contain(t)) return false;
149     
150     printf("Only the point (%.2f, %.2f) is for victory.\n", t.x, t.y);
151     return true;
152 }
153 
154 bool solve()
155 {
156     if(IsOnlyOnePoint()) return true;
157     memset(del, false, sizeof(del));
158     
159     for(int i = 0; i < n; ++i)
160         for(int j = 0; j < n; ++j)
161         {
162             if(del[j] || i == j) continue;
163             if(cir[i].contain(cir[j]))
164             {
165                 del[i] = true;
166                 break;
167             }
168         }
169         
170     double ans = 0.0;
171     for(int i = 0; i < n; ++i)
172     {
173         if(del[i]) continue;
174         last->clear();
175         Point p1, p2;
176         for(int j = 0; j < n; ++j)
177         {
178             if(del[j] || i == j) continue;
179             if(!cir[i].intersect(cir[j])) return false;
180             cur->clear();
181             IntersectionPoint(cir[i], cir[j], p1, p2);
182             double rs = Angle(p2 - cir[i].c);
183             double rt = Angle(p1 - cir[i].c);
184             if(dcmp(rs) < 0) rs += 2 * PI;
185             if(dcmp(rt) < 0) rt += 2 * PI;
186             if(last->n == 0)
187             {
188                 if(dcmp(rt - rs) < 0)
189                 {
190                     cur->add(Region(rs, 2*PI));
191                     cur->add(Region(0,  rt));
192                 }
193                 else cur->add(Region(rs, rt));
194             }
195             else
196             {
197                 for(int k = 0; k < last->n; ++k)
198                 {
199                     if(dcmp(rt - rs) < 0)
200                     {
201                         if(dcmp(last->v[k].st-rt) >= 0 && dcmp(last->v[k].ed-rs) <= 0) continue;
202                         if(dcmp(last->v[k].st-rt) < 0) cur->add(Region(last->v[k].st, std::min(last->v[k].ed, rt)));
203                         if(dcmp(last->v[k].ed-rs) > 0) cur->add(Region(std::max(last->v[k].st, rs), last->v[k].ed));
204                     }
205                     else
206                     {
207                         if(dcmp(rt-last->v[k].st <= 0 || dcmp(rs-last->v[k].ed) >= 0)) continue;
208                         cur->add(Region(std::max(rs, last->v[k].st), std::min(rt, last->v[k].ed)));
209                     }
210                 }
211             }
212             std::swap(cur, last);
213             if(last->n == 0) break;
214         }
215         for(int j = 0; j < last->n; ++j)
216         {
217             p1 = cir[i].get_point(last->v[j].st);
218             p2 = cir[i].get_point(last->v[j].ed);
219             ans += Cross(p1, p2) / 2;
220             double ang = last->v[j].ed - last->v[j].st;
221             ans += cir[i].r * cir[i].r * (ang - sin(ang)) / 2;
222         }
223     }
224     
225     if(dcmp(ans) == 0) return false;
226     printf("The total possible area is %.2f.\n", ans);
227     return true;
228 }
229 
230 int main(void)
231 {
232     //freopen("3467in.txt", "r", stdin);
233     last = new Region_vector;
234     cur =  new Region_vector;
235     while(scanf("%lf", &r) == 1)
236     {
237         Point t;
238         for(int i = 0; i < n; ++i)
239         {
240             t.read();
241             cir[i] = Circle(t, r);
242         }
243         if(!solve())
244             puts("Poor iSea, maybe 2012 is coming!");
245     }
246     
247     return 0;
248 }