Skip to content

Instantly share code, notes, and snippets.

@lychees

lychees/ds.cpp

Last active December 21, 2015 04:38
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save lychees/6250870 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save lychees/6250870 to your computer and use it in GitHub Desktop.
Download ZIP
HDU 4637 Rain on your Fat brother
Raw
ds.cpp
#include <functional>
#include <algorithm>
#include <iostream>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <numeric>
#include <cstring>
#include <cassert>
#include <cstdio>
#include <string>
#include <vector>
#include <bitset>
#include <queue>
#include <stack>
#include <cmath>
#include <ctime>
#include <list>
#include <set>
#include <map>
using namespace std;
#define DO(n) for ( int ____n ## __line__ = n; ____n ## __line__ -- ; )
#define ALL(A) A.begin(), A.end()
#define BSC(A, x) (lower_bound(ALL(A), x) - A.begin())
#define CTN(T, x) (T.find(x) != T.end())
#define SZ(A) int(A.size())
#define PB push_back
#define MP(A, B) make_pair(A, B)
#define fi first
#define se second
typedef long long LL;
typedef vector<int> VI;
typedef map<int, int> MII;
typedef pair<int, int> PII;
typedef pair<LL, LL> PLL;
template<class T> inline void RST(T &A){memset(A, 0, sizeof(A));}
template<class T> inline void FLC(T &A, int x){memset(A, x, sizeof(A));}
template<class T> inline void CLR(T &A){A.clear();}
//}
/** Constant List .. **/ //{
const int dx4[] = {-1, 0, 1, 0};
const int dy4[] = {0, 1, 0, -1};
const int dx8[] = {-1, 0, 1, 0 , -1 , -1 , 1 , 1};
const int dy8[] = {0, 1, 0, -1 , -1 , 1 , -1 , 1};
const int dxhorse[] = {-2 , -2 , -1 , -1 , 1 , 1 , 2 , 2};
const int dyhorse[] = {1 , -1 , 2 , -2 , 2 ,-2 , 1 ,-1};
const int MOD = 1000000007;
//int MOD = 99990001;
const int INF = 0x3f3f3f3f;
const LL INFF = 1LL << 60;
const double EPS = 1e-9;
const double OO = 1e15;
const double PI = acos(-1.0); //M_PI;
//}
template<class T> inline void checkMin(T &a,const T b){if (b<a) a=b;}
template<class T> inline void checkMax(T &a,const T b){if (a<b) a=b;}
//}
template<class T> inline T low_bit(T x) {return x & -x;}
/*****************************************************************/
#define mp make_pair
#define pb push_back
using namespace std;
namespace Geometry{
#define typec double
const typec eps=1e-12;
const typec pi=acos(-1.0);
const typec inf=1e20;
int dblcmp(double d){
if (fabs(d)<eps)return 0;
return d>eps?1:-1;
}
int sgn(double a) {return a<-eps?-1:a>eps;}
inline double sqr(double x){return x*x;}
struct Point2D{
typec x,y;
Point2D(){}
Point2D(typec _x,typec _y):x(_x),y(_y){};
void input(){
scanf("%lf%lf",&x,&y);
}
void output(){
printf("%.2f %.2f\n",x,y);
}
bool operator==(Point2D a)const{
return dblcmp(a.x-x)==0&&dblcmp(a.y-y)==0;
}
bool operator<(Point2D a)const{
return dblcmp(a.x-x)==0?dblcmp(y-a.y)<0:x<a.x;
}
typec len(){
return hypot(x,y);
}
typec len2(){
return x*x+y*y;
}
Point2D operator + (const Point2D &A) const{
return Point2D(x + A.x , y + A.y);
}
Point2D operator - (const Point2D &A) const{
return Point2D(x - A.x , y - A.y);
}
Point2D operator * (const typec _x) const{
return Point2D(x * _x , y * _x);
}
typec operator * (const Point2D &A) const{
return x * A.x + y * A.y;
}
typec operator ^ (const Point2D &A) const{
return x * A.y - y * A.x;
}
Point2D operator / (const typec _p) const{
return Point2D(x / _p , y / _p);
}
typec distance(Point2D p){
return hypot(x-p.x,y-p.y);
}
Point2D add(Point2D p){
return Point2D(x+p.x,y+p.y);
}
Point2D sub(Point2D p){
return Point2D(x-p.x,y-p.y);
}
Point2D mul(typec b){
return Point2D(x*b,y*b);
}
Point2D div(typec b){
return Point2D(x/b,y/b);
}
typec dot(Point2D p){
return x*p.x+y*p.y;
}
typec det(Point2D p){
return x*p.y-y*p.x;
}
typec rad(Point2D a,Point2D b){
Point2D p=*this;
return fabs(atan2(fabs(a.sub(p).det(b.sub(p))),a.sub(p).dot(b.sub(p))));
}
Point2D trunc(typec r){
typec l=len();
if (!dblcmp(l))return *this;
r/=l;
return Point2D(x*r,y*r);
}
Point2D rotleft(){
return Point2D(-y,x);
}
Point2D rotright(){
return Point2D(y,-x);
}
Point2D rotate(Point2D p,typec angle)//绕点p逆时针旋转angle角度
{
Point2D v=this->sub(p);
typec c=cos(angle),s=sin(angle);
return Point2D(p.x+v.x*c-v.y*s,p.y+v.x*s+v.y*c);
}
};
typec cross(Point2D a,Point2D b,Point2D c){
return (b.sub(a)).det(c.sub(a));
}
struct Line2D{
Point2D a,b;
double at;
Line2D(){}
Line2D(Point2D _a,Point2D _b){
a=_a;
b=_b;
}
bool operator==(Line2D v){
return (a==v.a)&&(b==v.b);
}
bool operator < (const Line2D &t) const {
return sgn(at-t.at)<0||
(sgn(at-t.at)==0&&sgn(cross(a,b,t.a))<0);
}
Point2D operator &(Line2D &A)const
{
Point2D res = a;
//注意：有些题目可能会有直线相交或者重合情况
//可以把返回值改成`pair<Point,int>`来返回两直线的状态。
double t = ((a - A.a) ^ (A.a - A.b)) / ((a - b) ^ (A.a - A.b));
res.x += (b.x - a.x) * t;
res.y += (b.y - a.y) * t;
return res;
}
//倾斜角angle
Line2D(Point2D p,double angle)
{
a=p;
if (dblcmp(angle-pi/2)==0)
{
b=a.add(Point2D(0,1));
}
else
{
b=a.add(Point2D(1,tan(angle)));
}
}
//矢量V以P为顶点逆时针旋转angle并放大scale倍
Point2D rotate(Point2D v,Point2D p,double angle,double scale){
Point2D ret=p;
v.x-=p.x,v.y-=p.y;
p.x=scale*cos(angle);
p.y=scale*sin(angle);
ret.x+=v.x*p.x-v.y*p.y;
ret.y+=v.x*p.y+v.y*p.x;
return ret;
}
//ax+by+c=0
Line2D(typec _a,typec _b,typec _c){
if (dblcmp(_a)==0)
{
a=Point2D(0,-_c/_b);
b=Point2D(1,-_c/_b);
}
else if (dblcmp(_b)==0)
{
a=Point2D(-_c/_a,0);
b=Point2D(-_c/_a,1);
}
else
{
a=Point2D(0,-_c/_b);
b=Point2D(1,(-_c-_a)/_b);
}
}
void input()
{
a.input();
b.input();
angle();
}
void adjust()
{
if (b<a)swap(a,b);
}
double length()
{
return a.distance(b);
}
double angle()//直线倾斜角 0<=angle<180
{
double k=atan2(b.y-a.y,b.x-a.x);
if (dblcmp(k)<0)k+=pi;
if (dblcmp(k-pi)==0)k-=pi;
at = k;
return k;
}
//点和线段关系
//1 在逆时针
//2 在顺时针
//3 平行
int relation(Point2D p)
{
int c=dblcmp(p.sub(a).det(b.sub(a)));
if (c<0)return 1;
if (c>0)return 2;
return 3;
}
bool pointonseg(Point2D p)
{
return dblcmp(p.sub(a).det(b.sub(a)))==0&&dblcmp(p.sub(a).dot(p.sub(b)))<=0;
}
bool parallel(Line2D v)
{
return dblcmp(b.sub(a).det(v.b.sub(v.a)))==0;
}
//2 规范相交
//1 非规范相交
//0 不相交
int segcrossseg(Line2D v)
{
int d1=dblcmp(b.sub(a).det(v.a.sub(a)));
int d2=dblcmp(b.sub(a).det(v.b.sub(a)));
int d3=dblcmp(v.b.sub(v.a).det(a.sub(v.a)));
int d4=dblcmp(v.b.sub(v.a).det(b.sub(v.a)));
if ((d1^d2)==-2&&(d3^d4)==-2)return 2;
return (d1==0&&dblcmp(v.a.sub(a).dot(v.a.sub(b)))<=0||
d2==0&&dblcmp(v.b.sub(a).dot(v.b.sub(b)))<=0||
d3==0&&dblcmp(a.sub(v.a).dot(a.sub(v.b)))<=0||
d4==0&&dblcmp(b.sub(v.a).dot(b.sub(v.b)))<=0);
}
int linecrossseg(Line2D v)//*this seg v line
{
int d1=dblcmp(b.sub(a).det(v.a.sub(a)));
int d2=dblcmp(b.sub(a).det(v.b.sub(a)));
if ((d1^d2)==-2)return 2;
return (d1==0||d2==0);
}
//0 平行
//1 重合
//2 相交
int linecrossline(Line2D v)
{
if ((*this).parallel(v))
{
return v.relation(a)==3;
}
return 2;
}
Point2D crosspoint(Line2D v)
{
double a1=v.b.sub(v.a).det(a.sub(v.a));
double a2=v.b.sub(v.a).det(b.sub(v.a));
return Point2D((a.x*a2-b.x*a1)/(a2-a1),(a.y*a2-b.y*a1)/(a2-a1));
}
double dispointtoline(Point2D p)
{
return fabs(p.sub(a).det(b.sub(a)))/length();
}
double dispointtoseg(Point2D p)
{
if (dblcmp(p.sub(b).dot(a.sub(b)))<0||dblcmp(p.sub(a).dot(b.sub(a)))<0)
{
return min(p.distance(a),p.distance(b));
}
return dispointtoline(p);
}
Point2D lineprog(Point2D p)
{
return a.add(b.sub(a).mul(b.sub(a).dot(p.sub(a))/b.sub(a).len2()));
}
Point2D symmetrypoint(Point2D p) //对称点
{
Point2D q=lineprog(p);
return Point2D(2*q.x-p.x,2*q.y-p.y);
}
};
struct Circle{
Point2D p;
typec r;
Circle(){}
Circle(Point2D _p,typec _r):
p(_p),r(_r){};
Circle(typec x,typec y,typec _r):p(Point2D(x,y)),r(_r){};
Circle(Point2D a,Point2D b,Point2D c){//三角形的外接圆
p=Line2D(a.add(b).div(2),a.add(b).div(2).add(b.sub(a).rotleft())).crosspoint(Line2D(c.add(b).div(2),c.add(b).div(2).add(b.sub(c).rotleft())));
r=p.distance(a);
}
Circle(Point2D a,Point2D b,Point2D c,bool t)//三角形的内切圆
{
Line2D u,v;
double m=atan2(b.y-a.y,b.x-a.x),n=atan2(c.y-a.y,c.x-a.x);
u.a=a;
u.b=u.a.add(Point2D(cos((n+m)/2),sin((n+m)/2)));
v.a=b;
m=atan2(a.y-b.y,a.x-b.x),n=atan2(c.y-b.y,c.x-b.x);
v.b=v.a.add(Point2D(cos((n+m)/2),sin((n+m)/2)));
p=u.crosspoint(v);
r=Line2D(a,b).dispointtoseg(p);
}
void input()
{
p.input();
scanf("%lf",&r);
}
void output()
{
printf("%.2lf %.2lf %.2lf\n",p.x,p.y,r);
}
bool operator==(Circle v)
{
return ((p==v.p)&&dblcmp(r-v.r)==0);
}
bool operator<(Circle v)const
{
return ((p<v.p)||(p==v.p)&&dblcmp(r-v.r)<0);
}
double area()
{
return pi*sqr(r);
}
double circumference()
{
return 2*pi*r;
}
//0 圆外
//1 圆上
//2 圆内
int relation(Point2D b)
{
double dst=b.distance(p);
if (dblcmp(dst-r)<0)return 2;
if (dblcmp(dst-r)==0)return 1;
return 0;
}
int relationseg(Line2D v)
{
double dst=v.dispointtoseg(p);
if (dblcmp(dst-r)<0)return 2;
if (dblcmp(dst-r)==0)return 1;
return 0;
}
int relationline(Line2D v)
{
double dst=v.dispointtoline(p);
if (dblcmp(dst-r)<0)return 2;
if (dblcmp(dst-r)==0)return 1;
return 0;
}
//过a b两点 半径r的两个圆
int getCircle(Point2D a,Point2D b,double r,Circle&c1,Circle&c2)
{
Circle x(a,r),y(b,r);
int t=x.pointcrossCircle(y,c1.p,c2.p);
if (!t)return 0;
c1.r=c2.r=r;
return t;
}
//与直线u相切 过点q 半径r1的圆
int getCircle(Line2D u,Point2D q,double r1,Circle &c1,Circle &c2)
{
double dis=u.dispointtoline(q);
if (dblcmp(dis-r1*2)>0)return 0;
if (dblcmp(dis)==0)
{
c1.p=q.add(u.b.sub(u.a).rotleft().trunc(r1));
c2.p=q.add(u.b.sub(u.a).rotright().trunc(r1));
c1.r=c2.r=r1;
return 2;
}
Line2D u1=Line2D(u.a.add(u.b.sub(u.a).rotleft().trunc(r1)),u.b.add(u.b.sub(u.a).rotleft().trunc(r1)));
Line2D u2=Line2D(u.a.add(u.b.sub(u.a).rotright().trunc(r1)),u.b.add(u.b.sub(u.a).rotright().trunc(r1)));
Circle cc=Circle(q,r1);
Point2D p1,p2;
if (!cc.pointcrossline(u1,p1,p2))cc.pointcrossline(u2,p1,p2);
c1=Circle(p1,r1);
if (p1==p2)
{
c2=c1;return 1;
}
c2=Circle(p2,r1);
return 2;
}
//同时与直线u,v相切 半径r1的圆
int getCircle(Line2D u,Line2D v,typec r1,Circle &c1,Circle &c2,Circle &c3,Circle &c4)
{
if (u.parallel(v))return 0;
Line2D u1=Line2D(u.a.add(u.b.sub(u.a).rotleft().trunc(r1)),u.b.add(u.b.sub(u.a).rotleft().trunc(r1)));
Line2D u2=Line2D(u.a.add(u.b.sub(u.a).rotright().trunc(r1)),u.b.add(u.b.sub(u.a).rotright().trunc(r1)));
Line2D v1=Line2D(v.a.add(v.b.sub(v.a).rotleft().trunc(r1)),v.b.add(v.b.sub(v.a).rotleft().trunc(r1)));
Line2D v2=Line2D(v.a.add(v.b.sub(v.a).rotright().trunc(r1)),v.b.add(v.b.sub(v.a).rotright().trunc(r1)));
c1.r=c2.r=c3.r=c4.r=r1;
c1.p=u1.crosspoint(v1);
c2.p=u1.crosspoint(v2);
c3.p=u2.crosspoint(v1);
c4.p=u2.crosspoint(v2);
return 4;
}
//同时与不相交圆cx,cy相切 半径为r1的圆
int getCircle(Circle cx,Circle cy,double r1,Circle&c1,Circle&c2)
{
Circle x(cx.p,r1+cx.r),y(cy.p,r1+cy.r);
int t=x.pointcrossCircle(y,c1.p,c2.p);
if (!t)return 0;
c1.r=c2.r=r1;
return t;
}
int pointcrossline(Line2D v,Point2D &p1,Point2D &p2)//求与线段交要先判断relationseg
{
if (!(*this).relationline(v))return 0;
Point2D a=v.lineprog(p);
double d=v.dispointtoline(p);
d=sqrt(r*r-d*d);
if (dblcmp(d)==0)
{
p1=a;
p2=a;
return 1;
}
p1=a.sub(v.b.sub(v.a).trunc(d));
p2=a.add(v.b.sub(v.a).trunc(d));
return 2;
}
//5 相离
//4 外切
//3 相交
//2 内切
//1 内含
int relationCircle(Circle v)
{
double d=p.distance(v.p);
if (dblcmp(d-r-v.r)>0)return 5;
if (dblcmp(d-r-v.r)==0)return 4;
double l=fabs(r-v.r);
if (dblcmp(d-r-v.r)<0&&dblcmp(d-l)>0)return 3;
if (dblcmp(d-l)==0)return 2;
if (dblcmp(d-l)<0)return 1;
}
int pointcrossCircle(Circle v,Point2D &p1,Point2D &p2)
{
int rel=relationCircle(v);
if (rel==1||rel==5)return 0;
double d=p.distance(v.p);
double l=(d+(sqr(r)-sqr(v.r))/d)/2;
double h=sqrt(sqr(r)-sqr(l));
p1=p.add(v.p.sub(p).trunc(l).add(v.p.sub(p).rotleft().trunc(h)));
p2=p.add(v.p.sub(p).trunc(l).add(v.p.sub(p).rotright().trunc(h)));
if (rel==2||rel==4)
{
return 1;
}
return 2;
}
//过一点做圆的切线 (先判断点和圆关系)
int tangentline(Point2D q,Line2D &u,Line2D &v){
int x=relation(q);
if (x==2)return 0;
if (x==1)
{
u=Line2D(q,q.add(q.sub(p).rotleft()));
v=u;
return 1;
}
double d=p.distance(q);
double l=sqr(r)/d;
double h=sqrt(sqr(r)-sqr(l));
u=Line2D(q,p.add(q.sub(p).trunc(l).add(q.sub(p).rotleft().trunc(h))));
v=Line2D(q,p.add(q.sub(p).trunc(l).add(q.sub(p).rotright().trunc(h))));
return 2;
}
double areaCircle(Circle v){
int rel=relationCircle(v);
if (rel>=4)return 0.0;
if (rel<=2)return min(area(),v.area());
double d=p.distance(v.p);
double hf=(r+v.r+d)/2.0;
double ss=2*sqrt(hf*(hf-r)*(hf-v.r)*(hf-d));
double a1=acos((r*r+d*d-v.r*v.r)/(2.0*r*d));
a1=a1*r*r;
double a2=acos((v.r*v.r+d*d-r*r)/(2.0*v.r*d));
a2=a2*v.r*v.r;
return a1+a2-ss;
}
double areatriangle(Point2D a,Point2D b){
if (dblcmp(p.sub(a).det(p.sub(b))==0))return 0.0;
Point2D q[5];
int len=0;
q[len++]=a;
Line2D l(a,b);
Point2D p1,p2;
if (pointcrossline(l,q[1],q[2])==2){
if (dblcmp(a.sub(q[1]).dot(b.sub(q[1])))<0)q[len++]=q[1];
if (dblcmp(a.sub(q[2]).dot(b.sub(q[2])))<0)q[len++]=q[2];
}
q[len++]=b;
if (len==4&&(dblcmp(q[0].sub(q[1]).dot(q[2].sub(q[1])))>0))swap(q[1],q[2]);
double res=0;
int i;
for (i=0;i<len-1;i++){
if (relation(q[i])==0||relation(q[i+1])==0){
double arg=p.rad(q[i],q[i+1]);
res+=r*r*arg/2.0;
}
else{
res+=fabs(q[i].sub(p).det(q[i+1].sub(p))/2.0);
}
}
return res;
}
};
}using namespace Geometry;
double v1 , v2 , v , t , X;
int Case;
const int N = 50009;
Line2D l;
Circle c;
struct Rain{
double x , y , r , h;
double in , out;
bool crash;
void input(){
in = INF;
out = -INF;
scanf("%lf%lf%lf%lf" , &x , &y ,&r , &h);
}
void intime(double s){
checkMax(out , s);
checkMin(in , s);
}
void gao(){
c = Circle(Point2D(x , y) , r);
l = Line2D(Point2D(X , 0) , Point2D(X - v1 * INF , v * INF));
Line2D AB , AC , BC;
AB = Line2D(Point2D(x - r , y) , Point2D(x , y + h));
AC = Line2D(Point2D(x + r , y) , Point2D(x , y + h));
BC = Line2D(Point2D(x - r , y) , Point2D(x + r , y));
if (AB.segcrossseg(l)){
Point2D p = AB.crosspoint(l);
intime((X - p.x) / v1);
// puts("AB");
// p.output();
}
if (AC.segcrossseg(l)){
Point2D p = AC.crosspoint(l);
intime((X - p.x) / v1);
// puts("AC");
// p.output();
}
if (BC.segcrossseg(l)){
Point2D p = BC.crosspoint(l);
intime((X - p.x) / v1);
// puts("BC");
// p.output();
}
Point2D p , q;
if (c.relationline(l) == 2){
c.pointcrossline(l , p , q);
if (dblcmp(y - p.y) >= 0) intime((X - p.x) / v1);
if (dblcmp(y - q.y) >= 0) intime((X - q.x) / v1);
// puts("Circle");
}
}
}R[N];
int n;
struct DB{
double x;
void in(){
scanf("%lf" , &x);
}
void out(){
printf("%.0lf" , x);
}
DB(){}
DB(double _x):x(_x){}
DB operator + (const DB & A) const{
return DB(x + A.x);
}
DB operator + (double A) const{
return DB(x + A);
}
DB operator - (const DB & A) const{
return DB(x + A.x);
}
DB operator - (double A) const{
return DB(x - A);
}
DB operator * (const DB & A) const{
return DB(x * A.x);
}
DB operator * (double A) const{
return DB(x * A);
}
bool operator == (const DB & A) const{
return dblcmp(x - A.x) == 0;
}
bool operator < (const DB & A) const{
return dblcmp(x - A.x) < 0;
}
bool operator < (double A) const{
return dblcmp(x - A) < 0;
}
bool operator > (const DB & A) const{
return dblcmp(x - A.x) > 0;
}
bool operator != (const DB & A) const{
return !(*this == A);
}
bool operator <= (const DB & A) const{
return *this < A || *this == A;
}
};
DB hashh[N];
int cnt[N];
void solve(){
scanf("%lf%lf%lf%lf%lf" , &v1 ,&v2 , &v , &t , &X);
t = t + v1 * t / (v2 - v1);
scanf("%d" , &n);
for(int i = 0 ; i < n ; ++i){
R[i].input();
R[i].gao();
}
int m = 0;
for (int i = 0 ; i < n ; ++i)
if (dblcmp(R[i].out - R[i].in) > 0){
hashh[m++] = R[i].in;
hashh[m++] = R[i].out;
}
hashh[m++] = t;
hashh[m++] = 0;
sort(hashh , hashh + m);
m = (unique(hashh , hashh + m) - hashh);
RST(cnt);
for (int i = 0 ; i < n ; ++i)
if (dblcmp(R[i].out - R[i].in) > 0){
// printf("%.5lf %.5lf\n" , R[i].in , R[i].out);
int Gin = lower_bound(hashh , hashh + m , DB(R[i].in)) - hashh;
int Gout = lower_bound(hashh , hashh + m , DB(R[i].out)) - hashh;
cnt[Gin]++;
cnt[Gout]--;
}
int now = 0;
double ans = 0;
for (int i = 0 ; i < m - 1 && hashh[i + 1] <= DB(t) ; ++i){
// printf(" %.5lf %.5lf" , hashh[i].x , hashh[i + 1].x);
now += cnt[i];
if (now && dblcmp(hashh[i].x) >= 0)
ans += hashh[i + 1].x - hashh[i].x;
}
// puts("");
printf("Case %d: %.4lf\n" , ++Case , ans + eps);
}
int main(){
#ifndef ONLINE_JUDGE
freopen("in.txt", "r", stdin);
//freopen("out.txt", "w", stdout);
#endif
int _;Case = 0;
cin >> _;
while(_--) solve();
}
Raw
ds_geometry.cpp
#include<vector>
#include<list>
#include<map>
#include<set>
#include<deque>
#include<queue>
#include<stack>
#include<bitset>
#include<algorithm>
#include<functional>
#include<numeric>
#include<utility>
#include<iostream>
#include<sstream>
#include<iomanip>
#include<cstdio>
#include<cmath>
#include<cstdlib>
#include<cctype>
#include<string>
#include<cstring>
#include<cstdio>
#include<cmath>
#include<cstdlib>
#include<ctime>
#include<climits>
#include<complex>
#define mp make_pair
#define pb push_back
using namespace std;
namespace Geometry{
#define typec double
const typec eps=1e-8;
const typec pi=acos(-1.0);
const typec inf=1e20;
const int maxp=100001;
int dblcmp(double d){
if (fabs(d)<eps)return 0;
return d>eps?1:-1;
}
int sgn(double a) {return a<-eps?-1:a>eps;}
inline double sqr(double x){return x*x;}
struct Point2D{
typec x,y;
Point2D(){}
Point2D(typec _x,typec _y):x(_x),y(_y){};
void input(){
scanf("%lf%lf",&x,&y);
}
void output(){
printf("%.2f %.2f\n",x,y);
}
bool operator==(Point2D a)const{
return dblcmp(a.x-x)==0&&dblcmp(a.y-y)==0;
}
bool operator<(Point2D a)const{
return dblcmp(a.x-x)==0?dblcmp(y-a.y)<0:x<a.x;
}
typec len(){
return hypot(x,y);
}
typec len2(){
return x*x+y*y;
}
Point2D operator + (const Point2D &A) const{
return Point2D(x + A.x , y + A.y);
}
Point2D operator - (const Point2D &A) const{
return Point2D(x - A.x , y - A.y);
}
Point2D operator * (const typec _x) const{
return Point2D(x * _x , y * _x);
}
typec operator * (const Point2D &A) const{
return x * A.x + y * A.y;
}
typec operator ^ (const Point2D &A) const{
return x * A.y - y * A.x;
}
Point2D operator / (const typec _p) const{
return Point2D(x / _p , y / _p);
}
typec distance(Point2D p){
return hypot(x-p.x,y-p.y);
}
Point2D add(Point2D p){
return Point2D(x+p.x,y+p.y);
}
Point2D sub(Point2D p){
return Point2D(x-p.x,y-p.y);
}
Point2D mul(typec b){
return Point2D(x*b,y*b);
}
Point2D div(typec b){
return Point2D(x/b,y/b);
}
typec dot(Point2D p){
return x*p.x+y*p.y;
}
typec det(Point2D p){
return x*p.y-y*p.x;
}
typec rad(Point2D a,Point2D b){
Point2D p=*this;
return fabs(atan2(fabs(a.sub(p).det(b.sub(p))),a.sub(p).dot(b.sub(p))));
}
Point2D trunc(typec r){
typec l=len();
if (!dblcmp(l))return *this;
r/=l;
return Point2D(x*r,y*r);
}
Point2D rotleft(){
return Point2D(-y,x);
}
Point2D rotright(){
return Point2D(y,-x);
}
Point2D rotate(Point2D p,typec angle)//绕点p逆时针旋转angle角度
{
Point2D v=this->sub(p);
typec c=cos(angle),s=sin(angle);
return Point2D(p.x+v.x*c-v.y*s,p.y+v.x*s+v.y*c);
}
};
typec cross(Point2D a,Point2D b,Point2D c){
return (b.sub(a)).det(c.sub(a));
}
struct Line2D{
Point2D a,b;
double at;
Line2D(){}
Line2D(Point2D _a,Point2D _b){
a=_a;
b=_b;
}
bool operator==(Line2D v){
return (a==v.a)&&(b==v.b);
}
bool operator < (const Line2D &t) const {
return sgn(at-t.at)<0||
(sgn(at-t.at)==0&&sgn(cross(a,b,t.a))<0);
}
Point2D operator &(Line2D &A)const
{
Point2D res = a;
//注意：有些题目可能会有直线相交或者重合情况
//可以把返回值改成`pair<Point,int>`来返回两直线的状态。
double t = ((a - A.a) ^ (A.a - A.b)) / ((a - b) ^ (A.a - A.b));
res.x += (b.x - a.x) * t;
res.y += (b.y - a.y) * t;
return res;
}
//倾斜角angle
Line2D(Point2D p,double angle)
{
a=p;
if (dblcmp(angle-pi/2)==0)
{
b=a.add(Point2D(0,1));
}
else
{
b=a.add(Point2D(1,tan(angle)));
}
}
//矢量V以P为顶点逆时针旋转angle并放大scale倍
Point2D rotate(Point2D v,Point2D p,double angle,double scale){
Point2D ret=p;
v.x-=p.x,v.y-=p.y;
p.x=scale*cos(angle);
p.y=scale*sin(angle);
ret.x+=v.x*p.x-v.y*p.y;
ret.y+=v.x*p.y+v.y*p.x;
return ret;
}
//ax+by+c=0
Line2D(typec _a,typec _b,typec _c){
if (dblcmp(_a)==0)
{
a=Point2D(0,-_c/_b);
b=Point2D(1,-_c/_b);
}
else if (dblcmp(_b)==0)
{
a=Point2D(-_c/_a,0);
b=Point2D(-_c/_a,1);
}
else
{
a=Point2D(0,-_c/_b);
b=Point2D(1,(-_c-_a)/_b);
}
}
void input()
{
a.input();
b.input();
angle();
}
void adjust()
{
if (b<a)swap(a,b);
}
double length()
{
return a.distance(b);
}
double angle()//直线倾斜角 0<=angle<180
{
double k=atan2(b.y-a.y,b.x-a.x);
if (dblcmp(k)<0)k+=pi;
if (dblcmp(k-pi)==0)k-=pi;
at = k;
return k;
}
//点和线段关系
//1 在逆时针
//2 在顺时针
//3 平行
int relation(Point2D p)
{
int c=dblcmp(p.sub(a).det(b.sub(a)));
if (c<0)return 1;
if (c>0)return 2;
return 3;
}
bool pointonseg(Point2D p)
{
return dblcmp(p.sub(a).det(b.sub(a)))==0&&dblcmp(p.sub(a).dot(p.sub(b)))<=0;
}
bool parallel(Line2D v)
{
return dblcmp(b.sub(a).det(v.b.sub(v.a)))==0;
}
//2 规范相交
//1 非规范相交
//0 不相交
int segcrossseg(Line2D v)
{
int d1=dblcmp(b.sub(a).det(v.a.sub(a)));
int d2=dblcmp(b.sub(a).det(v.b.sub(a)));
int d3=dblcmp(v.b.sub(v.a).det(a.sub(v.a)));
int d4=dblcmp(v.b.sub(v.a).det(b.sub(v.a)));
if ((d1^d2)==-2&&(d3^d4)==-2)return 2;
return (d1==0&&dblcmp(v.a.sub(a).dot(v.a.sub(b)))<=0||
d2==0&&dblcmp(v.b.sub(a).dot(v.b.sub(b)))<=0||
d3==0&&dblcmp(a.sub(v.a).dot(a.sub(v.b)))<=0||
d4==0&&dblcmp(b.sub(v.a).dot(b.sub(v.b)))<=0);
}
int linecrossseg(Line2D v)//*this seg v line
{
int d1=dblcmp(b.sub(a).det(v.a.sub(a)));
int d2=dblcmp(b.sub(a).det(v.b.sub(a)));
if ((d1^d2)==-2)return 2;
return (d1==0||d2==0);
}
//0 平行
//1 重合
//2 相交
int linecrossline(Line2D v)
{
if ((*this).parallel(v))
{
return v.relation(a)==3;
}
return 2;
}
Point2D crosspoint(Line2D v)
{
double a1=v.b.sub(v.a).det(a.sub(v.a));
double a2=v.b.sub(v.a).det(b.sub(v.a));
return Point2D((a.x*a2-b.x*a1)/(a2-a1),(a.y*a2-b.y*a1)/(a2-a1));
}
double dispointtoline(Point2D p)
{
return fabs(p.sub(a).det(b.sub(a)))/length();
}
double dispointtoseg(Point2D p)
{
if (dblcmp(p.sub(b).dot(a.sub(b)))<0||dblcmp(p.sub(a).dot(b.sub(a)))<0)
{
return min(p.distance(a),p.distance(b));
}
return dispointtoline(p);
}
Point2D lineprog(Point2D p)
{
return a.add(b.sub(a).mul(b.sub(a).dot(p.sub(a))/b.sub(a).len2()));
}
Point2D symmetrypoint(Point2D p) //对称点
{
Point2D q=lineprog(p);
return Point2D(2*q.x-p.x,2*q.y-p.y);
}
};
struct Circle{
Point2D p;
typec r;
Circle(){}
Circle(Point2D _p,typec _r):
p(_p),r(_r){};
Circle(typec x,typec y,typec _r):p(Point2D(x,y)),r(_r){};
Circle(Point2D a,Point2D b,Point2D c){//三角形的外接圆
p=Line2D(a.add(b).div(2),a.add(b).div(2).add(b.sub(a).rotleft())).crosspoint(Line2D(c.add(b).div(2),c.add(b).div(2).add(b.sub(c).rotleft())));
r=p.distance(a);
}
Circle(Point2D a,Point2D b,Point2D c,bool t)//三角形的内切圆
{
Line2D u,v;
double m=atan2(b.y-a.y,b.x-a.x),n=atan2(c.y-a.y,c.x-a.x);
u.a=a;
u.b=u.a.add(Point2D(cos((n+m)/2),sin((n+m)/2)));
v.a=b;
m=atan2(a.y-b.y,a.x-b.x),n=atan2(c.y-b.y,c.x-b.x);
v.b=v.a.add(Point2D(cos((n+m)/2),sin((n+m)/2)));
p=u.crosspoint(v);
r=Line2D(a,b).dispointtoseg(p);
}
void input()
{
p.input();
scanf("%lf",&r);
}
void output()
{
printf("%.2lf %.2lf %.2lf\n",p.x,p.y,r);
}
bool operator==(Circle v)
{
return ((p==v.p)&&dblcmp(r-v.r)==0);
}
bool operator<(Circle v)const
{
return ((p<v.p)||(p==v.p)&&dblcmp(r-v.r)<0);
}
double area()
{
return pi*sqr(r);
}
double circumference()
{
return 2*pi*r;
}
//0 圆外
//1 圆上
//2 圆内
int relation(Point2D b)
{
double dst=b.distance(p);
if (dblcmp(dst-r)<0)return 2;
if (dblcmp(dst-r)==0)return 1;
return 0;
}
int relationseg(Line2D v)
{
double dst=v.dispointtoseg(p);
if (dblcmp(dst-r)<0)return 2;
if (dblcmp(dst-r)==0)return 1;
return 0;
}
int relationline(Line2D v)
{
double dst=v.dispointtoline(p);
if (dblcmp(dst-r)<0)return 2;
if (dblcmp(dst-r)==0)return 1;
return 0;
}
//过a b两点 半径r的两个圆
int getCircle(Point2D a,Point2D b,double r,Circle&c1,Circle&c2)
{
Circle x(a,r),y(b,r);
int t=x.pointcrossCircle(y,c1.p,c2.p);
if (!t)return 0;
c1.r=c2.r=r;
return t;
}
//与直线u相切 过点q 半径r1的圆
int getCircle(Line2D u,Point2D q,double r1,Circle &c1,Circle &c2)
{
double dis=u.dispointtoline(q);
if (dblcmp(dis-r1*2)>0)return 0;
if (dblcmp(dis)==0)
{
c1.p=q.add(u.b.sub(u.a).rotleft().trunc(r1));
c2.p=q.add(u.b.sub(u.a).rotright().trunc(r1));
c1.r=c2.r=r1;
return 2;
}
Line2D u1=Line2D(u.a.add(u.b.sub(u.a).rotleft().trunc(r1)),u.b.add(u.b.sub(u.a).rotleft().trunc(r1)));
Line2D u2=Line2D(u.a.add(u.b.sub(u.a).rotright().trunc(r1)),u.b.add(u.b.sub(u.a).rotright().trunc(r1)));
Circle cc=Circle(q,r1);
Point2D p1,p2;
if (!cc.pointcrossline(u1,p1,p2))cc.pointcrossline(u2,p1,p2);
c1=Circle(p1,r1);
if (p1==p2)
{
c2=c1;return 1;
}
c2=Circle(p2,r1);
return 2;
}
//同时与直线u,v相切 半径r1的圆
int getCircle(Line2D u,Line2D v,typec r1,Circle &c1,Circle &c2,Circle &c3,Circle &c4)
{
if (u.parallel(v))return 0;
Line2D u1=Line2D(u.a.add(u.b.sub(u.a).rotleft().trunc(r1)),u.b.add(u.b.sub(u.a).rotleft().trunc(r1)));
Line2D u2=Line2D(u.a.add(u.b.sub(u.a).rotright().trunc(r1)),u.b.add(u.b.sub(u.a).rotright().trunc(r1)));
Line2D v1=Line2D(v.a.add(v.b.sub(v.a).rotleft().trunc(r1)),v.b.add(v.b.sub(v.a).rotleft().trunc(r1)));
Line2D v2=Line2D(v.a.add(v.b.sub(v.a).rotright().trunc(r1)),v.b.add(v.b.sub(v.a).rotright().trunc(r1)));
c1.r=c2.r=c3.r=c4.r=r1;
c1.p=u1.crosspoint(v1);
c2.p=u1.crosspoint(v2);
c3.p=u2.crosspoint(v1);
c4.p=u2.crosspoint(v2);
return 4;
}
//同时与不相交圆cx,cy相切 半径为r1的圆
int getCircle(Circle cx,Circle cy,double r1,Circle&c1,Circle&c2)
{
Circle x(cx.p,r1+cx.r),y(cy.p,r1+cy.r);
int t=x.pointcrossCircle(y,c1.p,c2.p);
if (!t)return 0;
c1.r=c2.r=r1;
return t;
}
int pointcrossline(Line2D v,Point2D &p1,Point2D &p2)//求与线段交要先判断relationseg
{
if (!(*this).relationline(v))return 0;
Point2D a=v.lineprog(p);
double d=v.dispointtoline(p);
d=sqrt(r*r-d*d);
if (dblcmp(d)==0)
{
p1=a;
p2=a;
return 1;
}
p1=a.sub(v.b.sub(v.a).trunc(d));
p2=a.add(v.b.sub(v.a).trunc(d));
return 2;
}
//5 相离
//4 外切
//3 相交
//2 内切
//1 内含
int relationCircle(Circle v)
{
double d=p.distance(v.p);
if (dblcmp(d-r-v.r)>0)return 5;
if (dblcmp(d-r-v.r)==0)return 4;
double l=fabs(r-v.r);
if (dblcmp(d-r-v.r)<0&&dblcmp(d-l)>0)return 3;
if (dblcmp(d-l)==0)return 2;
if (dblcmp(d-l)<0)return 1;
}
int pointcrossCircle(Circle v,Point2D &p1,Point2D &p2)
{
int rel=relationCircle(v);
if (rel==1||rel==5)return 0;
double d=p.distance(v.p);
double l=(d+(sqr(r)-sqr(v.r))/d)/2;
double h=sqrt(sqr(r)-sqr(l));
p1=p.add(v.p.sub(p).trunc(l).add(v.p.sub(p).rotleft().trunc(h)));
p2=p.add(v.p.sub(p).trunc(l).add(v.p.sub(p).rotright().trunc(h)));
if (rel==2||rel==4)
{
return 1;
}
return 2;
}
//过一点做圆的切线 (先判断点和圆关系)
int tangentline(Point2D q,Line2D &u,Line2D &v){
int x=relation(q);
if (x==2)return 0;
if (x==1)
{
u=Line2D(q,q.add(q.sub(p).rotleft()));
v=u;
return 1;
}
double d=p.distance(q);
double l=sqr(r)/d;
double h=sqrt(sqr(r)-sqr(l));
u=Line2D(q,p.add(q.sub(p).trunc(l).add(q.sub(p).rotleft().trunc(h))));
v=Line2D(q,p.add(q.sub(p).trunc(l).add(q.sub(p).rotright().trunc(h))));
return 2;
}
double areaCircle(Circle v){
int rel=relationCircle(v);
if (rel>=4)return 0.0;
if (rel<=2)return min(area(),v.area());
double d=p.distance(v.p);
double hf=(r+v.r+d)/2.0;
double ss=2*sqrt(hf*(hf-r)*(hf-v.r)*(hf-d));
double a1=acos((r*r+d*d-v.r*v.r)/(2.0*r*d));
a1=a1*r*r;
double a2=acos((v.r*v.r+d*d-r*r)/(2.0*v.r*d));
a2=a2*v.r*v.r;
return a1+a2-ss;
}
double areatriangle(Point2D a,Point2D b){
if (dblcmp(p.sub(a).det(p.sub(b))==0))return 0.0;
Point2D q[5];
int len=0;
q[len++]=a;
Line2D l(a,b);
Point2D p1,p2;
if (pointcrossline(l,q[1],q[2])==2){
if (dblcmp(a.sub(q[1]).dot(b.sub(q[1])))<0)q[len++]=q[1];
if (dblcmp(a.sub(q[2]).dot(b.sub(q[2])))<0)q[len++]=q[2];
}
q[len++]=b;
if (len==4&&(dblcmp(q[0].sub(q[1]).dot(q[2].sub(q[1])))>0))swap(q[1],q[2]);
double res=0;
int i;
for (i=0;i<len-1;i++){
if (relation(q[i])==0||relation(q[i+1])==0){
double arg=p.rad(q[i],q[i+1]);
res+=r*r*arg/2.0;
}
else{
res+=fabs(q[i].sub(p).det(q[i+1].sub(p))/2.0);
}
}
return res;
}
};
struct Polygon2D{
int n;
Point2D p[maxp];
Line2D l[maxp];
void copy(Polygon2D & A){
A.n=n;
for (int i=0;i<n;++i){
A.p[i]=p[i];
A.l[i]=l[i];
}
}
void input()
{
for (int i=0;i<n;i++)
{
p[i].input();
}
}
void add(Point2D q)
{
p[n++]=q;
}
void getline()
{
for (int i=0;i<n;i++)
{
l[i]=Line2D(p[i],p[(i+1)%n]);
}
}
struct cmp
{
Point2D p;
cmp(const Point2D &p0){p=p0;}
bool operator()(const Point2D &aa,const Point2D &bb)
{
Point2D a=aa,b=bb;
int d=dblcmp(a.sub(p).det(b.sub(p)));
if (d==0)
{
return dblcmp(a.distance(p)-b.distance(p))<0;
}
return d>0;
}
};
void norm()
{
Point2D mi=p[0];
for (int i=1;i<n;i++)mi=min(mi,p[i]);
sort(p,p+n,cmp(mi));
}
void getconvex(Polygon2D &convex)
{
int i,j,k;
sort(p,p+n);
convex.n=n;
for (i=0;i<min(n,2);i++)
{
convex.p[i]=p[i];
}
if (n<=2)return;
int &top=convex.n;
top=1;
for (i=2;i<n;i++)
{
while (top&&dblcmp(convex.p[top].sub(p[i]).det(convex.p[top-1].sub(p[i])))<=0)
top--;
convex.p[++top]=p[i];
}
int temp=top;
convex.p[++top]=p[n-2];
for (i=n-3;i>=0;i--)
{
while (top!=temp&& dblcmp(convex.p[top].sub(p[i]).det(convex.p[top-1].sub(p[i])))<=0)
top--;
convex.p[++top]=p[i];
}
}
bool isconvex()
{
bool s[3];
memset(s,0,sizeof(s));
int i,j,k;
for (i=0;i<n;i++)
{
j=(i+1)%n;
k=(j+1)%n;
s[dblcmp(p[j].sub(p[i]).det(p[k].sub(p[i])))+1]=1;
if (s[0]&&s[2])return 0;
}
return 1;
}
//3 点上
//2 边上
//1 内部
//0 外部
int relationpoint(Point2D q)
{
int i,j;
for (i=0;i<n;i++)
{
if (p[i]==q)return 3;
}
getline();
for (i=0;i<n;i++)
{
if (l[i].pointonseg(q))return 2;
}
int cnt=0;
for (i=0;i<n;i++)
{
j=(i+1)%n;
int k=dblcmp(q.sub(p[j]).det(p[i].sub(p[j])));
int u=dblcmp(p[i].y-q.y);
int v=dblcmp(p[j].y-q.y);
if (k>0&&u<0&&v>=0)cnt++;
if (k<0&&v<0&&u>=0)cnt--;
}
return cnt!=0;
}
//1 在多边形内长度为正
//2 相交或与边平行
//0 无任何交点
int relationline(Line2D u)
{
int i,j,k=0;
getline();
for (i=0;i<n;i++)
{
if (l[i].segcrossseg(u)==2)return 1;
if (l[i].segcrossseg(u)==1)k=1;
}
if (!k)return 0;
vector<Point2D>vp;
for (i=0;i<n;i++)
{
if (l[i].segcrossseg(u))
{
if (l[i].parallel(u))
{
vp.pb(u.a);
vp.pb(u.b);
vp.pb(l[i].a);
vp.pb(l[i].b);
continue;
}
vp.pb(l[i].crosspoint(u));
}
}
sort(vp.begin(),vp.end());
int sz=vp.size();
for (i=0;i<sz-1;i++)
{
Point2D mid=vp[i].add(vp[i+1]).div(2);
if (relationpoint(mid)==1)return 1;
}
return 2;
}
//直线u切割凸多边形左侧
//注意直线方向
void convexcut(Line2D u,Polygon2D &po)
{
int i,j,k;
int &top=po.n;
top=0;
for (i=0;i<n;i++)
{
int d1=dblcmp(p[i].sub(u.a).det(u.b.sub(u.a)));
int d2=dblcmp(p[(i+1)%n].sub(u.a).det(u.b.sub(u.a)));
if (d1>=0)po.p[top++]=p[i];
if (d1*d2<0)po.p[top++]=u.crosspoint(Line2D(p[i],p[(i+1)%n]));
}
}
double getcircumference()
{
double sum=0;
int i;
for (i=0;i<n;i++)
{
sum+=p[i].distance(p[(i+1)%n]);
}
return sum;
}
double getarea()
{
double sum=0;
int i;
for (i=0;i<n;i++)
{
sum+=p[i].det(p[(i+1)%n]);
}
return fabs(sum)/2;
}
bool getdir()//1代表逆时针 0代表顺时针
{
double sum=0;
int i;
for (i=0;i<n;i++)
{
sum+=p[i].det(p[(i+1)%n]);
}
if (dblcmp(sum)>0)return 1;
return 0;
}
Point2D getbarycentre()
{
Point2D ret(0,0);
double area=0;
int i;
for (i=1;i<n-1;i++)
{
double tmp=p[i].sub(p[0]).det(p[i+1].sub(p[0]));
if (dblcmp(tmp)==0)continue;
area+=tmp;
ret.x+=(p[0].x+p[i].x+p[i+1].x)/3*tmp;
ret.y+=(p[0].y+p[i].y+p[i+1].y)/3*tmp;
}
if (dblcmp(area))ret=ret.div(area);
return ret;
}
double areaintersection(Polygon2D po)
{
}
double areaunion(Polygon2D po)
{
return getarea()+po.getarea()-areaintersection(po);
}
void cg()
{
if (getdir())reverse(p,p+n);
}
int pointinpolygon(Point2D q)//点在凸多边形内部的判定
{
if (dblcmp(q.sub(p[0]).det(p[n-1].sub(p[0])))==0)
{
//if (line(p[n-1],p[0]).pointonseg(q))return n-1;
return -1;
}
if (dblcmp(q.sub(p[0]).det(p[1].sub(p[0])))==0)
{
//if (line(p[n-1],p[0]).pointonseg(q))return n-1;
return -1;
}
int low=1,high=n-2,mid;
static Polygon2D c;
while (low<=high)
{
mid=(low+high)>>1;
if (dblcmp(q.sub(p[0]).det(p[mid].sub(p[0])))>=0&&dblcmp(q.sub(p[0]).det(p[mid+1].sub(p[0])))<0)
{
c.p[0]=p[mid];
c.p[1]=p[mid+1];
c.p[2]=p[0];
c.n=3;
if (Line2D(p[mid],p[mid+1]).pointonseg(q))return -1;
if (c.relationpoint(q))return mid;
return -1;
}
if (dblcmp(q.sub(p[0]).det(p[mid].sub(p[0])))>0)
{
low=mid+1;
}
else
{
high=mid-1;
}
}
return -1;
}
/*
\subsection{直线与凸包求交点}
复杂度$O(\log{n})$。\\
需要先预处理几个东西。\\
\begin{lstlisting}[language=c++]
*/
//`二分[la,lb]这段区间那条边与Line2D相交`
int Gao(int la,int lb,Line2D line)
{
if (la > lb)
lb += n;
int l = la,r = lb,mid;
while (l < r)
{
mid = l+r+1>>1;
if (dblcmp((line.b-line.a)^(p[la]-line.a))*dblcmp((line.b-line.a)^(p[mid]-line.a)) >= 0)
l = mid;
else
r = mid-1;
}
return l%n;
}
//`求l与凸包的交点`
//`先调用Gettheta预处理出凸包每条边的斜率，然后处理成升序排列`
double theta[maxp];
void Gettheta()
{
for (int i = 0;i < n;i++)
{
Point2D v = p[(i+1)%n]-p[i];
theta[i] = atan2(v.y,v.x);
}
for (int i = 1;i < n;i++)
if (theta[i-1] > theta[i]+eps)
theta[i] += 2*pi;
}
double lineInterConvexHull(Line2D l , Point2D & pa , Point2D & pb)
{
double tnow;
Point2D v = l.b-l.a;
tnow = atan2(v.y,v.x);
if (dblcmp(tnow-theta[0]) < 0) tnow += 2*pi;
int pl = lower_bound(theta,theta+n,tnow)-theta;
tnow = atan2(-v.y,-v.x);
if (dblcmp(tnow-theta[0]) < 0) tnow += 2*pi;
int pr = lower_bound(theta,theta+n,tnow)-theta;
//`pl和pr是在l方向上距离最远的点对`
pl = pl%n;
pr = pr%n;
if (dblcmp(v^(p[pl]-l.a))*dblcmp(v^(p[pr]-l.a)) >= 0)
return 0.0;
int xa = Gao(pl,pr,l);
int xb = Gao(pr,pl,l);
if (xa > xb) swap(xa,xb);
//`与[xa,xa+1]和[xb,xb+1]`这两条线段相交
if (dblcmp(v^(p[xa+1]-p[xa])) == 0) return 0.0;
if (dblcmp(v^(p[xb+1]-p[xb])) == 0) return 0.0;
pa = Line2D(p[xa],p[xa+1])&l;
pb = Line2D(p[xb],p[xb+1])&l;
//return MP(pa , pb);
//`题目：求直线切凸包得到的两部分的面积`
//double area0 = sum[xb]-sum[xa+1]+(pa*p[xa+1])/2.0+(p[xb]*pb)/2.0+(pb*pa)/2.0;
//double area1 = sum[xa+n]-sum[xb+1]+(pb*p[xb+1])/2.0+(p[xa]*pa)/2.0+(pa*pb)/2.0;
//return min(area0,area1);
}
};
struct Point3D {
typec x,y,z;
Point3D(typec _x=0,typec _y=0,typec _z=0):x(_x),y(_y),z(_z){}
Point3D operator - (const Point3D &p) const {
return Point3D(x-p.x,y-p.y,z-p.z);
}
Point3D operator + (const Point3D &p) const {
return Point3D(x+p.x,y+p.y,z+p.z);
}
Point3D operator * (const Point3D &p) const {
return Point3D(y*p.z-z*p.y,z*p.x-x*p.z,x*p.y-y*p.x);
}
Point3D operator * (const typec &k) const {
return Point3D(x*k,y*k,z*k);
}
typec operator ^ (const Point3D &p) const {
return x*p.x+y*p.y+z*p.z;
}
typec len2(){return x*x+y*y+z*z; }
typec len() {return sqrt(len2()); }
Point3D normalize() {
double l=len();
return Point3D(x/l,y/l,z/l);
}
void input() {scanf("%lf%lf%lf",&x,&y,&z);}
};
struct Plane3D {
typec A,B,C,D;
bool setPlane(Point3D p1,Point3D p2,Point3D p3) {
Point3D norm=(p1-p2)*(p2-p3);
if (norm.len()<eps) return false;
A=norm.x;B=norm.y;C=norm.z;
D=-(norm^p1);
return true;
}
};
struct Line3D {
Point3D s,d;
void SetLine(Point3D a,Point3D b) {
s=a;d=b-a;
}
};
struct Polygon3D{
const static int MAXV = 505;
struct fac
{
int a, b, c;
bool ok;
};
int n;
Point3D P[MAXV];
int cnt;
fac F[MAXV*8];
int to[MAXV][MAXV];
double vlen(Point3D a)
{
return sqrt(a.x*a.x+a.y*a.y+a.z*a.z);
}
double area(Point3D a,Point3D b, Point3D c)
{
return vlen((b-a)*(c-a));
}
double volume(Point3D a, Point3D b, Point3D c, Point3D d)
{
return (b-a)*(c-a)^(d-a);
}
double ptof(Point3D &p, fac &f)
{
Point3D m = P[f.b]-P[f.a], n = P[f.c]-P[f.a], t = p-P[f.a];
return (m * n) ^ t;
}
void deal(int p, int a, int b)
{
int f = to[a][b];
fac add;
if (F[f].ok)
{
if (ptof(P[p], F[f]) > eps)
dfs(p, f);
else
{
add.a = b, add.b = a, add.c = p, add.ok = 1;
to[p][b] = to[a][p] = to[b][a] = cnt;
F[cnt++] = add;
}
}
}
void dfs(int p, int cur)
{
F[cur].ok = 0;
deal(p, F[cur].b, F[cur].a);
deal(p, F[cur].c, F[cur].b);
deal(p, F[cur].a, F[cur].c);
}
bool same(int s, int t)
{
Point3D &a = P[F[s].a], &b = P[F[s].b], &c = P[F[s].c];
return fabs(volume(a, b, c, P[F[t].a])) < eps && fabs(volume(a, b, c,
P[F[t].b])) < eps && fabs(volume(a, b, c, P[F[t].c])) < eps;
}
void construct()
{
cnt = 0;
if (n < 4)
return;
bool sb = 1;
for (int i = 1; i < n; i++)
{
if (vlen(P[0] - P[i]) > eps)
{
swap(P[1], P[i]);
sb = 0;
break;
}
}
if (sb)return;
sb = 1;
for (int i = 2; i < n; i++)
{
if (vlen((P[0] - P[1]) * (P[1] - P[i])) > eps)
{
swap(P[2], P[i]);
sb = 0;
break;
}
}
if (sb)return;
sb = 1;
for (int i = 3; i < n; i++)
{
if (fabs((P[0] - P[1]) * (P[1] - P[2]) ^ (P[0] - P[i])) > eps)
{
swap(P[3], P[i]);
sb = 0;
break;
}
}
if (sb)return;
fac add;
for (int i = 0; i < 4; i++)
{
add.a = (i+1)%4, add.b = (i+2)%4, add.c = (i+3)%4, add.ok = 1;
if (ptof(P[i], add) > 0)
swap(add.b, add.c);
to[add.a][add.b] = to[add.b][add.c] = to[add.c][add.a] = cnt;
F[cnt++] = add;
}
for (int i = 4; i < n; i++)
{
for (int j = 0; j < cnt; j++)
{
if (F[j].ok && ptof(P[i], F[j]) > eps)
{
dfs(i, j);
break;
}
}
}
int tmp = cnt;
cnt = 0;
for (int i = 0; i < tmp; i++)
{
if (F[i].ok)
{
F[cnt++] = F[i];
}
}
}
//表面积
double area()
{
double ret = 0.0;
for (int i = 0; i < cnt; i++)
{
ret += area(P[F[i].a], P[F[i].b], P[F[i].c]);
}
return ret / 2.0;
}
//体积
double volume()
{
Point3D O(0, 0, 0);
double ret = 0.0;
for (int i = 0; i < cnt; i++)
{
ret += volume(O, P[F[i].a], P[F[i].b], P[F[i].c]);
}
return fabs(ret / 6.0);
}
//表面三角形数
int facetCnt_tri()
{
return cnt;
}
//表面多边形数
int facetCnt()
{
int ans = 0;
for (int i = 0; i < cnt; i++)
{
bool nb = 1;
for (int j = 0; j < i; j++)
{
if (same(i, j))
{
nb = 0;
break;
}
}
ans += nb;
}
return ans;
}
Point3D Fc[MAXV*8];
double V[MAXV*8];
Point3D Center()//重心
{
Point3D O(0,0,0);
for (int i = 0; i < cnt; i++)
{
Fc[i].x = (O.x+P[F[i].a].x+P[F[i].b].x+P[F[i].c].x)/4.0;
Fc[i].y = (O.y+P[F[i].a].y+P[F[i].b].y+P[F[i].c].y)/4.0;
Fc[i].z = (O.z+P[F[i].a].z+P[F[i].b].z+P[F[i].c].z)/4.0;
V[i] = volume(O,P[F[i].a],P[F[i].b],P[F[i].c]);
}
Point3D res = Fc[0],tmp;
double m = V[0];
for (int i = 1; i < cnt; i++)
{
if (fabs(m+V[i]) < eps)
V[i] += eps;
tmp.x = (m*res.x+V[i]*Fc[i].x)/(m+V[i]);
tmp.y = (m*res.y+V[i]*Fc[i].y)/(m+V[i]);
tmp.z = (m*res.z+V[i]*Fc[i].z)/(m+V[i]);
m += V[i];
res = tmp;
}
return res;
}
};
bool LineCrossPlane(Plane3D pla,Line3D l,Point3D &p) {
Point3D norm=Point3D(pla.A,pla.B,pla.C);
Point3D d=l.d.normalize();
typec rt=(norm^d);
if (sgn(rt)==0) return false;
double t=(-pla.D-(norm^l.s))/rt;
p=l.s+(d*t);
return true;
}
struct Square{ //正方形给对角线顶点求另两顶点
Point2D A,B,C,D;
Square(){}
Square(Point2D a, Point2D b,Point2D c,Point2D d):A(a),B(b),C(c),D(d){}
Square(Point2D a,Point2D c){
A=a,C=c;
DoneSq(A,C,B,D);
}
void DoneSq(Point2D a, Point2D c,Point2D & b,Point2D &d){
typec x,y,mx, my;
mx = (a.x+c.x)/2.0, my = (a.y+c.y)/2.0;
x = a.x - mx; y = a.y - my;
b.x = -y + mx; b.y = x + my;
x = c.x - mx; y = c.y - my;
d.x = - y + mx; d.y = x + my;
}
};
struct Tetra{ //三棱锥体积计算
typec a,b,c,d,e,f;
Tetra(){}
Tetra(typec _a,typec _b,typec _c,typec _d,typec _e,typec _f):a(_a),b(_b),c(_c),d(_d),e(_e),f(_f){}
void input(){
scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&e,&f);
}
typec getVolume(){
double tr1,tr2,tr3,tr4,temp,ans;
tr1=acos( (c*c+b*b-f*f) / (2.0*b*c) );
tr2=acos( (a*a+c*c-e*e) / (2.0*a*c) );
tr3=acos( (a*a+b*b-d*d) / (2.0*a*b) );
tr4=(tr1+tr2+tr3)/2.0;
temp=sqrt( sin(tr4) * sin(tr4-tr1) * sin(tr4-tr2) * sin(tr4-tr3));
ans=a*b*c*temp/3.0;
return ans;
}
};
struct Triangle{
Point2D a,b,c;
Triangle(){}
Triangle(Point2D _a,Point2D _b,Point2D _c):a(_a),b(_b),c(_c){}
//外心
Point2D circumcenter(){
Line2D u,v;
u.a.x=(a.x+b.x)/2;
u.a.y=(a.y+b.y)/2;
u.b.x=u.a.x-a.y+b.y;
u.b.y=u.a.y+a.x-b.x;
v.a.x=(a.x+c.x)/2;
v.a.y=(a.y+c.y)/2;
v.b.x=v.a.x-a.y+c.y;
v.b.y=v.a.y+a.x-c.x;
return u.crosspoint(v);
}
//内心
Point2D incenter(){
Line2D u,v;
double m,n;
u.a=a;
m=atan2(b.y-a.y,b.x-a.x);
n=atan2(c.y-a.y,c.x-a.x);
u.b.x=u.a.x+cos((m+n)/2);
u.b.y=u.a.y+sin((m+n)/2);
v.a=b;
m=atan2(a.y-b.y,a.x-b.x);
n=atan2(c.y-b.y,c.x-b.x);
v.b.x=v.a.x+cos((m+n)/2);
v.b.y=v.a.y+sin((m+n)/2);
return u.crosspoint(v);
}
//垂心
Point2D perpencenter(){
Line2D u,v;
u.a=c;
u.b.x=u.a.x-a.y+b.y;
u.b.y=u.a.y+a.x-b.x;
v.a=b;
v.b.x=v.a.x-a.y+c.y;
v.b.y=v.a.y+a.x-c.x;
return u.crosspoint(v);
}
//重心
//到三角形三顶点距离的平方和最小的点
//三角形内到三边距离之积最大的点
Point2D barycenter(){
Line2D u,v;
u.a.x=(a.x+b.x)/2;
u.a.y=(a.y+b.y)/2;
u.b=c;
v.a.x=(a.x+c.x)/2;
v.a.y=(a.y+c.y)/2;
v.b=b;
return u.crosspoint(v);
}
//费马点
//到三角形三顶点距离之和最小的点
Point2D fermentpoint(){
Point2D u,v;
double step=fabs(a.x)+fabs(a.y)+fabs(b.x)+fabs(b.y)+fabs(c.x)+fabs(c.y);
int i,j,k;
u.x=(a.x+b.x+c.x)/3;
u.y=(a.y+b.y+c.y)/3;
while (step>1e-10)
for (k=0;k<10;step/=2,k++)
for (i=-1;i<=1;i++)
for (j=-1;j<=1;j++){
v.x=u.x+step*i;
v.y=u.y+step*j;
if (u.distance(a)+u.distance(b)+u.distance(c)>v.distance(a)+v.distance(b)+v.distance(c))
u=v;
}
return u;
}
};
struct Ball{
//计算圆心角lat表示纬度,-90<=w<=90,lng表示经度
//返回两点所在大圆劣弧对应圆心角,0<=angle<=pi
double angle(double lng1,double lat1,double lng2,double lat2){
double dlng=fabs(lng1-lng2)*pi/180;
while (dlng>=pi+pi)
dlng-=pi+pi;
if (dlng>pi)
dlng=pi+pi-dlng;
lat1*=pi/180,lat2*=pi/180;
return acos(cos(lat1)*cos(lat2)*cos(dlng)+sin(lat1)*sin(lat2));
}
//计算距离,r为球半径
double line_dist(double r,double lng1,double lat1,double lng2,double lat2){
double dlng=fabs(lng1-lng2)*pi/180;
while (dlng>=pi+pi)
dlng-=pi+pi;
if (dlng>pi)
dlng=pi+pi-dlng;
lat1*=pi/180,lat2*=pi/180;
return r*sqrt(2-2*(cos(lat1)*cos(lat2)*cos(dlng)+sin(lat1)*sin(lat2)));
}
//计算球面距离,r为球半径
inline double sphere_dist(double r,double lng1,double lat1,double lng2,double lat2){
return r*angle(lng1,lat1,lng2,lat2);
}
};
} ;
Raw
xiaodao.cpp
/** Micro Mezz Macro Flation -- Overheated Economy ., Last Update: Aug. 4th 2013 **/ //{
/** Header .. **/ //{
#pragma comment(linker, "/STACK:36777216")
//#pragma GCC optimize ("O2")
#define LOCAL
//#include "testlib.h"
#include <functional>
#include <algorithm>
#include <iostream>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <numeric>
#include <cstring>
#include <climits>
#include <cassert>
#include <complex>
#include <cstdio>
#include <string>
#include <vector>
#include <bitset>
#include <queue>
#include <stack>
#include <cmath>
#include <ctime>
#include <list>
#include <set>
#include <map>
//#include <tr1/unordered_set>
//#include <tr1/unordered_map>
//#include <array>
using namespace std;
#define REP(i, n) for (int i=0;i<int(n);++i)
#define FOR(i, a, b) for (int i=int(a);i<int(b);++i)
#define DWN(i, b, a) for (int i=int(b-1);i>=int(a);--i)
#define REP_1(i, n) for (int i=1;i<=int(n);++i)
#define FOR_1(i, a, b) for (int i=int(a);i<=int(b);++i)
#define DWN_1(i, b, a) for (int i=int(b);i>=int(a);--i)
#define REP_C(i, n) for (int n____=int(n),i=0;i<n____;++i)
#define FOR_C(i, a, b) for (int b____=int(b),i=a;i<b____;++i)
#define DWN_C(i, b, a) for (int a____=int(a),i=b-1;i>=a____;--i)
#define REP_N(i, n) for (i=0;i<int(n);++i)
#define FOR_N(i, a, b) for (i=int(a);i<int(b);++i)
#define DWN_N(i, b, a) for (i=int(b-1);i>=int(a);--i)
#define REP_1_C(i, n) for (int n____=int(n),i=1;i<=n____;++i)
#define FOR_1_C(i, a, b) for (int b____=int(b),i=a;i<=b____;++i)
#define DWN_1_C(i, b, a) for (int a____=int(a),i=b;i>=a____;--i)
#define REP_1_N(i, n) for (i=1;i<=int(n);++i)
#define FOR_1_N(i, a, b) for (i=int(a);i<=int(b);++i)
#define DWN_1_N(i, b, a) for (i=int(b);i>=int(a);--i)
#define REP_C_N(i, n) for (int n____=(i=0,int(n));i<n____;++i)
#define FOR_C_N(i, a, b) for (int b____=(i=0,int(b);i<b____;++i)
#define DWN_C_N(i, b, a) for (int a____=(i=b-1,int(a));i>=a____;--i)
#define REP_1_C_N(i, n) for (int n____=(i=1,int(n));i<=n____;++i)
#define FOR_1_C_N(i, a, b) for (int b____=(i=1,int(b);i<=b____;++i)
#define DWN_1_C_N(i, b, a) for (int a____=(i=b,int(a));i>=a____;--i)
#define ECH(it, A) for (__typeof(A.begin()) it=A.begin(); it != A.end(); ++it)
#define REP_S(i, str) for (char*i=str;*i;++i)
#define REP_L(i, hd, nxt) for (int i=hd;i;i=nxt[i])
#define REP_G(i, u) REP_L(i,hd[u],suc)
#define REP_SS(x, s) for (int x=s;x;x=(x-1)&s)
#define DO(n) for ( int ____n = n; ____n-->0; )
#define REP_2(i, j, n, m) REP(i, n) REP(j, m)
#define REP_2_1(i, j, n, m) REP_1(i, n) REP_1(j, m)
#define REP_3(i, j, k, n, m, l) REP(i, n) REP(j, m) REP(k, l)
#define REP_3_1(i, j, k, n, m, l) REP_1(i, n) REP_1(j, m) REP_1(k, l)
#define REP_4(i, j, k, ii, n, m, l, nn) REP(i, n) REP(j, m) REP(k, l) REP(ii, nn)
#define REP_4_1(i, j, k, ii, n, m, l, nn) REP_1(i, n) REP_1(j, m) REP_1(k, l) REP_1(ii, nn)
#define ALL(A) A.begin(), A.end()
#define LLA(A) A.rbegin(), A.rend()
#define CPY(A, B) memcpy(A, B, sizeof(A))
#define INS(A, P, B) A.insert(A.begin() + P, B)
#define ERS(A, P) A.erase(A.begin() + P)
#define BSC(A, x) (lower_bound(ALL(A), x) - A.begin())
#define CTN(T, x) (T.find(x) != T.end())
#define SZ(A) int((A).size())
#define PB push_back
#define MP(A, B) make_pair(A, B)
#define PTT pair<T, T>
#define fi first
#define se second
#define re real()
#define im imag()
#define Rush for(int ____T=RD(); ____T--;)
#define Display(A, n, m) { \
REP(i, n){ \
REP(j, m-1) cout << A[i][j] << " "; \
cout << A[i][m-1] << endl; \
} \
}
#define Display_1(A, n, m) { \
REP_1(i, n){ \
REP_1(j, m-1) cout << A[i][j] << " "; \
cout << A[i][m] << endl; \
} \
}
string __file__(){
string res = __FILE__;
int r = SZ(res) - 1; while (res[r] != '.') --r;
int l = r - 1; while (res[l] != '\\') --l; ++l;
return res.substr(l, r-l);
}
void Exec(string a, string b, string c){
if (b.empty()) b = __file__();
string cmd = a + ' ' + b + '.' + c;
system(cmd.c_str());
}
void Ruby(string file = ""){Exec("ruby", file, "rb");}
void Python(string file = ""){Exec("python", file, "py");}
void Haskell(string file = ""){Exec("runghc", file, "hs");}
void Pascal(string file = ""){Exec("pascal", file, "pas");}
void Ocaml(string file = ""){Exec("ocaml", file, "ml");}
typedef long long LL;
//typedef long double DB;
typedef double DB;
typedef unsigned UINT;
typedef unsigned long long ULL;
typedef vector<int> VI;
typedef vector<char> VC;
typedef vector<string> VS;
typedef vector<LL> VL;
typedef vector<DB> VF;
typedef set<int> SI;
typedef set<string> SS;
typedef map<int, int> MII;
typedef map<string, int> MSI;
typedef pair<int, int> PII;
typedef pair<LL, LL> PLL;
typedef vector<PII> VII;
typedef vector<VI> VVI;
typedef vector<VII> VVII;
template<class T> inline T& RD(T &);
template<class T> inline void OT(const T &);
//inline int RD(){int x; return RD(x);}
inline LL RD(){LL x; return RD(x);}
inline DB& RF(DB &);
inline DB RF(){DB x; return RF(x);}
inline char* RS(char *s);
inline char& RC(char &c);
inline char RC();
inline char& RC(char &c){scanf(" %c", &c); return c;}
inline char RC(){char c; return RC(c);}
//inline char& RC(char &c){c = getchar(); return c;}
//inline char RC(){return getchar();}
template<class T> inline T& RDD(T &x){
char c; for (c = getchar(); c < '-'; c = getchar());
if (c == '-'){x = '0' - getchar(); for (c = getchar(); '0' <= c && c <= '9'; c = getchar()) x = x * 10 + '0' - c;}
else {x = c - '0'; for (c = getchar(); '0' <= c && c <= '9'; c = getchar()) x = x * 10 + c - '0';}
return x;
}
inline LL RDD(){LL x; return RDD(x);}
template<class T0, class T1> inline T0& RD(T0 &x0, T1 &x1){RD(x0), RD(x1); return x0;}
template<class T0, class T1, class T2> inline T0& RD(T0 &x0, T1 &x1, T2 &x2){RD(x0), RD(x1), RD(x2); return x0;}
template<class T0, class T1, class T2, class T3> inline T0& RD(T0 &x0, T1 &x1, T2 &x2, T3 &x3){RD(x0), RD(x1), RD(x2), RD(x3); return x0;}
template<class T0, class T1, class T2, class T3, class T4> inline T0& RD(T0 &x0, T1 &x1, T2 &x2, T3 &x3, T4 &x4){RD(x0), RD(x1), RD(x2), RD(x3), RD(x4); return x0;}
template<class T0, class T1, class T2, class T3, class T4, class T5> inline T0& RD(T0 &x0, T1 &x1, T2 &x2, T3 &x3, T4 &x4, T5 &x5){RD(x0), RD(x1), RD(x2), RD(x3), RD(x4), RD(x5); return x0;}
template<class T0, class T1, class T2, class T3, class T4, class T5, class T6> inline T0& RD(T0 &x0, T1 &x1, T2 &x2, T3 &x3, T4 &x4, T5 &x5, T6 &x6){RD(x0), RD(x1), RD(x2), RD(x3), RD(x4), RD(x5), RD(x6); return x0;}
template<class T0, class T1> inline void OT(const T0 &x0, const T1 &x1){OT(x0), OT(x1);}
template<class T0, class T1, class T2> inline void OT(const T0 &x0, const T1 &x1, const T2 &x2){OT(x0), OT(x1), OT(x2);}
template<class T0, class T1, class T2, class T3> inline void OT(const T0 &x0, const T1 &x1, const T2 &x2, const T3 &x3){OT(x0), OT(x1), OT(x2), OT(x3);}
template<class T0, class T1, class T2, class T3, class T4> inline void OT(const T0 &x0, const T1 &x1, const T2 &x2, const T3 &x3, const T4 &x4){OT(x0), OT(x1), OT(x2), OT(x3), OT(x4);}
template<class T0, class T1, class T2, class T3, class T4, class T5> inline void OT(const T0 &x0, const T1 &x1, const T2 &x2, const T3 &x3, const T4 &x4, const T5 &x5){OT(x0), OT(x1), OT(x2), OT(x3), OT(x4), OT(x5);}
template<class T0, class T1, class T2, class T3, class T4, class T5, class T6> inline void OT(const T0 &x0, const T1 &x1, const T2 &x2, const T3 &x3, const T4 &x4, const T5 &x5, const T6 &x6){OT(x0), OT(x1), OT(x2), OT(x3), OT(x4), OT(x5), OT(x6);}
inline char& RC(char &a, char &b){RC(a), RC(b); return a;}
inline char& RC(char &a, char &b, char &c){RC(a), RC(b), RC(c); return a;}
inline char& RC(char &a, char &b, char &c, char &d){RC(a), RC(b), RC(c), RC(d); return a;}
inline char& RC(char &a, char &b, char &c, char &d, char &e){RC(a), RC(b), RC(c), RC(d), RC(e); return a;}
inline char& RC(char &a, char &b, char &c, char &d, char &e, char &f){RC(a), RC(b), RC(c), RC(d), RC(e), RC(f); return a;}
inline char& RC(char &a, char &b, char &c, char &d, char &e, char &f, char &g){RC(a), RC(b), RC(c), RC(d), RC(e), RC(f), RC(g); return a;}
inline DB& RF(DB &a, DB &b){RF(a), RF(b); return a;}
inline DB& RF(DB &a, DB &b, DB &c){RF(a), RF(b), RF(c); return a;}
inline DB& RF(DB &a, DB &b, DB &c, DB &d){RF(a), RF(b), RF(c), RF(d); return a;}
inline DB& RF(DB &a, DB &b, DB &c, DB &d, DB &e){RF(a), RF(b), RF(c), RF(d), RF(e); return a;}
inline DB& RF(DB &a, DB &b, DB &c, DB &d, DB &e, DB &f){RF(a), RF(b), RF(c), RF(d), RF(e), RF(f); return a;}
inline DB& RF(DB &a, DB &b, DB &c, DB &d, DB &e, DB &f, DB &g){RF(a), RF(b), RF(c), RF(d), RF(e), RF(f), RF(g); return a;}
inline void RS(char *s1, char *s2){RS(s1), RS(s2);}
inline void RS(char *s1, char *s2, char *s3){RS(s1), RS(s2), RS(s3);}
template<class T0,class T1>inline void RDD(T0&a, T1&b){RDD(a),RDD(b);}
template<class T0,class T1,class T2>inline void RDD(T0&a, T1&b, T2&c){RDD(a),RDD(b),RDD(c);}
template<class T> inline void RST(T &A){memset(A, 0, sizeof(A));}
template<class T> inline void FLC(T &A, int x){memset(A, x, sizeof(A));}
template<class T> inline void CLR(T &A){A.clear();}
template<class T0, class T1> inline void RST(T0 &A0, T1 &A1){RST(A0), RST(A1);}
template<class T0, class T1, class T2> inline void RST(T0 &A0, T1 &A1, T2 &A2){RST(A0), RST(A1), RST(A2);}
template<class T0, class T1, class T2, class T3> inline void RST(T0 &A0, T1 &A1, T2 &A2, T3 &A3){RST(A0), RST(A1), RST(A2), RST(A3);}
template<class T0, class T1, class T2, class T3, class T4> inline void RST(T0 &A0, T1 &A1, T2 &A2, T3 &A3, T4 &A4){RST(A0), RST(A1), RST(A2), RST(A3), RST(A4);}
template<class T0, class T1, class T2, class T3, class T4, class T5> inline void RST(T0 &A0, T1 &A1, T2 &A2, T3 &A3, T4 &A4, T5 &A5){RST(A0), RST(A1), RST(A2), RST(A3), RST(A4), RST(A5);}
template<class T0, class T1, class T2, class T3, class T4, class T5, class T6> inline void RST(T0 &A0, T1 &A1, T2 &A2, T3 &A3, T4 &A4, T5 &A5, T6 &A6){RST(A0), RST(A1), RST(A2), RST(A3), RST(A4), RST(A5), RST(A6);}
template<class T0, class T1> inline void FLC(T0 &A0, T1 &A1, int x){FLC(A0, x), FLC(A1, x);}
template<class T0, class T1, class T2> inline void FLC(T0 &A0, T1 &A1, T2 &A2, int x){FLC(A0, x), FLC(A1, x), FLC(A2, x);}
template<class T0, class T1, class T2, class T3> inline void FLC(T0 &A0, T1 &A1, T2 &A2, T3 &A3, int x){FLC(A0, x), FLC(A1, x), FLC(A2, x), FLC(A3, x);}
template<class T0, class T1, class T2, class T3, class T4> inline void FLC(T0 &A0, T1 &A1, T2 &A2, T3 &A3, T4 &A4, int x){FLC(A0, x), FLC(A1, x), FLC(A2, x), FLC(A3, x), FLC(A4, x);}
template<class T0, class T1, class T2, class T3, class T4, class T5> inline void FLC(T0 &A0, T1 &A1, T2 &A2, T3 &A3, T4 &A4, T5 &A5, int x){FLC(A0, x), FLC(A1, x), FLC(A2, x), FLC(A3, x), FLC(A4, x), FLC(A5, x);}
template<class T0, class T1, class T2, class T3, class T4, class T5, class T6> inline void FLC(T0 &A0, T1 &A1, T2 &A2, T3 &A3, T4 &A4, T5 &A5, T6 &A6, int x){FLC(A0, x), FLC(A1, x), FLC(A2, x), FLC(A3, x), FLC(A4, x), FLC(A5, x), FLC(A6, x);}
template<class T> inline void CLR(priority_queue<T, vector<T>, less<T> > &Q){while (!Q.empty()) Q.pop();}
template<class T> inline void CLR(priority_queue<T, vector<T>, greater<T> > &Q){while (!Q.empty()) Q.pop();}
template<class T> inline void CLR(stack<T> &S){while (!S.empty()) S.pop();}
template<class T0, class T1> inline void CLR(T0 &A0, T1 &A1){CLR(A0), CLR(A1);}
template<class T0, class T1, class T2> inline void CLR(T0 &A0, T1 &A1, T2 &A2){CLR(A0), CLR(A1), CLR(A2);}
template<class T0, class T1, class T2, class T3> inline void CLR(T0 &A0, T1 &A1, T2 &A2, T3 &A3){CLR(A0), CLR(A1), CLR(A2), CLR(A3);}
template<class T0, class T1, class T2, class T3, class T4> inline void CLR(T0 &A0, T1 &A1, T2 &A2, T3 &A3, T4 &A4){CLR(A0), CLR(A1), CLR(A2), CLR(A3), CLR(A4);}
template<class T0, class T1, class T2, class T3, class T4, class T5> inline void CLR(T0 &A0, T1 &A1, T2 &A2, T3 &A3, T4 &A4, T5 &A5){CLR(A0), CLR(A1), CLR(A2), CLR(A3), CLR(A4), CLR(A5);}
template<class T0, class T1, class T2, class T3, class T4, class T5, class T6> inline void CLR(T0 &A0, T1 &A1, T2 &A2, T3 &A3, T4 &A4, T5 &A5, T6 &A6){CLR(A0), CLR(A1), CLR(A2), CLR(A3), CLR(A4), CLR(A5), CLR(A6);}
template<class T> inline void CLR(T &A, int n){REP(i, n) CLR(A[i]);}
template<class T> inline bool EPT(T &a){return a.empty();}
template<class T> inline T& SRT(T &A){sort(ALL(A)); return A;}
template<class T> inline T& RVS(T &A){reverse(ALL(A)); return A;}
template<class T> inline T& UNQ(T &A){A.resize(unique(ALL(SRT(A)))-A.begin());return A;}
template<class T, class C> inline T& SRT(T &A, C B){sort(ALL(A), B); return A;}
//}
/** Constant List .. **/ //{
const int MOD = int(1e9) + 7;
//int MOD = 99990001;
const int INF = 0x3f3f3f3f;
const LL INFF = 0x3f3f3f3f3f3f3f3fLL;
const DB EPS = 1e-9;
const DB OO = 1e20;
const DB PI = acos(-1.0); //M_PI;
const int dx[] = {-1, 0, 1, 0};
const int dy[] = {0, 1, 0, -1};
//}
/** Add On .. **/ //{
// <<= '0. Nichi Joo ., //{
template<class T> inline void checkMin(T &a,const T b){if (b<a) a=b;}
template<class T> inline void checkMax(T &a,const T b){if (a<b) a=b;}
template<class T> inline void checkMin(T &a, T &b, const T x){checkMin(a, x), checkMin(b, x);}
template<class T> inline void checkMax(T &a, T &b, const T x){checkMax(a, x), checkMax(b, x);}
template <class T, class C> inline void checkMin(T& a, const T b, C c){if (c(b,a)) a = b;}
template <class T, class C> inline void checkMax(T& a, const T b, C c){if (c(a,b)) a = b;}
template<class T> inline T min(T a, T b, T c){return min(min(a, b), c);}
template<class T> inline T max(T a, T b, T c){return max(max(a, b), c);}
template<class T> inline T min(T a, T b, T c, T d){return min(min(a, b), min(c, d));}
template<class T> inline T max(T a, T b, T c, T d){return max(max(a, b), max(c, d));}
template<class T> inline T sqr(T a){return a*a;}
template<class T> inline T cub(T a){return a*a*a;}
template<class T> inline T ceil(T x, T y){return (x - 1) / y + 1;}
inline int sgn(DB x){return x < -EPS ? -1 : x > EPS;}
inline int sgn(DB x, DB y){return sgn(x - y);}
inline DB cot(DB x){return 1./tan(x);};
inline DB sec(DB x){return 1./cos(x);};
inline DB csc(DB x){return 1./sin(x);};
//}
// <<= '9. Comutational Geometry .,//{
namespace CG{
struct Po{
DB x,y;Po(DB x=0,DB y=0):x(x),y(y){}
inline void input(){RF(x,y);}
inline friend istream& operator >>(istream& in, Po &p){return in >> p.x >> p.y;}
inline friend ostream& operator <<(ostream& out, Po p){return out << "(" << p.x << ", " << p.y << ")";}
inline bool operator ==(const Po& r)const{return!sgn(x-r.x)&&!sgn(y-r.y);};
inline bool operator !=(const Po& r)const{return sgn(x-r.x)||sgn(y-r.y);}
inline bool operator <(const Po &r) const{return sgn(x,r.x)<0||(!sgn(x,r.x)&&sgn(y,r.y)<0);}
inline bool operator <=(const Po &r) const{return *this<r||*this==r;}
inline bool operator >(const Po &r) const{return !(*this<=r);} //#
inline bool operator >=(const Po &r) const{return !(*this<r);} //#
inline Po operator -()const{return Po(-x,-y);}
inline Po& operator +=(const Po &r){x+=r.x,y+=r.y;return *this;}
inline Po& operator -=(const Po &r){x-=r.x,y-=r.y;return *this;}
inline Po& operator *=(DB k){x*=k,y*=k;return*this;}
inline Po& operator /=(DB k){x/=k,y/=k;return*this;}
inline Po operator +(const Po& r)const{return Po(x+r.x,y+r.y);}
inline Po operator -(const Po& r)const{return Po(x-r.x,y-r.y);}
inline Po operator *(DB k)const{return Po(x*k,y*k);}
inline Po operator /(DB k)const{return Po(x/k,y/k);}
inline DB operator *(const Po&)const; //dot
inline DB operator ^(const Po&)const; //det
inline DB length_sqr()const{return sqr(x)+sqr(y);}
inline DB length()const{return sqrt(length_sqr());}
inline Po unit()const{return *this/length();}
inline bool dgt()const{return !sgn(x)&&!sgn(y);}
inline DB atan()const{return atan2(y,x);}
inline Po sym(){swap(x, y);return*this;}
inline Po symm(){x=-x,y=-y;return*this;}
inline Po lt90(){sym(),x=-x; return*this;}
inline Po rt90(){sym(),y=-y; return*this;}
inline Po rotate(DB alpha, const Po& o = Po()){
x -= o.x, y -= o.y;
return *this = Po(x * cos(alpha) - y * sin(alpha), y * cos(alpha) + x * sin(alpha)) + o;
}
};
inline Po sym(Po p){p.sym();return p;}
inline Po symm(Po p){p.symm();return p;}
inline Po lt90(Po p){p.lt90();return p;}
inline Po rt90(Po p){p.rt90();return p;}
typedef vector<Po> VP;
inline Po operator *(DB k,Po a){return a*k;}
inline DB dot(const DB &x1, const DB &y1, const DB &x2, const DB &y2){return x1 * x2 + y1 * y2;}
inline DB dot(const Po &a, const Po &b){return dot(a.x, a.y, b.x, b.y);}
inline DB dot(const Po &p0, const Po &p1, const Po &p2){return dot(p1 - p0, p2 - p0);}
inline DB det(const DB &x1, const DB &y1, const DB &x2, const DB &y2){return x1 * y2 - x2 * y1;}
inline DB det(const Po &a, const Po &b){return det(a.x, a.y, b.x, b.y);}
inline DB det(const Po &p0, const Po &p1, const Po &p2){return det(p1 - p0, p2 - p0);}
template<class T1, class T2> inline int dett(const T1 &x, const T2 &y){return sgn(det(x, y));}
template<class T1, class T2> inline int dott(const T1 &x, const T2 &y){return sgn(dot(x, y));}
template<class T1, class T2, class T3> inline int dett(const T1 &x, const T2 &y, const T3 &z){return sgn(det(x, y, z));}
template<class T1, class T2, class T3> inline int dott(const T1 &x, const T2 &y, const T3 &z){return sgn(dot(x, y, z));}
template<class T1, class T2, class T3, class T4> inline int dett(const T1 &x, const T2 &y, const T3 &z, const T4 &w){return sgn(det(x, y, z, w));}
template<class T1, class T2, class T3, class T4> inline int dott(const T1 &x, const T2 &y, const T3 &z, const T4 &w){return sgn(dot(x, y, z, w));}
inline DB dist_sqr(const DB &x, const DB &y){return sqr(x) + sqr(y);}
inline DB dist_sqr(const DB &x, const DB &y, const DB &z){return sqr(x) + sqr(y) + sqr(z);}
inline DB dist_sqr(const Po &a, const Po &b){return dist_sqr(a.x - b.x, a.y - b.y);}
template<class T1, class T2> inline DB dist(const T1 &x, const T2 &y){return sqrt(dist_sqr(x, y));}
template<class T1, class T2, class T3> inline DB dist(const T1 &x, const T2 &y, const T3 &z){return sqrt(dist_sqr(x, y, z));}
DB Po::operator *(const Po &r)const{return dot(*this, r);}
DB Po::operator ^(const Po &r)const{return det(*this, r);}
struct Line{
Po a,b;Line(const Po&a=Po(),const Po&b=Po()):a(a),b(b){}
Line(DB x0,DB y0,DB x1,DB y1):a(Po(x0,y0)),b(Po(x1,y1)){}
Line(const Line &l):a(l.a),b(l.b){}
Line(const Po &a,DB alpha):a(a),b(a+Po(cos(alpha),sin(alpha))){}
//Ax+By+C=0
Line(DB A,DB B,DB C){
C=-C;if(!::sgn(A))a=Po(0,C/B),b=Po(1,C/B);
else if(!::sgn(B))a=Po(C/A,0),b=Po(C/A,1);
else a=Po(0,C/B),b=Po(1,(C-A)/B);
}
void input(){a.input(), b.input();}
inline friend istream& operator >>(istream& in, Line& p){return in >> p.a >> p.b;}
inline friend ostream& operator <<(ostream& out, Line p){return out << p.a << "-" << p.b;}
inline Line operator +(const Po& x)const{return Line(a+x,b+x);}
inline Line operator -(const Po& x)const{return Line(a-x,b-x);}
inline Line operator *(DB k)const{return Line(a*k,b*k);}
inline Po operator *(const Line&)const;
inline Po operator *(const Po&)const;
inline DB length_sqr()const{return (b-a).length_sqr();}
inline DB length()const{return (b-a).length();}
inline bool dgt()const{return (b-a).dgt();}
inline int sgn(const Po& p)const{return dett(a, b, p);}
inline int sgn(const Line&)const;
inline bool same_sgn(const Po& p1,const Po& p2)const{return sgn(p1)==sgn(p2);}
inline void get_equation(DB&A,DB&B,DB&C)const{A=a.y-b.y,B=b.x-a.x,C=det(a, b);}
inline void get_intersect(const Line&,Po&)const;
inline Po foot(const Po&)const;
};
inline DB dot(const Line &l1, const Line &l2){return dot(l1.b - l1.a, l2.b - l2.a);}
inline DB det(const Line &l1, const Line &l2){return det(l1.b - l1.a, l2.b - l2.a);}
inline int Line::sgn(const Line&l)const{return dett(*this, l);}
inline Po intersect(const Line &l1, const Line &l2){return l1.a+(l1.b-l1.a)*(det(l2.a,l1.a,l2.b)/det(l2,l1));}
inline Po Line::operator*(const Line& l)const{return intersect(*this, l);}
inline Po Line::foot(const Po& p)const{return *this*Line(p, p+lt90(b-a));}
inline Po Line::operator*(const Po& p)const{return this->foot(p);}
inline void Line::get_intersect(const Line&l,Po&p)const{p=intersect(*this,l);}
struct Seg: public Line{
Seg(const Po&a=Po(),const Po&b=Po()):Line(a,b){}
Seg(DB x0,DB y0,DB x1,DB y1):Line(x0,y0,x1,y1){}
Seg(const Line &l):Line(l){}
Seg(const Po &a,DB alpha):Line(a,alpha){}
Seg(DB A,DB B,DB C):Line(A,B,C){}
inline int sgn(const Po&p)const;
inline bool qrt(const Seg&l)const;
inline int sgn(const Seg&l)const;
};
// quick_rejection_test
inline bool Seg::qrt(const Seg&l)const{
return
min(a.x,b.x) <= max(l.a.x,l.b.x) && min(l.a.x,l.b.x) <= max(a.x,b.x) &&
min(a.y,b.y) <= max(l.a.y,l.b.y) && min(l.a.y,l.b.y) <= max(a.y,b.y);
}
// -1不相交 0相交（不规范） 1相交（规范）
inline int Seg::sgn(const Seg&l)const{
if (!qrt(l)) return -1;
return
(dett(a,b,l.a)*dett(a,b,l.b)<=0 &&
dett(l.a,l.b,a)*dett(l.a,l.b,b)<=0)?1:-1;
int d1=dett(a,b,l.a),d2=dett(a,b,l.b);
int d3=dett(l.a,l.b,a),d4=dett(l.a,l.b,b);
if ((d1^d2)==-2 && (d3^d4)==-2) return 1;
return ((!d1 && dott(l.a-a, l.a-b) <= 0)||(!d2 && dott(l.b-a, l.b-b) <= 0)||
(!d3 && dott(a-l.a, a-l.b) <= 0)||(!d4 && dott(b-l.a, b-l.b) <= 0))?0:-1;
}
inline int Seg::sgn(const Po&p)const{
if (dett(p, a, b)) return -1; if (p == a || p == b) return 0;
return (::sgn(a.x, p.x) * ::sgn(b.x, p.x) < 0 && ::sgn(a.y, p.y) * ::sgn(b.y, p.y) < 0) ? 1 : -1; //#
}
inline DB dist_sqr(const Line &l,const Po &p){Po v0 = l.b - l.a, v1 = p - l.a; return sqr(fabs(det(v0, v1))) / v0.length_sqr();}
inline DB dist_sqr(const Line&l1,const Line&l2){return dett(l1,l2)?0:dist_sqr(l1, l2.a);}
inline DB dist_sqr(const Seg &l,const Po &p){
Po v0 = l.b - l.a, v1 = p - l.a, v2 = p - l.b;
if (sgn(dot(v0, v1)) * sgn(dot(v0, v2)) <= 0) return dist_sqr(Line(l),p);
else return min(v1.length_sqr(), v2.length_sqr());
}
inline DB dist_sqr(const Seg&l1,const Line&l2){
Po v0 = l2.b-l2.a, v1=l1.a-l2.a,v2=l1.b-l2.a;DB c1=det(v0,v1),c2=det(v0,v2);
return sgn(c1)!=sgn(c2) ? 0 : sqr(min(fabs(c1), fabs(c2)))/v0.length_sqr();
}
inline DB dist_sqr(const Seg&l1,const Seg&l2){
if (l1.sgn(l2)) return 0;
else return min(dist_sqr(l2,l1.a), dist_sqr(l2,l1.b), dist_sqr(l1,l2.a), dist_sqr(l1,l2.b));
}
inline Po rotate(Po p,DB alpha,const Po&o = Po()){return p.rotate(alpha, o);}
const int Disjoint = -2, Exscribe = -1, Cross = 0, Inscribe = 1, Contain = 2;
struct Circle{
Po o; DB r; Circle(const Po&o,DB r):o(o),r(r){}
Circle(DB x,DB y,DB r):o(x,y),r(r){}
Circle(const Po&a,const Po&b,const Po&c){ // 外接圆
o = Line((a-b)/2,(a+b)/2+(b-a).lt90())*Line((a-c)/2,(a+c)/2+(c-a).lt90());
r = dist(o, a);
}
void input(){RF(o.x, o.y, r);}
bool operator <(const Circle&c)const{return r<c.r;}
//-1相离 0圆上 1包含
inline int sgn(const Po&p)const{return ::sgn(sqr(r), dist_sqr(o, p));}
//-1相离 0相切 1包含
inline int sgn(const Line&l)const{return ::sgn(sqr(r), dist_sqr(l, o));}
// -2外离 -1外切 0相交 1内切 2包含
inline int sgn(const Circle&c)const{
DB d=dist_sqr(o,c.o);
if (::sgn(sqr(r+c.r),d)<0) return Disjoint;
if (!::sgn(sqr(r+c.r), d)) return Exscribe;
if (!::sgn(sqr(r-c.r), d)) return Inscribe;
if (::sgn(sqr(r-c.r), d)>0) return Contain;
return Cross;
}
void get_intersect(const Line&l,Po&p0,Po&p1){
Po m = l.foot(o), d = (l.b-l.a).unit() * sqrt(sqr(r)-dist_sqr(l, o));
p0 = m + d, p1 = m - d;
}
};
//circum circle 外接圆
//inscribed circle
} using namespace CG;//}
//}
/** I/O Accelerator Interface .. **/ //{
template<class T> inline T& RD(T &x){
//cin >> x;
//scanf("%d", &x);
char c; for (c = getchar(); c < '0'; c = getchar()); x = c - '0'; for (c = getchar(); '0' <= c && c <= '9'; c = getchar()) x = x * 10 + c - '0';
//char c; c = getchar(); x = c - '0'; for (c = getchar(); c >= '0'; c = getchar()) x = x * 10 + c - '0';
return x;
}
inline DB& RF(DB &x){
//cin >> x;
//x = RD();
x = RDD();
//scanf("%lf", &x);
/*char t; while ((t=getchar())==' '||t=='\n'); x = t - '0';
while ((t=getchar())!=' '&&t!='\n'&&t!='.')x*=10,x+=t-'0';
if (t=='.'){DB l=1; while ((t=getchar())!=' '&&t!='\n')l*=0.1,x += (t-'0')*l;}*/
return x;
}
inline char* RS(char *s){
//gets(s);
scanf("%s", s);
return s;
}
LL last_ans; int Case; template<class T> inline void OT(const T &x){
printf("Case %d: ", ++Case);
//printf("%lld\n", x);
printf("%.4f\n", x);
//printf("%d\n", x);
//cout << x << endl;
//last_ans = x;
}
//}
//}/* .................................................................................................................................. */
vector<pair<DB, int> > I; DB bg, ed;
DB v1, v2, v, T, x0;
int n;
inline void add(const Po&p){
checkMin(bg, p.y), checkMax(ed, p.y);
}
int main(){
#ifndef ONLINE_JUDGE
freopen("in.txt", "r", stdin);
//freopen("out.txt", "w", stdout);
#endif
Rush{
RF(v1, v2, v, T, x0), T += (v1*T) / (v2-v1);
CLR(I); Seg l(x0, 0, x0-T*v1, T*v);
Rush{
DB x, y, r, h; Po p0, p1; bg = OO, ed = -OO;
RF(x, y, r, h); Po a(x-r,y),b(x+r,y),c(x,y+h);
Seg ac(a,c),bc(b,c); Circle C(x,y,r);
if (~C.sgn(l)){
C.get_intersect(l, p0, p1);
if (sgn(p0.y, y) <= 0) add(p0);
if (sgn(p1.y, y) <= 0) add(p1);
}
if (~ac.sgn(l)) add(ac*l);
if (~bc.sgn(l)) add(bc*l);
checkMax(bg, 0.0), checkMin(ed, T*v);
if (sgn(bg, ed) < 0) I.PB(MP(bg, 1)), I.PB(MP(ed, -1));
}
SRT(I); DB ans = 0; int cov = 0; REP(i, SZ(I)){
if (cov) ans += I[i].fi - I[i-1].fi;
cov += I[i].se;
}
OT(ans/v);
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment