计算几何小结

好像又鸽子了几天博客

poj真是神奇的网站……

基础知识

点积

(a·b=a.x*b.x+a.y*b.y)

$a$在$b$上的投影乘以$b$的模长

叉积

(a×b=a.x*b.y-a.y*b.x)

$a,b$围成的平行四边形的有向面积

判断线段相交（跨立实验）

由一条线段的端点向另一条线段两端点做叉积，判断符号是否一致

向量旋转

逆时针旋转$k$度：

inline vector turn (vector a,double k)
{
	return vector(a.x*cos(k)-a.y*sin(k),a.x*sin(k)+a.y*cos(k));
}

判断点在多边形内部

由一点引一条射线，若与多边形有奇数个交点则在多边形内部，否则在多边形外

二维凸包

求凸包

二维平面上给定$n$个点求凸包

显然$y$最小的点应该在凸包上

我们把$y$最小的点先选出来，将其他点极角排序

按照极角序将点入栈，利用斜率单调性判断是否弹栈

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
namespace red{
#define y1 qwq
#define int long long
#define eps (1e-10)
	inline int read()
	{
		int x=0;char ch,f=1;
		for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
		if(ch=='-') f=0,ch=getchar();
		while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
		return f?x:-x;
	}
	const int N=1e5+10;
	int n;
	int st[N],top;
	double ret;
	struct node
	{
		double x,y;
		node (double tx=0,double ty=0){x=tx,y=ty;}
		inline double operator ^ (const node &t) const
		{
			return x*t.y-y*t.x;
		}
		inline node operator - (const node &t) const
		{
			return node(x-t.x,y-t.y);
		}
	}a[N];
	inline double sqr(double x){return x*x;}
	inline double dis(node a,node b)
	{
		return sqrt(sqr(b.x-a.x)+sqr(b.y-a.y));
	}
	inline bool cmp1(node a,node b)
	{
		return a.x==b.x?a.y<b.y:a.x<b.x;
	}
	inline bool cmp2(node q,node b)
	{
		double tmp=(q-a[1])^(b-a[1]);
		if(!tmp) return q.x==b.x?q.y>b.y:q.x<b.x;
		return tmp>0;
	}
	inline void main()
	{
		n=read();
		for(int i=1;i<=n;++i)
		{
			scanf("%lf%lf",&a[i].x,&a[i].y);
		}
		sort(a+1,a+n+1,cmp1);
		int tot=0;
		a[0].x=-1e9+7;
		for(int i=1;i<=n;++i)
		{
			if(a[i].x!=a[tot].x||a[i].y!=a[tot].y) a[++tot]=a[i];
		}
		sort(a+2,a+tot+1,cmp2);
		st[++top]=1;st[++top]=2;
		for(int i=3;i<=tot;++i)
		{
			while(top>2&&((a[st[top-1]]-a[st[top]])^(a[i]-a[st[top]]))>=0) --top;
			st[++top]=i;
		}
		for(int i=1;i<top;++i)
		{
			ret+=dis(a[st[i]],a[st[i+1]]);
		}
		ret+=dis(a[st[top]],a[st[1]]);
		printf("%.2lf",ret);
	}
}
signed main()
{
	red::main();
	return 0;
}

二维凸包面积

利用叉积

inline double cross(node a,node b,node c)
{
	return (b-a)^(c-a);//叉积
}
inline double are()
{
	double ret=0;
	for(int i=1;i<n;++i) ret+=fabs(cross(a[1],a[i],a[i+1]));
	return ret*2;
}

多次判断点在凸包内部

将凸包分成$n-2$个三角形，每次取中间的三角形，用叉积判断当前点在三角形的内部还是左侧或右侧，二分求解

没写过

动态二维凸包

初始有三个点，每次加入一个点，维护凸包

由于凸包只会扩大，所以以初始的三点围成的三角形的重心为基准点求其他点的极角，将点插入合适的位置，再判断两侧的点是否删除

用平衡树/$set$维护

注意每次插入一个点不一定只弹出相邻的两个点……

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<set>
using namespace std;
namespace red{
#define y1 qwq
#define int long long
#define double long double
#define iter set<node>::iterator
#define eps (1e-8)
	inline int read()
	{
		int x=0;char ch,f=1;
		for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
		if(ch=='-') f=0,ch=getchar();
		while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
		return f?x:-x;
	}
	const int N=1e5+10;
	int n,opt;
	struct node
	{
		double x,y,ang;
		node(double tx=0,double ty=0,double tang=0){x=tx,y=ty,ang=tang;}
		friend bool operator < (node a,node b)
		{
			return (a.ang<b.ang)||(a.ang==b.ang&&a.x<b.x);
		}
		inline node operator + (const node b) const
		{
			return node(x+b.x,y+b.y,0);
		}
		inline double operator ^ (const node &b) const
		{
			return x*b.y-y*b.x;
		}
		inline node operator - (const node &b) const
		{
			return node(x-b.x,y-b.y,0);
		}
	}a[N],zx;
	set<node> q;
	iter it1,it2;
	iter pre(iter x)
	{
		return x==q.begin()?--q.end():--x;
	}
	iter nxt(iter x)
	{
		return ++x==q.end()?q.begin():x;
	}
	inline void main()
	{
		n=read();
		for(int i=1;i<=3;++i)
		{
			opt=read();
			a[i].x=read(),a[i].y=read();
			zx=zx+a[i];
			a[i].x*=3,a[i].y*=3;
		}
		for(int i=1;i<=3;++i)
		{
			a[i].ang=atan2(a[i].y-zx.y,a[i].x-zx.x);
			q.insert(a[i]);
		}
		for(int i=4;i<=n;++i)
		{
			opt=read();
			a[i].x=read()*3,a[i].y=read()*3;
			a[i].ang=atan2(a[i].y-zx.y,a[i].x-zx.x);
			it2=q.lower_bound(a[i]);
			if(it2==q.end()) it2=q.begin();
			it1=pre(it2);
			if(opt==1)
			{
				if(((*it2-a[i])^(*it1-a[i]))<=0) continue;
				q.insert(a[i]);
				iter now=q.find(a[i]);
				iter tx=pre(now),ty=pre(tx);
				
				while(((*ty-*now)^(*tx-*now))<=0) q.erase(tx),tx=ty,ty=pre(tx);
				tx=nxt(now),ty=nxt(tx);
				while(((*ty-*now)^(*tx-*now))>=0) q.erase(tx),tx=ty,ty=nxt(tx);
			}
			else
			{
				if(((*it2-a[i])^(*it1-a[i]))<=0) puts("YES");
				else puts("NO");
			}
		}
	}
}
signed main()
{
	red::main();
	return 0;
}

旋转卡壳

组合数学入门题

众所周知，旋转卡壳有$232*2=24$种读法（雾

其实就是单调性优化的枚举？

给定凸包求最大直径

然后我们发现枚举一个点的话另一个点有单调性，所以这个过程是$O(n)$的……

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
namespace red{
#define y1 qwq
#define int long long
#define eps (1e-10)
	inline int read()
	{
		int x=0;char ch,f=1;
		for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
		if(ch=='-') f=0,ch=getchar();
		while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
		return f?x:-x;
	}
	const int N=1e5+10;
	int n;
	int st[N],top;
	double ret;
	struct node
	{
		double x,y;
		node (double tx=0,double ty=0){x=tx,y=ty;}
		inline double operator ^ (const node &t) const
		{
			return x*t.y-y*t.x;
		}
		inline node operator - (const node &t) const
		{
			return node(x-t.x,y-t.y);
		}
	}a[N];
	inline double sqr(double x){return x*x;}
	inline double dis(node a,node b)
	{
		return sqr(b.x-a.x)+sqr(b.y-a.y);
	}
	inline bool cmp1(node a,node b)
	{
		return a.x==b.x?a.y<b.y:a.x<b.x;
	}
	inline bool cmp2(node q,node b)
	{
		double tmp=(q-a[1])^(b-a[1]);
		if(!tmp) return q.x==b.x?q.y>b.y:q.x<b.x;
		return tmp>0;
	}
	inline double are(node a,node b,node c)
	{
		return fabs((a.x*b.y+b.x*c.y+c.x*a.y-a.x*c.y-b.x*a.y-c.x*b.y)/2);
	}
	inline void main()
	{
		n=read();
		for(int i=1;i<=n;++i)
		{
			scanf("%lf%lf",&a[i].x,&a[i].y);
		}
		sort(a+1,a+n+1,cmp1);
		int tot=0;
		a[0].x=-1e9+7;
		for(int i=1;i<=n;++i)
		{
			if(a[i].x!=a[tot].x||a[i].y!=a[tot].y) a[++tot]=a[i];
		}
		sort(a+2,a+tot+1,cmp2);
		st[++top]=1;st[++top]=2;
		for(int i=3;i<=tot;++i)
		{
			while(top>2&&((a[st[top-1]]-a[st[top]])^(a[i]-a[st[top]]))>=0) --top;
			st[++top]=i;
		}
		tot=top;
		for(int i=1;i<=tot;++i) st[++top]=st[i];
		int tmp=3;
		for(int i=1;i<=tot;++i)
		{
			while(are(a[st[i]],a[st[i+1]],a[st[tmp+1]])>are(a[st[i]],a[st[i+1]],a[st[tmp]])) ++tmp;
			ret=max(ret,max(dis(a[st[i]],a[st[tmp]]),dis(a[st[i+1]],a[st[tmp]])));
		}
		printf("%lld",(int)ret);
	}
}
signed main()
{
	red::main();
	return 0;
}

另一个经典应用：给定两凸包求两凸包上点的最近距离

取其中一个凸包的$y$最小的点，另一凸包$y$最大的点开始枚举，方向相反，这样就有单调性了

#include<cstdio>
#include<iostream>
#include<cmath>
#include<algorithm>
using namespace std;
namespace red{
#define int long long
#define ls(p) (p<<1)
#define rs(p) (p<<1|1)
#define eps (1e-9)
	inline int read()
	{
		int x=0;char ch,f=1;
		for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
		if(ch=='-') f=0,ch=getchar();
		while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
		return f?x:-x;
	}
	const int N=5e4+10;
	int n,m;
	struct node
	{
		double x,y;
		node() {};
    	node (double _x,double _y):x(_x),y(_y){}
		inline double operator ^ (const node &t) const {return (x*t.y-y*t.x);}
		inline double operator * (const node &t) const {return (x*t.x+y*t.y);}
		inline node operator - (const node &t) const {return node(x-t.x,y-t.y);}
		inline bool operator < (const node &b) const {return atan2(y,x)<atan2(b.y,b.x);}
	}a[N],b[N];
	inline double sqr(double x) {return x*x;}
	inline double dist(node a,node b)
	{
		return sqrt(sqr(a.x-b.x)+sqr(a.y-b.y));
	}
	double multi(node A,node B,node C)//点积
	{
		return (B-A)*(C-A);
	}
	double cross(node A,node B,node C)//叉积
	{
    	return (B-A)^(C-A);
	}
	inline double disline(node A,node B,node C)
	{
   		if(dist(A,B)<eps) return dist(B,C);
   		if(multi(A,B,C)<-eps) return dist(A,C);//用点积判断是否在线段内部
   		if(multi(B,A,C)<-eps) return dist(B,C);
    	return fabs(cross(A,B,C)/dist(A,B));//点到线段距离
	}
	double mindist(node A,node B,node C,node D)
	{
    	return min(min(disline(A,B,C),disline(A,B,D)),min(disline(C,D,A),disline(C,D,B)));
	}
	inline void check(node a[],int n)
	{
		for(int i=1;i<n-1;++i)
		{
			double tmp=cross(a[i],a[i+1],a[i+2]);
			if(tmp<-eps) return;
			else if(tmp>eps)
			{
				reverse(a+1,a+n+1);
				return;
			}
		}
	}
	inline double solve(node a[],node b[],int n,int m)
	{
		int ymina=1,ymaxb=1;
		for(int i=1;i<=n;++i) if(a[i].y<a[ymina].y) ymina=i;
		for(int i=1;i<=m;++i) if(a[i].y>a[ymaxb].y) ymaxb=i;
		double tmp,ret=1e9+7;
		a[n+1]=a[1],b[m+1]=b[1];
		for(int i=1;i<=n;++i)
		{
			while(tmp=cross(a[ymina+1],b[ymaxb+1],a[ymina])-cross(a[ymina+1],b[ymaxb],a[ymina])<eps)
			{
				++ymaxb;
				if(ymaxb>m) ymaxb=1;
			}
			ret=min(ret,mindist(a[ymina],a[ymina+1],b[ymaxb],b[ymaxb+1]));
			++ymina;
			if(ymina>n) ymina=1;
		}
		return ret;
	}
	inline void main()
	{
		while("haku")
		{
			n=read(),m=read();
			if(!(n|m)) return ;
			for(int i=1;i<=n;++i) scanf("%lf%lf",&a[i].x,&a[i].y);
			for(int i=1;i<=m;++i) scanf("%lf%lf",&b[i].x,&b[i].y);
			check(a,n);check(b,m);
			printf("%.5lf ",min(solve(a,b,n,m),solve(b,a,m,n)));
		}
	}
}
signed main()
{
	red::main();
	return 0;
}

半平面交

有点恶心的东西

给定$n$条直线，每条直线保留左侧部分，求最后留下的面积

首先我们有$n^2$的暴力算法：每加入一条直线，就枚举以前的直线

如果在左侧则保留

如果在右侧把它扔掉

如果相交保留左侧端点，右侧改为交点

当然一般来说$n^2$是不够用的

我们可以先将所有直线按照与$x$轴的夹角排序，夹角相同的保留靠左的

用双端队列维护半平面交

维护方式：先保证队列中至少有两个元素

如果队尾的两个元素的交点在当前直线右侧，把队尾弹出

队首同理

但是记住一定要先弹出队尾，因为新加入的元素与一开始加入的元素围成的面积更小

如在这张图中，如果标号就是加入队列的顺序，那么当$3$加入的时候，先判断$1$就会踢$1$，先判断$2$就会踢$2$，但是新加入的元素和队首围成的面积会更小，所以应该先踢队尾

#include<bits/stdc++.h>
using namespace std;
namespace red{
#define int long long
#define ls(p) (p<<1)
#define rs(p) (p<<1|1)
#define eps (1e-8)
	inline int read()
	{
		int x=0;char ch,f=1;
		for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
		if(ch=='-') f=0,ch=getchar();
		while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
		return f?x:-x;
	}
	const int N=5000;
	int n,m,tot,sum,head,tail;
	struct node
	{
		double x,y;
		inline double operator ^ (const node &t) const{return x*t.y-y*t.x;}
		inline node operator - (const node &t) const{return (node){x-t.x,y-t.y};}
		inline node operator + (const node &t) const{return (node){x+t.x,y+t.y};}
		inline node operator * (const double &t) const{return (node){x*t,y*t};}
	}a[N],c[N];
	struct seg
	{
		node a,b;
		double k;
		seg(){}
		seg(const node &aa,const node &bb):a(aa),b(bb){k=atan2(b.y,b.x);}
		inline bool operator < (const seg &t) const
		{
			return k<t.k;
		}
	}q[N],que[N];
	inline double cross(node a,node b,node c)
	{
		return (b-a)^(c-a);
	}
	node get_node(seg A,seg B)//求交点
	{
    	node C=A.a-B.a;
    	double t=(B.b^C)/(A.b^B.b);
    	return A.a+A.b*t; 
	}
	inline int dcmp(double x)
	{
		if(fabs(x)<eps) return 0;
		return x>0?1:-1;
	}
	inline double are()
	{
		double are=0;
		for(int i=head;i<tail;++i){ are+=fabs(cross(c[head],c[i],c[i+1])); }
		return are/2;
	}
	inline void work()
	{
		head=tail=1;
		que[tail]=q[1];
		for(int i=2;i<=sum;++i)
		{
			while(head<tail&&(q[i].b^(c[tail-1]-q[i].a))<=eps) --tail;
			while(head<tail&&(q[i].b^(c[head]-q[i].a))<=eps) ++head;
			que[++tail]=q[i];
			if(fabs(que[tail].b^que[tail-1].b)<=eps) //平行保留较左的
			{
				--tail;
				if((que[tail].b^(q[i].a-que[tail].a))>eps) que[tail]=q[i];
			}
			if(head<tail) c[tail-1]=get_node(que[tail-1],que[tail]);//c[i]表示c[i]和c[i+1]的交点
		}
		while(head<tail&&(que[head].b^(c[tail-1]-que[head].a))<=eps) --tail;
		if(tail-head<=1) return;
		c[tail]=get_node(que[head],que[tail]);
	}
	inline void main()
	{
		n=read();
		for(int i=1;i<=n;++i)
		{
			m=read();
			for(int j=1;j<=m;++j)
			{
				a[j].x=read(),a[j].y=read();
			}
			for(int j=1;j<=m;++j)
			{
				++sum;
				int t=j+1;
				if(j==m) t=1;
				q[sum]=seg(a[j],a[t]-a[j]);
			}
		}
		sort(q+1,q+sum+1);
		work();
		printf("%.3lf",are());
	}
}
signed main()
{
	red::main();
	return 0;
}

自适应辛普森积分法

$simple$积分法

求积分$int_^frac{cx+d}{ax+b}dx$

先来说什么是辛普森积分法

假设有函数$g(x)=Ax^2+Bx+C$

(int_{a}^{b}g(x)dx)

微积分基本定理：

(=int_{a}^{b}frac{A}{3}(b^3-a^3)+frac{B}{2}(b^2-a^2)+C(b-a))

大力化简

(=frac{A}{3}(b-a)(a^2+ab+b^2)+frac{B}{2}(b+a)(b-a)+C(b-a))

(=frac{b-a}{6}(2A(a^2+ab+b^2)+3B(b+a)+6C))

(=frac{b-a}{6}(2Aa^2+2Aab+2Ab^2+3Ba+3Bb+6C))

(=frac{b-a}{6}(Aa^2+Ba+C+Ab^2+Bb+C+4A(frac{a+b}{2})^2+4B(frac{a+b}{2})+4C))

(=frac{b-a}{6}(g(a)+g(b)+4g(frac{a+b}{2})))

得到辛普森积分公式

(int_{a}^{b}g(x)dx=frac{b-a}{6}(g(a)+g(b)+4g(frac{a+b}{2})))

不过这是二次函数的积分公式，和我们上面那个式子有啥关系？

我们知道，一个复杂的函数可以近似拟合为一个二次函数

如果$f(x)$的拟合二次函数为$g(x)$，那么

(int_{a}^{b}f(x)dx≈frac{b-a}{6}(f(a)+f(b)+4f(frac{a+b}{2})))

当然拟合肯定存在误差，不过当定义域越小的时候，这个误差也就越小

我们可以把原函数分段，对于每一段求拟合函数$g(x)$的值再加起来，只要分得的段足够小就可以得到正确答案

但是分的太小会导致答案不正确， 太多会降低程序效率

这个时候我们可以让程序自适应

二分递归答案的段数和区间长度，精度达到要求就停止递归返回答案

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<set>
using namespace std;
namespace red{
#define y1 qwq
	inline int read()
	{
		int x=0;char ch,f=1;
		for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
		if(ch=='-') f=0,ch=getchar();
		while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
		return f?x:-x;
	}
	double a,b,c,d,l,r;
	inline double f(double x)
	{
		return (c*x+d)/(a*x+b);
	}
	inline double simpson(double l,double r)
	{
		double mid=(l+r)/2;
		return (f(l)+4*f(mid)+f(r))*(r-l)/6;
	}
	inline double asr(double l,double r,double eps,double ans)
	{
		double mid=(l+r)/2;
		double tl=simpson(l,mid),tr=simpson(mid,r);
		if(fabs(tl+tr-ans)<=15*eps) return tl+tr+(tl+tr-ans)/15;
		return asr(l,mid,eps/2,tl)+asr(mid,r,eps/2,tr);
	}
	inline double asr(double l,double r,double eps)
	{
		return asr(l,r,eps,simpson(l,r));
	}
	inline void main()
	{
		scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&l,&r);
		printf("%.6f
",asr(l,r,1e-6));
	}
}
signed main()
{
	red::main();
	return 0;
}