J、Distance to Work

C++(g++ 7.5.0) 解法, 执行用时: 248ms, 内存消耗: 488K, 提交时间: 2022-10-22 13:21:34

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef long double db;
const int N=305;
const db eps=1e-15,pi=acos(-1);
int sgn(db x){return x<-eps?-1:x>eps;}//负数:-1,零:0,正数:1
int cmp(db a,db b){return sgn(a-b);}//小于:-1,等于:0,大于:1
bool solve_eqution(db a,db b,db c,db &x1,db &x2){//解方程 axx+bx+c=0
    db delta=b*b-4*a*c;
    if(sgn(delta)==-1)return 0;
    db t=sqrt(delta);
    db q;
    if(sgn(b)==-1)q=-0.5*(b-t);
    else q=-0.5*(b+t);
    x1=q/a;x2=c/q;if(cmp(x1,x2)==1)swap(x1,x2);
    return 1;
}
struct P{//点、向量类
    db x,y;
    P(db a=0,db b=0):x(a),y(b){}
	P operator + (const P&p)const{return P(x+p.x,y+p.y);}
    P operator - (const P&p)const{return P(x-p.x,y-p.y);}
    P operator * (db a)const{return P(x*a,y*a);}
    P operator / (db a)const{return P(x/a,y/a);}
    db len2(){return x*x+y*y;}
};
db dot(P a,P b){return a.x*b.x+a.y*b.y;}//点积
db cross(P a,P b){return a.x*b.y-a.y*b.x;}//叉积
P lerp(P a,P b,db t){return a*(1-t)+b*t;}//直线上截取某一点连线
db get_angle(P a,P b){return abs(atan2(abs(cross(a,b)),dot(a,b)));}//获得向量OA和OB的夹角或其余角0[0,pi/2]
db getS(P a[],int n){//求多边形面积
    db s=0;a[n+1]=a[1]; 
    for(int i=1;i<=n;++i)s+=cross(a[i],a[i+1]);
    return abs(s)/2;
}
struct circle{//圆类
    P c;db r;
    circle(P p=P(0,0),db x=0):c(p),r(x){}
    P point(db a){return P(c.x+r*cos(a),c.y+r*sin(a));}
};
bool circle_int_line(circle c,P a,P b,db &x1,db &x2){//圆与直线交点,lerp(a,b,x1),lerp(a,b,x2)
    P d=b-a;
    db A=d.len2(),B=dot(d,a-c.c)*2;
	db C=(a-c.c).len2()-c.r*c.r;
    return solve_eqution(A,B,C,x1,x2);
}
db circle_int_triangle(circle c,P a,P b){//圆和三角形(c+a,c+b,o)相交的面积，o是圆心
    if( sgn(cross(a-c.c,b-c.c)) == 0 ) return 0;
    P q[5];int cnt=0;
    db t0,t1; q[cnt++]=a;
    if( circle_int_line(c,a,b,t0,t1) ) {
        if(0<=t0&&t0<=1)q[cnt++]=lerp(a,b,t0);
        if(0<=t1&&t1<=1)q[cnt++]=lerp(a,b,t1);
    }
    q[cnt++]=b;
    db s=0;
    for(int i=1;i<cnt;++i) {
        P z=(q[i-1]+q[i])/2;
        if((z-c.c).len2()<=c.r*c.r)
            s+=abs(cross(q[i-1]-c.c,q[i]-c.c))/2;
        else
            s+=c.r*c.r*get_angle(q[i-1]-c.c,q[i]-c.c)/2;
    }
    return s;
}
db circle_int_poly(circle C,P p[],int n){//圆与多边形相交的面积
    db s = 0;p[n+1]=p[1];
    for(int i=1;i<=n;++i)
        s+=circle_int_triangle(C,p[i],p[i+1])*sgn(cross(p[i]-C.c,p[i+1]-C.c));
    return abs(s);
}
int n;
P a[N];
int main(){
	scanf("%d",&n);
	for(int i=1;i<=n;++i)scanf("%Lf%Lf",&a[i].x,&a[i].y);
	a[n+1]=a[1];
	db S=getS(a,n);
	int q;
	scanf("%d",&q);
	while(q--){
		db x1,y1,p,q;
		scanf("%Lf%Lf%Lf%Lf",&x1,&y1,&p,&q);
		circle C(P(x1,y1),0);
		db l=0,r=5000;
		while(abs(r-l)>=1e-7){
			db m=(l+r)/2;
			C.r=m;
			db now=circle_int_poly(C,a,n);
			if(cmp(q*now,S*(q-p))>=0)r=m;
			else l=m;
		}
		printf("%.10Lf\n",r);
	}
}

C++14(g++5.4) 解法, 执行用时: 1800ms, 内存消耗: 608K, 提交时间: 2018-07-26 16:51:52

