题解：蓝书习题·第4章·几何问题

二维计算几何模板

Coord getCoord()
{
	lf x,y;
	return scanf("%lf%lf",&x,&y),Coord(x,y);
}

三维计算几何模板

基础题目

Triangle Fun

按照题目意思的直接解法。

int main()
{
	int n;
	for(scanf("%d",&n); n--;)
	{
		Coord A=getCoord(),B=getCoord(),C=getCoord(),
		      D=B+(C-B)/3.0,E=C+(A-C)/3.0,F=A+(B-A)/3.0,
		      P=getLineIntersection(A,D,B,E),
		      Q=getLineIntersection(B,E,C,F),
		      R=getLineIntersection(C,F,A,D);
		printf("%.0f\n",fabs(Area2(P,Q,R)/2));
	}
}

可证明$△ABC$的面积是$△PQR$的7倍，于是有下面的解法。

int main()
{
	int n;
	for(scanf("%d",&n);
	        n--;
	        printf("%.0f\n",
	               fabs(Area2(getCoord(),getCoord(),getCoord())/14)));
}

Determine the Shape

int main()
{
	char s[6][32]=
	{
		"Square",
		"Rectangle",
		"Rhombus",
		"Parallelogram",
		"Trapezium",
		"Ordinary Quadrilateral"
	};
	int t,kase=0,ans;
	for(scanf("%d",&t); t--; printf("Case %d: %s\n",++kase,s[ans]))
	{
		ans=5;
		Coord p[4]= {getCoord(),getCoord(),getCoord(),getCoord()};
		sort(p,p+4,cmpCoord);
		do
		{
			Coord &a=p[0],&b=p[1],&c=p[2],&d=p[3];
			if(!sgn(Cross(b-a,d-c))||!sgn(Cross(a-d,c-b)))ans=min(ans,4);
			if(!sgn(Cross(b-a,d-c))&&!sgn(Cross(a-d,c-b)))
			{
				ans=min(ans,sgn(Dot(c-a,b-d))?3:2);
				if(!sgn(Dot(b-a,c-b)))
					ans=min(ans,sgn(Dot(c-a,b-d))?1:0);
			}
		}
		while(next_permutation(p,p+4,cmpCoord));
	}
}

Athletics Track

int main()
{
	for(lf a,b,c,kase=0;
	        ~scanf("%lf : %lf",&a,&b);
	        printf("Case %.0f: %.5f %.5f\n",kase+=1,a*c,b*c))
		c=200/(sqrt(a*a+b*b)*atan(b/a)+a);
}

Tunnelling the Earth

int main()
{
	int n;
	Sphere earth(Coord3(0,0,0),6371009);
	for(scanf("%d",&n); n--;)
	{
		lf a,b,c,d;
		scanf("%lf%lf%lf%lf",&a,&b,&c,&d);
		Coord3 A(earth.point(a,b)),B(earth.point(c,d));
		printf("%.0f\n",fabs(Angle(A,B)*earth.r)-abs(A-B));
	}
}

二维几何计算

Morley’s Theorem

int main()
{
	int t;
	for(scanf("%d",&t); t--; printf("\n"))
	{
		Coord A=getCoord(),B=getCoord(),C=getCoord(),D;
		for(int i=0; i!=3; ++i)
		{
			lf b=Angle(A-B,C-B),c=Angle(A-C,B-C);
			Line BD(B,Rotate(C-B,b/3)),CD(C,Rotate(B-C,-c/3));
			D=getLineIntersection(BD,CD);
			printf("%.6f %.6f ",D.X,D.Y);
			swap(A,B),swap(B,C);
		}
	}
}

That Nice Euler Circuit

struct CoordCmp
{
	bool operator()(Coord A,Coord B)
	{
		return cmpCoord(A,B);
	}
};
int main()
{
    for(int e,kase=0; scanf("%d",&e)&&e; printf("Case %d: There are %d pieces.\n",++kase,--e))//读入时e本身大1个（最后一笔回到起点），因此上面得到的边数要减去一个
    {
        vector<Coord> p;
        set<Coord,CoordCmp> v;
        for(int i=0; i<e; ++i)
        {
            p.push_back(getCoord());
            v.insert(p.back());
            for(int j=1; j<i; ++j)
                if(SegmentProperIntersection(p[i-1],p[i],p[j-1],p[j]))//不包含端点的相交 
                    v.insert(getLineIntersection(p[i-1],p[i],p[j-1],p[j]));
        }
        for(set<Coord,CoordCmp>::iterator i=v.begin(); i!=v.end(); ++i)
            for(int j=1; j<p.size(); ++j)
                if(onSegment(*i,p[j-1],p[j]))
                    ++e;
        e+=2-v.size();//欧拉定理：平面图的点数V、边数E和面数F满足F=E+2-V
    }
}

Dog Distance

int main()
{
	int I,kase=0,a[2];
	lf maxd,mind,t,len[2];
	Coord p[2];
	vector<Coord> q[2];
	for(scanf("%d",&I); I--; printf("Case %d: %.0f\n",++kase,maxd-mind))
	{
		scanf("%d%d",&a[0],&a[1]);
		for(int i=0; i<2; ++i)
		{
			while(a[i]--)q[i].push_back(getCoord());
			for(int j=len[i]=0; j+1<q[i].size(); ++j)
				len[i]+=abs(q[i][j]-q[i][j+1]);
			p[i]=q[i].back(),q[i].pop_back();
		}
		for(maxd=mind=abs(p[0]-p[1]); !q[0].empty()&&!q[1].empty();)
		{
			Coord v[2]= {q[0].back()-p[0],q[1].back()-p[1]};
			t=min(abs(v[0])/len[0],abs(v[1])/len[1]);
			for(int i=0; i<2; ++i)
				if(q[i].back()==(p[i]+=v[i]*=t*len[i]/abs(v[i])))
					q[i].pop_back();
			maxd=max(maxd,abs(p[0]-p[1]));
			mind=min(mind,DistanceToSegment(p[0],p[1]-v[1]+v[0],p[1]));
		}
	}
}

2D Geometry 110 in 1!

EPS调成1e-6，再小判定圆和直线相切受影响。

lf getlf()
{
	lf x;
	return scanf("%lf",&x),x;
}
Coord getCoord()
{
	lf x=getlf(),y=getlf();
	return Coord(x,y);
}
void printf(lf sol)
{
	printf("%.6lf",sol);
}
void printf(Coord sol)
{
	printf("("),printf(sol.X),printf(","),
	       printf(sol.Y),printf(")");
}
void printf(Circle sol)
{
	printf("("),printf(sol.c.X),printf(","),
	       printf(sol.c.Y),printf(","),
	       printf(sol.r),printf(")");
}
template<typename T>
void printf(vector<T> sol)
{
	sort(sol.begin(),sol.end(),cmpCoord);
	printf("[");
	for(int i=0; i!=sol.size(); ++i)
		printf(i?",":""),printf(sol[i]);
	printf("]");
}
int main()
{
	for(char s[128]; ~scanf("%s",s); printf("\n"))
	{
		if(!strcmp(s,"CircumscribedCircle"))
		{
			Coord p[3]= {getCoord(),getCoord(),getCoord()},
			            p0=(p[1]+p[2])/2.0,
			            p2=(p[1]+p[0])/2.0;
			Line l0(p0,Normal(p[2]-p[1])),
			     l2(p2,Normal(p[0]-p[1]));

			Circle sol(getLineIntersection(l0,l2));
			sol.r=abs(sol.c-p[1]);
			printf(sol);
		}
		else if(!strcmp(s,"InscribedCircle"))
		{
			Coord p[3]= {getCoord(),getCoord(),getCoord()};
			lf arg0=(arg(p[1]-p[0])+arg(p[2]-p[0]))/2,
			       arg2=(arg(p[1]-p[2])+arg(p[0]-p[2]))/2;
			Line l0(p[0],polar(1.0,arg0)),
			     l2(p[2],polar(1.0,arg2));

			Circle sol(getLineIntersection(l0,l2));
			sol.r=fabs(DistanceToLine(sol.c,Line(p[0],p[2]-p[0])));
			printf(sol);
		}
		else if(!strcmp(s,"TangentLineThroughPoint"))
		{
			Circle C(getCoord());
			C.r=getlf();
			Coord P=getCoord();
			vector<Coord> ans;
			getTangents(P,C,ans);

			vector<lf> sol;
			for(int i=0; i!=ans.size(); ++i)
			{
				sol.push_back(arg(ans[i]-P)*180/PI);
				if(sol.back()<0)sol.back()+=180;
			}
			printf(sol);
		}
		else if(!strcmp(s,"CircleThroughAPointAndTangentToALineWithRadius"))
		{
			Circle C(getCoord());
			Coord p[2]= {getCoord(),getCoord()};
			C.r=getlf();
			Coord mov=Normal(p[1]-p[0])*C.r;
			Line L0(p[0]+mov,p[1]-p[0]),L1(p[0]-mov,p[1]-p[0]);

			vector<Coord> sol;
			getLineCircleIntersection(L0,C,sol);
			getLineCircleIntersection(L1,C,sol);
			printf(sol);
		}
		else if(!strcmp(s,"CircleTangentToTwoLinesWithRadius"))
		{
			Coord p[4]= {getCoord(),getCoord(),getCoord(),getCoord()};
			lf r=getlf();
			Coord mov=Normal(p[1]-p[0])*r;
			Line L0(p[0]+mov,p[1]-p[0]),L1(p[0]-mov,p[1]-p[0]);
			mov=Normal(p[3]-p[2])*r;
			Line L2(p[2]+mov,p[3]-p[2]),L3(p[2]-mov,p[3]-p[2]);

			vector<Coord> sol;
			sol.push_back(getLineIntersection(L0,L2));
			sol.push_back(getLineIntersection(L0,L3));
			sol.push_back(getLineIntersection(L1,L2));
			sol.push_back(getLineIntersection(L1,L3));
			printf(sol);
		}
		else if(!strcmp(s,"CircleTangentToTwoDisjointCirclesWithRadius"))
		{
			Circle C[2];
			for(int i=0; i!=2; ++i)
				C[i].c=getCoord(),C[i].r=getlf();
			lf r=getlf();
			C[0].r+=r,C[1].r+=r;

			vector<Coord> sol;
			getCircleIntersection(C[0],C[1],sol);
			printf(sol);
		}
	}
}

几何算法

Board Wrapping

凸包模板题。

int main()
{
	int t,n;
	for(scanf("%d",&t); t--;)
	{
		vector<Coord> p;
		lf x,y,w,h,j,s=0;
		for(scanf("%d",&n); n--; s+=w*h)
		{
			scanf("%lf%lf%lf%lf%lf",&x,&y,&w,&h,&j);
			p.push_back(Coord(x,y)+Rotate(Coord(-w/2,-h/2),-j/180*PI));
			p.push_back(Coord(x,y)+Rotate(Coord(w/2,-h/2),-j/180*PI));
			p.push_back(Coord(x,y)+Rotate(Coord(-w/2,h/2),-j/180*PI));
			p.push_back(Coord(x,y)+Rotate(Coord(w/2,h/2),-j/180*PI));
		}
		printf("%.1f %%\n",s*100/PolygonArea(ConvexHull(p)));
	}
}

Airport

int main()
{
	int t,n,kase=0;
	for(scanf("%d",&t); t--;)
	{
		vector<Coord> p;
		for(scanf("%d",&n); n--;)
		{
			lf x,y;
			scanf("%lf%lf",&x,&y);
			p.push_back(Coord(x,y));
		}
		if(p.size()<3)
		{
			printf("Case #%d: 0.000\n",++kase);
			continue;
		}
		vector<Coord> ch(ConvexHull(p));
		lf ans=INF;
		for(int i=0; i<ch.size(); ++i)
		{
			lf sum=0,len=abs(ch[i]-ch[(i+1)%ch.size()]);
			for(int j=0; j<p.size(); ++j)
				sum+=Cross(ch[i]-p[j],ch[(i+1)%ch.size()]-p[j])/len;
			ans=min(ans,sum);
		}
		printf("Case #%d: %.3f\n",++kase,ans/p.size());
	}
}

The Great Divide

判断两凸包是否相交。

bool ConvexPolygonDisjoint(const vector<Coord> &ch1,const vector<Coord> &ch2)
{
	for(int i=0; i<ch1.size(); ++i)
		if(inPolygon(ch1[i],ch2))
			return 0;
	for(int i=0; i<ch2.size(); ++i)
		if(inPolygon(ch2[i], ch1))
			return 0;
	for(int i=0; i<ch1.size(); ++i)
		for(int j=0; j<ch2.size(); ++j)
			if(SegmentProperIntersection(ch1[i],ch1[(i+1)%ch1.size()],ch2[j],ch2[(j+1)%ch2.size()]))
				return 0;
	return 1;
}
int main()
{
	for(int n[2]; ~scanf("%d%d",&n[0],&n[1])&&n[0]&&n[1];)
	{
		vector<Coord> P[2];
		for(int i=0; i<2; ++i)
			for(int j=0; j<n[i]; ++j)
			{
				lf x,y;
				scanf("%lf%lf",&x,&y);
				P[i].push_back(Coord(x,y));
			}
		printf(ConvexPolygonDisjoint(ConvexHull(P[0]),ConvexHull(P[1]))?"Yes\n":"No\n");
	}
}

Squares

int main()
{
	int t,n;
	for(scanf("%d",&t); t--;)
	{
		scanf("%d",&n);
		vector<Coord> p;
		for(int i=0; i<n; ++i)
		{
			lf x,y,w;
			scanf("%lf%lf%lf",&x,&y,&w);
			p.push_back(Coord(x,y));
			p.push_back(Coord(x,y+w));
			p.push_back(Coord(x+w,y));
			p.push_back(Coord(x+w,y+w));
		}
		vector<Coord> ch(ConvexHull(p));
		lf ans;
		if(ch.size()==1)ans=0;
		else if(ch.size()==2)ans=norm(ch[0]-ch[1]);
		else
			for(int u=ans=0,v=1,k; u<ch.size(); ++u)
				for(;; v=(v+1)%ch.size())
					if(k=sgn(Cross(ch[(u+1)%ch.size()]-ch[u],ch[(v+1)%ch.size()]-ch[v])),k<=0)
					{
						ans=max(ans,norm(ch[u]-ch[v]));
						if(k==0)ans=max(ans,norm(ch[u]-ch[(v+1)%ch.size()]));
						break;
					}
		printf("%.0f\n",ans);
	}
}

Most Distant Point from the Sea

int main()
{
	for(int n; ~scanf("%d", &n)&&n;)
	{
		lf l,r;
		vector<Coord> p;
		for(int i=0; i<n; ++i)
		{
			scanf("%lf%lf",&l,&r);
			p.push_back(Coord(l,r));
		}
		for(l=0,r=INF; r-l>EPS;)
		{
			lf m=(l+r)/2;
			vector<Line> L;
			for(int i=0; i<n; ++i)
			{
				L.push_back(Line(p[i],p[(i+1)%n]-p[i]));
				L.back().p+=m*Normal(L.back().v);
			}
			if(getHalfPlaneIntersection(L).empty())r=m;
			else l=m;
		}
		printf("%.6f\n",l);
	}
}

三维几何

3D Triangles

int main()
{
	int t;
	for(scanf("%d",&t); t--;)
	{
		Coord3 p[6];
		for(int i=0; i<6; ++i)
			scanf("%lf%lf%lf",&p[i].X,&p[i].Y,&p[i].Z);
		printf("%d\n",TriTriIntersection(p,p+3));
	}
}

Ardenia

这一题太特殊了，所有运算都要在分数下进行，因此把完整代码贴出来。
原有的模板几乎无虚修改即可应用，分数类敲sqrt比较麻烦，因此所有距离计算的是是平方后的值（见注释中的原模板和修改之后的代码）。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
struct Fraction
{
	ll num,den;
	Fraction(ll n=0,ll d=1):num(n),den(d)
	{
		d=__gcd(num,den),num/=d,den/=d;
		if(den<0)num=-num,den=-den;
	}
	friend Fraction operator+(const Fraction& A,const Fraction& B)
	{
		ll d=__gcd(A.den,B.den);
		return Fraction(B.den/d*A.num+A.den/d*B.num,A.den/d*B.den);
	}
	Fraction& operator+=(const Fraction &c)
	{
		return *this=*this+c;
	}
	Fraction operator-()const
	{
		Fraction r(*this);
		return r.num=-r.num,r;
	}
	friend Fraction operator-(const Fraction &a,const Fraction &c)
	{
		return -c+a;
	}
	Fraction& operator-=(const Fraction &c)
	{
		return *this=*this-c;
	}
	friend Fraction operator*(const Fraction& A,const Fraction& B)
	{
		return Fraction(A.num*B.num,A.den*B.den);
	}
	Fraction& operator*=(const Fraction &c)
	{
		return *this=*this*c;
	}
	friend Fraction operator/(const Fraction& A,const Fraction& B)
	{
		return Fraction(A.num*B.den,A.den*B.num);
	}
	Fraction& operator/=(const Fraction &c)
	{
		return *this=*this/c;
	}
	friend Fraction operator%(const Fraction &a,const Fraction &c)
	{
		return Fraction(a.num*c.den%(c.num*a.den),a.den*c.den);
	}
	Fraction& operator%=(const Fraction &c)
	{
		return *this=*this%c;
	}
	friend bool operator==(const Fraction &a,const Fraction &b)
	{
		return a.num*b.den==a.den*b.num;
	}
	friend bool operator!=(const Fraction &a,const Fraction &b)
	{
		return !(a==b);
	}
	friend bool operator<(const Fraction &a,const Fraction &b)
	{
		return a.num*b.den<a.den*b.num;
	}
	friend bool operator>(const Fraction &a,const Fraction &b)
	{
		return b<a;
	}
	friend bool operator<=(const Fraction &a,const Fraction &b)
	{
		return !(a>b);
	}
	friend bool operator>=(const Fraction &a,const Fraction &b)
	{
		return !(a<b);
	}
	friend Fraction abs(Fraction f)
	{
		if(f.num<0)f.num=-f.num;
		return f;
	}
	friend ostream& operator<<(ostream &os,const Fraction &f)
	{
		return !f.num?os<<0:
		       f.den==1?os<<f.num:
		       os<<f.num<<'/'<<f.den;
	}
};
typedef Fraction lf;
const lf EPS=1e-9,INF=1e9;
int sgn(lf d)
{
	return (d>EPS)-(d<-EPS);
}
struct Coord3
{
	lf X,Y,Z;
	Coord3(lf X=0,lf Y=0,lf Z=0):X(X),Y(Y),Z(Z) {}
	friend bool operator!=(const Coord3 &a,const Coord3 &b)
	{
		return sgn(a.X-b.X)||sgn(a.Y-b.Y)||sgn(a.Z-b.Z);
	}
	friend bool operator==(const Coord3 &a,const Coord3 &b)
	{
		return !(a!=b);
	}
	Coord3& operator+=(const Coord3 &b)
	{
		return X+=b.X,Y+=b.Y,Z+=b.Z,*this;
	}
	friend Coord3 operator+(Coord3 a,const Coord3 &b)
	{
		return a+=b;
	}
	Coord3& operator-=(const Coord3 &b)
	{
		return X-=b.X,Y-=b.Y,Z-=b.Z,*this;
	}
	friend Coord3 operator-(Coord3 a,const Coord3 &b)
	{
		return a-=b;
	}
	Coord3& operator*=(lf d)
	{
		return X*=d,Y*=d,Z*=d,*this;
	}
	friend Coord3 operator*(Coord3 a,lf d)
	{
		return a*=d;
	}
	friend Coord3 operator*(lf d,Coord3 a)
	{
		return a*=d;
	}
	Coord3& operator/=(lf d)
	{
		return X/=d,Y/=d,Z/=d,*this;
	}
	friend Coord3 operator/(Coord3 a,lf d)
	{
		return a/=d;
	}
};
struct Line3
{
	Coord3 p,v;
	Line3(Coord3 p=Coord3(),Coord3 v=Coord3()):p(p),v(v) {}
	Coord3 point(lf t)const
	{
		return p+v*t;
	}
};
lf Dot(const Coord3& A,const Coord3& B)
{
	return A.X*B.X+A.Y*B.Y+A.Z*B.Z;
}
Coord3 Cross(const Coord3& A,const Coord3& B)
{
	return Coord3(A.Y*B.Z-A.Z*B.Y,A.Z*B.X-A.X*B.Z,A.X*B.Y-A.Y*B.X);
}
lf norm(const Coord3& A)
{
	return Dot(A,A);
}
/*
lf DistanceToLine(Coord3 P,Coord3 A,Coord3 B)//点P到直线AB的距离
{
	Coord3 v1=B-A,v2=P-A;
	return abs(Cross(v1,v2))/abs(v1);
}
*/
lf Distance2ToLine(Coord3 P,Coord3 A,Coord3 B)//点P到直线AB的距离
{
	Coord3 v1=B-A,v2=P-A;
	return norm(Cross(v1,v2))/norm(v1);
}
/*
lf DistanceToSeg(Coord3 P,Coord3 A,Coord3 B)//点到线段的距离
{
	if(A==B)return abs(P-A);
	Coord3 v1=B-A,v2=P-A,v3=P-B;
	if(sgn(Dot(v1,v2))<0)return abs(v2);
	if(sgn(Dot(v1,v3))>0)return abs(v3);
	return fabs(DistanceToLine(P,A,B));
}
*/
lf Distance2ToSeg(Coord3 P,Coord3 A,Coord3 B)//点到线段的距离
{
	if(A==B)return norm(P-A);
	Coord3 v1=B-A,v2=P-A,v3=P-B;
	if(sgn(Dot(v1,v2))<0)return norm(v2);
	if(sgn(Dot(v1,v3))>0)return norm(v3);
	return abs(Distance2ToLine(P,A,B));
}
bool LineDistance3D(Coord3 p1,Coord3 u,Coord3 p2,Coord3 v,lf& s)//求异面直线 p1+s*u与p2+t*v的公垂线对应的s，如果平行|重合，返回0
{
	lf b=Dot(u,u)*Dot(v,v)-Dot(u,v)*Dot(u,v);
	if(!sgn(b))return 0;
	lf a=Dot(u,v)*Dot(v,p1-p2)-Dot(v,v)*Dot(u,p1-p2);
	return s=a/b,1;
}
Coord3 getCoord3()
{
	ll x,y,z;
	return scanf("%lld%lld%lld",&x,&y,&z),Coord3(x,y,z);
}
int main()
{
	int t;
	for(scanf("%d",&t); t--;)
	{
		Coord3 A=getCoord3(),B=getCoord3(),C=getCoord3(),D=getCoord3();
		lf s,t,ans=INF;
		if(LineDistance3D(A,B-A,C,D-C,s)&&0<s&&s<1&&
		        LineDistance3D(C,D-C,A,B-A,t)&&0<t&&t<1)
			ans=norm(Line3(A,B-A).point(s)-Line3(C,D-C).point(t));
		else
		{
			ans=min(ans,Distance2ToSeg(A,C,D));
			ans=min(ans,Distance2ToSeg(B,C,D));
			ans=min(ans,Distance2ToSeg(C,A,B));
			ans=min(ans,Distance2ToSeg(D,A,B));
		}
		printf("%lld %lld\n",ans.num,ans.den);
	}
}

Asteroids

int main()
{
	lf ans=0,x,y,z;
	for(int n,m,flag=1; ~scanf("%d",&n);)
	{
		vector<Coord3> p;
		for(; n--; p.push_back(Coord3(x,y,z)))
			scanf("%lf%lf%lf",&x,&y,&z);
		ConvexPolyhedron c(p);
		ans+=c.dist(c.centroid());
		if(flag^=1)printf("%.5f\n",ans),ans=0;
	}
}

Paperweight