hdu 5134 Highway

题意是：你当前距离高速公路D米，你在高速公路的行驶速度为v1,在其他地方行驶速度是v0  (v1>v0)，然后问你在T秒能到达的区域的面积

我们可以在高速公路上找到一点p，只考虑上半个平面（下半个是一样的）,使得 从 now 到 p，然后在从 p 沿着高速公路向上走到达一个最远距离，可以用三分得到这个点

而且此时 now -> p的时间会小于等于 now -> O + O ->p 证明如下：

因为肯定存在一点ST,使得 now->ST 的时间等于 now->O + O->ST

假设 p点在ST的上面，now->p 的时间会大于 now-ST + ST->p 的时间，所以走now - >ST +ST->p 的路径会更优，所以 p点一定在 O到ST 这段里面

所以在从now 到  O到p  这一段中的点 p1，直接从now->p1所花费的时间会最少，然后在以p1为出发点所能到达的区域就是以 p1为圆心的圆，并且这个圆会与大圆相内切

对于p点上面部分的点，以它们为圆心的圆 会形成一个三角形（证明略）

所以我们的答案就是 ans = 大圆面积 + 2*三角形面积 - 2*三角形与圆的面积交

  1 //#include<bits/stdc++.h>
  2 #include<iostream>
  3 #include<cstdio>
  4 #include<cmath>
  5 #include<cstring>
  6 #include<vector>
  7 #include<algorithm>
  8 using namespace std;
  9 typedef long long LL;
 10
 11 const double PI = acos(-1.0);
 12 const double EPS = 1e-7;
 13
 14 inline int sgn(double x) {
 15     return (x > EPS) - (x < -EPS);
 16 }
 17
 18 struct Point {
 19     double x, y;
 20     Point() {}
 21     Point(double x, double y): x(x), y(y) {}
 22     void read() {
 23         scanf("%lf%lf", &x, &y);
 24     }
 25     double angle() {
 26         return atan2(y, x);
 27     }
 28     Point operator + (const Point &rhs) const {
 29         return Point(x + rhs.x, y + rhs.y);
 30     }
 31     Point operator - (const Point &rhs) const {
 32         return Point(x - rhs.x, y - rhs.y);
 33     }
 34     Point operator * (double t) const {
 35         return Point(x * t, y * t);
 36     }
 37     Point operator / (double t) const {
 38         return Point(x / t, y / t);
 39     }
 40     double operator *(const Point &b)const
 41     {
 42         return x*b.x + y*b.y;
 43     }
 44     double length() const {
 45         return sqrt(x * x + y * y);
 46     }
 47     Point unit() const {            //单位向量
 48         double l = length();
 49         return Point(x / l, y / l);
 50     }
 51 };
 52 double cross(Point a,Point b) {
 53     return a.x * b.y - a.y * b.x;
 54 }
 55 double cross(Point a,Point b,Point c){
 56     return (a.x-c.x)*(b.y-c.y)-(a.y-c.y)*(b.x-c.x);
 57 }
 58 double xml(Point a,Point b,Point c){
 59     return (a.x-c.x)*(b.x-c.x)+(a.y-c.y)*(b.y-c.y);
 60 }
 61 double sqr(double x) {
 62     return x * x;
 63 }
 64 double dist(const Point &p1, const Point &p2) {
 65     return (p1 - p2).length();
 66 }
 67 double sdist(Point a,Point b){
 68     return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y);
 69 }
 70 //向量 op 逆时针旋转 angle
 71 Point rotate(const Point &p, double angle, const Point &o = Point(0, 0)) {
 72     Point t = p - o;
 73     double x = t.x * cos(angle) - t.y * sin(angle);
 74     double y = t.y * cos(angle) + t.x * sin(angle);
 75     return Point(x, y) + o;
 76 }
 77 Point line_inter(Point A,Point B,Point C,Point D){ //直线相交交点
 78         Point ans;
 79         double a1=A.y-B.y;
 80         double b1=B.x-A.x;
 81         double c1=A.x*B.y-B.x*A.y;
 82
 83         double a2=C.y-D.y;
 84         double b2=D.x-C.x;
 85         double c2=C.x*D.y-D.x*C.y;
 86
 87         ans.x=(b1*c2-b2*c1)/(a1*b2-a2*b1);
 88         ans.y=(a2*c1-a1*c2)/(a1*b2-a2*b1);
 89         return ans;
 90 }
 91 Point p_to_seg(Point p,Point a,Point b){        //点到线段的最近点
 92     Point tmp=p;
 93     tmp.x+=a.y-b.y;
 94     tmp.y+=b.x-a.x;
 95     if(cross(a-p,tmp-p)*cross(b-p,tmp-p)>0) return dist(p,a)<dist(p,b)?a:b;
 96     return line_inter(p,tmp,a,b);
 97 }
 98 void line_circle(Point c,double r,Point a,Point b,Point &p1,Point &p2){
 99     Point tmp=c;
100     double t;
101     tmp.x+=(a.y-b.y);//求垂直于ab的直线
102     tmp.y+=(b.x-a.x);
103     tmp=line_inter(tmp,c,a,b);
104     t=sqrt(sqr(r)-sqr( dist(c,tmp)))/dist(a,b); //比例
105     p1.x=tmp.x+(b.x-a.x)*t;
106     p1.y=tmp.y+(b.y-a.y)*t;
107     p2.x=tmp.x-(b.x-a.x)*t;
108     p2.y=tmp.y-(b.y-a.y)*t;
109 }
110 struct Region {
111     double st, ed;
112     Region() {}
113     Region(double st, double ed): st(st), ed(ed) {}
114     bool operator < (const Region &rhs) const {
115         if(sgn(st - rhs.st)) return st < rhs.st;
116         return ed < rhs.ed;
117     }
118 };
119 struct Circle {
120     Point c;
121     double r;
122     vector<Region> reg;
123     Circle() {}
124     Circle(Point c, double r): c(c), r(r) {}
125     void read() {
126         c.read();
127         scanf("%lf", &r);
128     }
129     void add(const Region &r) {
130         reg.push_back(r);
131     }
132     bool contain(const Circle &cir) const {
133         return sgn(dist(cir.c, c) + cir.r - r) <= 0;
134     }
135     bool intersect(const Circle &cir) const {
136         return sgn(dist(cir.c, c) - cir.r - r) < 0;
137     }
138 };
139 void intersection(const Circle &cir1, const Circle &cir2, Point &p1, Point &p2) {   //两圆相交 交点
140     double l = dist(cir1.c, cir2.c);                            //两圆心的距离
141     double d = (sqr(l) - sqr(cir2.r) + sqr(cir1.r)) / (2 * l);  //cir1圆心到交点直线的距离
142     double d2 = sqrt(sqr(cir1.r) - sqr(d));                     //交点到 两圆心所在直线的距离
143     Point mid = cir1.c + (cir2.c - cir1.c).unit() * d;
144     Point v = rotate(cir2.c - cir1.c, PI / 2).unit() * d2;
145     p1 = mid + v, p2 = mid - v;
146 }
147 Point calc(const Circle &cir, double angle) {
148     Point p = Point(cir.c.x + cir.r, cir.c.y);
149     return rotate(p, angle, cir.c);
150 }
151 const int MAXN = 1010;
152 Circle cir[MAXN],cir2[MAXN];
153 bool del[MAXN];
154 int n;
155 double get_area(Circle* cir,int n) {            //多个圆的相交面积
156     double ans = 0;
157     memset(del,0,sizeof(del));
158     for(int i = 0; i < n; ++i) {
159         for(int j = 0; j < n; ++j) if(!del[j]) {                //删除被包含的圆
160             if(i == j) continue;
161             if(cir[j].contain(cir[i])) {
162                 del[i] = true;
163                 break;
164             }
165         }
166     }
167     for(int i = 0; i < n; ++i) if(!del[i]) {
168         Circle &mc = cir[i];
169         Point p1, p2;
170         bool flag = false;
171         for(int j = 0; j < n; ++j) if(!del[j]) {
172             if(i == j) continue;
173             if(!mc.intersect(cir[j])) continue;
174             flag = true;
175             intersection(mc, cir[j], p1, p2);                   //求出两圆的交点
176             double rs = (p2 - mc.c).angle(), rt = (p1 - mc.c).angle();
177             if(sgn(rs) < 0) rs += 2 * PI;
178             if(sgn(rt) < 0) rt += 2 * PI;
179             if(sgn(rs - rt) > 0) mc.add(Region(rs, PI * 2)), mc.add(Region(0, rt)); //添加相交区域
180             else mc.add(Region(rs, rt));
181         }
182         if(!flag) {
183             ans += PI * sqr(mc.r);
184             continue;
185         }
186         sort(mc.reg.begin(), mc.reg.end());                 //对相交区域进行排序
187         int cnt = 1;
188         for(int j = 1; j < int(mc.reg.size()); ++j) {
189             if(sgn(mc.reg[cnt - 1].ed - mc.reg[j].st) >= 0) {   //如果有区域可以合并，则合并
190                 mc.reg[cnt - 1].ed = max(mc.reg[cnt - 1].ed, mc.reg[j].ed);
191             } else mc.reg[cnt++] = mc.reg[j];
192         }
193         mc.add(Region());
194         mc.reg[cnt] = mc.reg[0];
195         for(int j = 0; j < cnt; ++j) {
196             p1 = calc(mc, mc.reg[j].ed);
197             p2 = calc(mc, mc.reg[j + 1].st);
198             ans += cross(p1, p2) / 2;                           //
199             double angle = mc.reg[j + 1].st - mc.reg[j].ed;
200             if(sgn(angle) < 0) angle += 2 * PI;
201             ans += 0.5 * sqr(mc.r) * (angle - sin(angle));      //弧所对应的的面积
202         }
203     }
204     return ans;
205 }
206 double two_cir(Circle t1,Circle t2){            //两个圆的相交面积
207     if(t1.contain(t2)||t2.contain(t1))    return PI * sqr(min(t2.r,t1.r));
208     if(!t1.intersect(t2)) return 0;
209     double ans=0,len=dist(t1.c,t2.c);
210     double x=(sqr(t1.r)+sqr(len)-sqr(t2.r))/(2*len);
211     double angle1=acos(x/t1.r),angle2=acos((len-x)/t2.r);
212     ans=sqr(t1.r)*angle1+sqr(t2.r)*angle2-len*t1.r*sin(angle1);    // 两个扇形 减去一个四边形面积
213     return ans;
214 }
215 double EP=0;
216 double triangle_circle(Point a,Point b,Point c,double r){//三角形与圆交
217     double A,B,C,x,y,tS;
218     A=dist(b,c);
219     B=dist(a,c);
220     C=dist(b,a);
221     if(A<r&&B<r)
222     return cross(a,b,c)/2;
223     else if(A<r&&B>=r){
224         x=(xml(a,c,b)+sqrt(r*r*C*C-cross(a,c,b)*cross(a,c,b)))/C;
225         tS=cross(a,b,c)/2;
226         return asin(tS*(1-x/C)*2/r/B*(1-EP))*r*r/2+tS*x/C;
227     }
228     else if(A>=r&&B<r){
229         y=(xml(b,c,a)+sqrt(r*r*C*C-cross(b,c,a)*cross(b,c,a)))/C;
230         tS=cross(a,b,c)/2;
231         return asin(tS*(1-y/C)*2/r/A*(1-EP))*r*r/2+tS*y/C;
232     }
233     else if(fabs(cross(a,b,c))>=r*C||xml(b,c,a)<=0||xml(a,c,b)<=0){
234         if(xml(a,b,c)<0)
235             if(cross(a,b,c)<0)
236                 return (-acos(-1.0)-asin(cross(a,b,c)/A/B*(1-EP)))*r*r/2;
237             else return (acos(-1.0)-asin(cross(a,b,c)/A/B*(1-EP)))*r*r/2;
238         else return asin(cross(a,b,c)/A/B*(1-EP))*r*r/2;
239     }
240     else{
241         x=(xml(a,c,b)+sqrt(r*r*C*C-cross(a,c,b)*cross(a,c,b)))/C;
242         y=(xml(b,c,a)+sqrt(r*r*C*C-cross(b,c,a)*cross(b,c,a)))/C;
243         tS=cross(a,b,c)/2;
244         return (asin(tS*(1-x/C)*2/r/B*(1-EP))+asin(tS*(1-y/C)*2/r/A*(1-EP)))*r*r/2+tS*((y+x)/C-1);
245     }
246 }
247 Point pt1[5100],cter;
248 double r;
249 Point three_P_mincover(Point a,Point b,Point c){ //三角形的外接圆
250     Point ret;
251     double a1=b.x-a.x,b1=b.y-a.y,c1=(a1*a1+b1*b1)/2;
252     double a2=c.x-a.x,b2=c.y-a.y,c2=(a2*a2+b2*b2)/2;
253     double d=a1*b2-a2*b1;
254     ret.x=a.x+(c1*b2-c2*b1)/d;
255     ret.y=a.y+(a1*c2-a2*c1)/d;
256     return ret;
257 }
258 void min_circle_cover(){                //点的最小覆盖
259     random_shuffle(pt1,pt1+n);
260      cter=pt1[0];
261      r=0;
262      for(int i=1;i<n;i++){
263         if(dist(cter,pt1[i])-r>EPS){
264             cter=pt1[i];
265             r=0;
266             for(int j=0;j<i;j++){
267                 if(dist(cter,pt1[j])-r>EPS){
268                     cter=(pt1[i]+pt1[j])/2.0;
269                     r=dist(cter,pt1[i]);
270                     for(int k=0;k<j;k++){
271                         if(dist(cter,pt1[k])-r>EPS){
272                             cter=three_P_mincover(pt1[i],pt1[j],pt1[k]);
273                             r=dist(cter,pt1[i]);
274                         }
275                     }
276                 }
277             }
278         }
279      }
280 }
281 int main(){
282     #ifndef ONLINE_JUDGE
283     freopen("input.txt","r",stdin);
284     #endif // ONLINE_JUDGE
285     double v0,v1,D,T;
286     int cas=0;
287     Point cter;
288     while(scanf("%lf%lf%lf%lf",&v0,&v1,&D,&T)!=EOF){
289         double p_t,l=D/v0,r=T,mid,mmid;
290         cter=Point(-D,0);
291         if(D>=v0*T){
292             printf("Case %d: %.8f\n",++cas,PI*v0*v0*T*T);
293             continue;
294         }
295         for(int i=0;i<100;i++){
296             mid=(l+r)/2;
297             mmid=(r+mid)/2;
298             double len1,len2;
299             len1=sqrt(v0*mid*v0*mid-D*D)+v1*(T-mid);
300             len2=sqrt(v0*mmid*v0*mmid-D*D)+v1*(T-mmid);
301             if(len1>len2) r=mmid;
302             else l=mid;
303         }
304         p_t=(l+r)/2;
305         double x1,y1,y;
306         x1=T/p_t*D-D;
307         y1=sqrt(v0*T*v0*T-(x1+D)*(x1+D));
308         y=sqrt(v0*p_t*v0*p_t-D*D)+v1*(T-p_t);
309         double ans=0,tmp=0;
310         Point a,b,c;
311         a=Point(x1,y1);
312         b=Point(0,y);
313         c=Point(-x1,y1);
314         tmp+=triangle_circle(a,b,cter,v0*T);
315         tmp+=triangle_circle(b,c,cter,v0*T);
316         tmp+=triangle_circle(c,a,cter,v0*T);
317         tmp=fabs(tmp);
318         ans=PI*v0*v0*T*T+cross(a-c,b-c);
319         ans-=tmp*2;
320         printf("Case %d: %.8f\n",++cas,ans);
321     }
322 }