#include<bits/stdc++.h>
#include<cmath>
using namespace std;
typedef long double ld;
const ld eps=1e-8;
const ld pi=acos(-1);
const int N=300;
struct point{
    ld x,y;
    point(ld X=0,ld Y=0){x=X,y=Y;}
}a[N],o;
ld all;
int n,m;
ld sqr(ld x){
    return x*x;
}
point operator-(point a,point b){
    return point(a.x-b.x,a.y-b.y);
}
point operator+(point a,point b){
    return point(a.x+b.x,a.y+b.y);
}
point operator*(point a,ld b){
    return point(a.x*b,a.y*b);
}
ld operator*(point a,point b){
    return a.x*b.y-a.y*b.x;
}
point rotate(point a,ld w){
    return point(a.x*cos(w)-a.y*sin(w),a.x*sin(w)+a.y*cos(w));
}
ld cross(point a,point b,point c){
    return (b-a)*(c-a);
}
ld area(point a,point b,point c){
    return fabs((b-a)*(c-a));
}
ld dis(point a,point b){
    return sqrt(sqr(a.x-b.x)+sqr(a.y-b.y));
}
ld dis2(point a,point b){
    return sqr(a.x-b.x)+sqr(a.y-b.y);
}
ld get(point o,point a,point b,ld r){
    if (cross(o,a,b)<0)swap(a,b);
    ld A=dis(o,a),B=dis(o,b),v=area(o,a,b)/2,C=dis(a,b);
    ld w=acos((sqr(A)+sqr(B)-sqr(C))/(2*A*B));
    ld ww=max(acos((sqr(C)+sqr(B)-sqr(A))/(2*C*B)),acos((sqr(A)+sqr(C)-sqr(B))/(2*A*C)));
    ld S=r*r*w/2;
    ld h=v*2/C;
    if (ww>=pi/2){
        if (min(A,B)>=r) return S;
    }else
        if (h>=r)return S;
    point p1,p2;
    if (A<=r)p1=a;
    else{
        ld w1=acos(h/A)-acos(h/r);
        p1=rotate(a-o,w1)*(r/A)+o;
    }
    if (B<=r)p2=b;
    else{
        ld w1=acos(h/B)-acos(h/r);
        p2=rotate(b-o,-w1)*(r/B)+o;
    }
    ld W=acos((dis2(o,p1)+dis2(o,p2)-dis2(p1,p2))/(2*dis(o,p1)*dis(o,p2)));
    return S-r*r*W/2+area(o,p1,p2)/2;
}
ld get(ld r){
    ld S=0;
    for (int i=0;i<n;i++){
        ld v=cross(o,a[i],a[(i+1)%n]);
        if (fabs(v)<eps)continue;
        if (v<0)S-=get(o,a[i],a[(i+1)%n],r); 
        else
            S+=get(o,a[i],a[(i+1)%n],r);
    }
    return fabs(S)/all;
}
int main(){
    //freopen("j.in","r",stdin);
    scanf("%d",&n);
    for (int i=0;i<n;i++){
        int x,y;
        scanf("%d %d",&x,&y);
        a[i]=point(x,y);
    }
    for (int i=1;i<n-1;i++)
        all+=cross(a[0],a[i],a[(i+1)%n])/2;
    all=fabs(all);
    scanf("%d",&m);
    while (m--){
        ld x,y;
        int xx,yy;
        scanf("%d %d",&xx,&yy);
        o=point(xx,yy);
        scanf("%d %d",&xx,&yy);
        x=xx,y=yy;
        x=y-x;
        x/=y;
        ld l=0,r=1e4,mid;
        for(int kk=1;kk<100;++kk){
            if (get(mid=(l+r)/2)<=x)l=mid;
            else
                r=mid;
        }
        printf("%.12f\n",(double)l);
    }
    return 0;
}

C++ 解法, 执行用时: 97ms, 内存消耗: 460K, 提交时间: 2021-10-04 16:45:25

#include<bits/stdc++.h>
#define rep(i,a,b) for(int i=a;i<(int)b;i++)
using namespace std;
struct P{double x,y;};
P operator+(P a,P b){return {a.x+b.x,a.y+b.y};};
P operator-(P a,P b){return {a.x-b.x,a.y-b.y};};
double operator*(P a,P b){return a.x*b.x+a.y*b.y;};
double operator^(P a,P b){return a.x*b.y-a.y*b.x;}

P operator*(P a,double b){return {a.x*b,a.y*b};}
double poly_circle(P c,double r,vector<P>&poly){
	#define len2(a) (a.x*a.x+a.y*a.y)//圆和多边形相交面积
	#define tan(a,b) atan2(a^b,a*b)
	auto tri=[=](P p,P q)->double{//有向三角形和圆相交面积
		P d=q-p;
		double a=d*p/len2(d),b=(len2(p)-r*r)/len2(d),
			r2=r*r/2,det=a*a-b;
		if(det<=0) return tan(p,q)*r2;
		double s=max(0.,-a-sqrt(det)),t=min(1.,-a+sqrt(det));
		if(t<0||s>=1) return tan(p,q)*r2;
		P u=p+d*s,v=p+d*t;
		return tan(p,u)*r2+(u^v)/2+tan(v,q)*r2;
	};
	double sum=0;
	for(int i=0,s=poly.size();i<s;i++)
		sum+=tri(poly[i]-c,poly[(i+1)%s]-c);
	return abs(sum);
}
vector<P> poly;
void solve(P c,double area){
	double lb=0,ub=1e4,mid;
	rep(i,0,60){
		mid=(lb+ub)/2;
		if(poly_circle(c,mid,poly)>area) ub=mid;
		else lb=mid;
	}
	printf("%.9lf\n",ub);
}
int main(){
	int n,m;
	scanf("%d",&n);
	poly.resize(n);
	for(auto&[x,y]:poly) scanf("%lf%lf",&x,&y);
	double s=poly[n-1]^poly[0];
	rep(i,0,n-1) s+=poly[i]^poly[i+1];
	s=abs(s)*.5;
	scanf("%d",&m);
	while(m--){
		P c;double p,q;
		scanf("%lf%lf%lf%lf",&c.x,&c.y,&p,&q);
		solve(c,s*(q-p)/q);
	}
}

详情