【数学】计算几何相关
LiYueLin06 · · 个人记录
提醒：
- 1.叉乘没有交换率，注意前后顺序，点乘叉乘不要搞错！
- 2.下标看清是什么，
i 或j ！p[stk[]] 有没有加stk ！u,v,w ！
T1
CF1025F
几何的一点点结论，极角，（反）三角函数，向量
(直接全文引用算了
如果两个三角形不相交，那么一定存在两条内公切线
可以枚举内公切线
枚举每个点与其他n−1个点组成直线，极角排序
逆时针方向枚举这些直线，同时维护在直线右侧的点的数量
对于每一条直线，利用乘法原理计算方案
还有一些小细节：每一条直线实际上可以组成两个三角形
如上图中，直线CE既可以组成△ABC和△DEF，也可以组成△ABE和△CDF
所以每次计算的时候要乘2
但是每一个三角形又会被两条直线同时计算，所以最后答案又要减半
那么实际上就完全抵消了
同时，我们只枚举极角小于0的直线，防止共线的情况算重。#include<iostream>
#include<stdio.h>
#include<cmath>
#include<algorithm>
using namespace std;
typedef long long LL;
const int N=2005;
const double pi=acos(-1);//π的表示方法
struct vec{
double x,y;//向量就当坐标存了
}a[N];
double b[N];
LL ans;
LL c2(LL x)
{
return x*(x-1)/2;
}
int main()
{
int n;
scanf("%d",&n);
for(int i=1;i<=n;i++)
{
scanf("%lf%lf",&a[i].x,&a[i].y);
}
for(int i=1;i<=n;i++)
{
int cnt=0;
for(int j=1;j<=n;j++)
{
if(i!=j)
{
b[++cnt]=atan2(a[j].y-a[i].y,a[j].x-a[i].x);
//atan2是先纵坐标再横坐标，听说比atan稳定
//变向量用弧度排极角序
}
}
sort(b+1,b+cnt+1);
for(int j=1,k=1;j<=cnt&&b[j]<=0;j++)
{
while(k<=cnt&&b[k]-b[j]<pi)//在直线右侧
{
k++;
}
LL l=k-j-1,r=cnt-l-1;
ans+=c2(l)*c2(r); //左右两边分别选两个点与切线形成三角形
}
}
printf("%lld\n",ans);
return 0;
}T2
P1355 神秘大三角
点在线上，点在形内
能用向量点积叉积就用嘛，简洁明了
#include<bits/stdc++.h>
using namespace std;
#define pii pair<int,int>
#define x first
#define y second
pii a[4],vec[3];
int read(){
int x=0,flag=1;
char c=getchar();
while(c<'0'||c>'9'){
if(c=='-') flag=-1;
c=getchar();
}
while(c>='0'&&c<='9'){
x=x*10+c-'0';
c=getchar();
}
return x*flag;
}
int cha(pii b,pii c){
return b.x*c.y-b.y*c.x;
}
int dian(pii b,pii c){
return b.x*c.x+b.y*c.y;
}
int main(){
int i,j;
for(i=0;i<4;++i){
// a[i].x=read();
// a[i].y=read();
scanf(" (%d,%d)",&a[i].x,&a[i].y);//nb啊但是为什么必须要前面空一格呢
}
for(i=0;i<3;++i){//顶点
if(a[i].x==a[3].x&&a[i].y==a[3].y){
cout<<"4"<<endl;
return 0;
}
}
for(i=0;i<3;++i){//点在线段上
vec[0]={a[i].x-a[3].x,a[i].y-a[3].y};
j=(i==2)?0:i+1;
vec[1]={a[j].x-a[3].x,a[j].y-a[3].y};
vec[2]={a[j].x-a[i].x,a[j].y-a[i].y};
if(cha(vec[2],vec[1])==0&&dian(vec[0],vec[1])<0){
cout<<"3"<<endl;
return 0;
}
}
int op=0,now;
for(i=0;i<3;++i){//点在凸多边形内
vec[0]={a[3].x-a[i].x,a[3].y-a[i].y};
j=(i==2)?0:i+1;
vec[1]={a[j].x-a[3].x,a[j].y-a[3].y};
now=cha(vec[1],vec[0]);
if(op==0){
if(now>0) op=1;
else op=-1;
}
else{
if((now>0&&op==-1)||(now<0&&op==1)){
cout<<"2"<<endl;
return 0;
}
}
}
cout<<"1"<<endl;
return 0;
}T3
P2742 圈奶牛 /【模板】二维凸包
求凸包，先按xy坐标从小到大排序，再用栈分别维护上下凸壳，用叉积判断顺逆时针
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N=2e5+5;//*2,可能是1441,线段？！
#define pii pair<double,double>
#define x first
#define y second
pii p[N],vec[3];
int stk[N],top;
pii operator -(pii a,pii b){
return {a.x-b.x,a.y-b.y};
}
pii operator +(pii a,pii b){
return {a.x+b.x,a.y+b.y};
}
double operator ^(pii a,pii b){//叉积
return (a.x*b.y-a.y*b.x);
}
double operator *(pii a,pii b){//点积
return (a.x*b.x+a.y*b.y);
}
double get_len(pii a){
return sqrt(a*a);
}
int main(){
int n,i,m;
scanf("%d",&n);
for(i=1;i<=n;++i){
scanf("%lf%lf",&p[i].x,&p[i].y);
}
sort(p+1,p+1+n);//默认先x后y从小到大
stk[++top]=1;
for(i=2;i<=n;++i){
while(top>=2&&((p[stk[top]]-p[stk[top-1]])^(p[i]-p[stk[top]]))<0) top--;
stk[++top]=i;
}
m=top;//不要把搞好的下凸壳弹了
for(i=n-1;i>=1;--i){
while(top>=m&&((p[stk[top]]-p[stk[top-1]])^(p[i]-p[stk[top]]))<0) top--;
stk[++top]=i;
}
double ans=0;
for(i=1;i<top;++i){
ans+=get_len(p[stk[i+1]]-p[stk[i]]);
}
printf("%.2lf\n",ans);
return 0;
} T4
CF598C Nearest vectors
1.极角排序，atan2(y,x)被卡了，用叉积，逆时针转，记得要不同象限强行定序，随便从哪里开始都行，不然就成a>b,b>c,c>a了
2.比较角的大小，acos被卡了，用叉积点积，对于向量ab夹的角，把b都转到x轴正半轴上，算得a'的坐标，a'的y值全取正，就可以忽视ab谁前谁后的问题，就可以对两个a'叉积判大小比角度了
#include<bits/stdc++.h>
using namespace std;
const int N=2e5+5;
typedef long long LL;//两次叉积成4次方，还得LL
#define pii pair<LL,LL>
#define Pii pair<int,pair<LL,LL> >
#define x first
#define y second
pii vec[2];
Pii p[N];
LL operator *(pii a,pii b){
return a.x*b.x+a.y*b.y;
}
LL operator ^(pii a,pii b){
return a.x*b.y-a.y*b.x;
}
int whe(pii a){//四个象限四个半轴
if(a.x>0&&a.y>0) return 1;
else if(a.x==0&&a.y>0) return 2;
else if(a.x<0&&a.y>0) return 3;
else if(a.x<0&&a.y==0) return 4;
else if(a.x<0&&a.y<0) return 5;
else if(a.x==0&&a.y<0) return 6;
else if(a.x>0&&a.y<0) return 7;
else if(a.x>0&&a.y==0) return 8;
else return 0;
}
bool cmp(Pii a,Pii b){
if(whe(a.y)!=whe(b.y)){
return whe(a.y)<whe(b.y);
}
else{
return (a.y^b.y)>0ll;
}
}
/*bool blw(pii a){
return (a.y>0||(a.y==0&&a.x>0));
}
bool cmp(Pii a,Pii b){
if(blw(a.y)!=blw(b.y)){//强行按上下分开也行
return blw(a.y);
}
else{
return (a.y^b.y)>0ll;
}
}*/
int main(){
int n,i,j;
scanf("%d",&n);
for(i=1;i<=n;++i){
scanf("%lld%lld",&p[i].y.x,&p[i].y.y);
p[i].x=i;
}
sort(p+1,p+1+n,cmp);
i=n,j=1;
int mni=p[n].x,mnj=p[1].x;
vec[0]={p[i].y*p[j].y,abs(p[i].y^p[j].y)};//y取abs!!!
for(i=1;i<n;++i){
j=i+1;
vec[1]={p[i].y*p[j].y,abs(p[i].y^p[j].y)};//点前叉后!!!
if((vec[1]^vec[0])>0ll){//>0!!!
vec[0]=vec[1];
mni=p[i].x;mnj=p[j].x;
}
}
printf("%d %d\n",mni,mnj);
return 0;
}T5
P3829 信用卡凸包
把信用卡缩一圈变矩形求凸包再加一个圆的周长即可
因为有重复点，踢出凸包栈的时候要写<=0才能踢掉不合法的点！两个相同点向量叉乘另一个会变0！如果这两个相同点不行就踢不掉了！或者写个去重也行
#include<bits/stdc++.h>
using namespace std;
typedef double lb;
const int N=20005;
const lb pi=acos(-1);
#define pii pair<lb,lb>
#define x first
#define y second
pii p[N*4];
int dx[4]={1,1,-1,-1};
int dy[4]={1,-1,1,-1};
int stk[N*4],top;
pii operator +(pii a,pii b){
return {a.x+b.x,a.y+b.y};
}
pii operator -(pii a,pii b){
return {a.x-b.x,a.y-b.y};
}
lb operator *(pii a,pii b){
return a.x*b.x+a.y*b.y;
}
lb operator ^(pii a,pii b){
return a.x*b.y-a.y*b.x;
}
lb get_len(pii a){
return sqrt(a*a);
}
int main(){
int n,i,j,m=0,st;
lb a,b,r,ag,cx,cy,nx,ny,res=0;
cin>>n>>a>>b>>r;
a=a/2.0-r;b=b/2.0-r;
for(i=1;i<=n;++i){
cin>>cx>>cy>>ag;
for(j=0;j<4;++j){
nx=b*dx[j];
ny=a*dy[j];
p[++m]={(nx*cos(ag)-ny*sin(ag))+cx,(ny*cos(ag)+nx*sin(ag))+cy};//向量+点！
}
}
sort(p+1,p+1+m);
stk[++top]=1;
for(i=2;i<=m;++i){
while(top>=2&&((p[stk[top]]-p[stk[top-1]])^(p[i]-p[stk[top]]))<=0) top--;
stk[++top]=i;
}
st=top;
for(i=m-1;i>=1;--i){
while(top>st&&((p[stk[top]]-p[stk[top-1]])^(p[i]-p[stk[top]]))<=0) top--;
stk[++top]=i;
}
res=2.0*pi*r;
for(i=2;i<=top;++i){
res+=get_len(p[stk[i]]-p[stk[i-1]]);
}
printf("%.2lf\n",res);
return 0;
}T6
U130828 【模板】计算几何1
判两线段是否相交，判一条线段两端点是否在另一条线段两侧，用叉积的积小于等于0，等于0是点在线段上，注意精度问题
#include<bits/stdc++.h>
using namespace std;
typedef long double lb;
const int N=10005;
const lb eps=1e-18;
#define pii pair<lb,lb>
#define x first
#define y second
struct lin{
pii p[2];
}l[N];
pii operator +(pii a,pii b){
return {a.x+b.x,a.y+b.y};
}
pii operator -(pii a,pii b){
return {a.x-b.x,a.y-b.y};
}
lb operator *(pii a,pii b){
return a.x*b.x+a.y*b.y;
}
lb operator ^(pii a,pii b){
return a.x*b.y-a.y*b.x;
}
pii vec[3];
bool chk(int i,int j){//判线段i两端点是否在直线j两侧
vec[0]=l[j].p[0]-l[j].p[1];
vec[1]=l[i].p[0]-l[j].p[1];
vec[2]=l[i].p[1]-l[j].p[1];
lb ans=(vec[0]^vec[1])*(vec[0]^vec[2]);
if(abs(ans)<=eps) return 1;
if(ans<0) return 1;
else return 0;
}
int main(){
int n,i,j,q;
cin>>n;
for(i=1;i<=n;++i){
cin>>l[i].p[0].x>>l[i].p[0].y>>l[i].p[1].x>>l[i].p[1].y;
}
cin>>q;
while(q--){
cin>>i>>j;
if(chk(i,j)&&chk(j,i)){
cout<<"Y"<<endl;
}
else{
cout<<"N"<<endl;
}
}
return 0;
}T7
P1452 Beauty Contest G /【模板】旋转卡壳
旋转卡壳求凸包直径
边到点的距离先增后减，点到点的不行！枚举边，双指针找点
如果栈里只有三个点，表示是直线，要特判，不然会跑死！
#include<bits/stdc++.h>
using namespace std;
const int N=200005;
typedef long long LL;
#define pii pair<LL,LL>
#define x first
#define y second
pii p[N],vec[2];
int stk[N],top;
bool us[N];
pii operator -(pii a,pii b){
return {a.x-b.x,a.y-b.y};
}
LL operator *(pii a,pii b){
return a.x*b.x+a.y*b.y;
}
LL operator ^(pii a,pii b){
return a.x*b.y-a.y*b.x;
}
LL len2(pii a){
return a*a;
}
LL dis(int i,int j){
return max(len2(p[stk[i]]-p[stk[j]]),len2(p[stk[i+1]]-p[stk[j]]));
}
double dis_to_l(int i,int j){//j点到i,i+1线距离
vec[0]=p[stk[(i==top)?1:i+1]]-p[stk[i]];
vec[1]=p[stk[j]]-p[stk[i]];
return abs(vec[0]^vec[1])/sqrt(len2(vec[0]));
}
int main(){
int n,i,j,m;
scanf("%d",&n);
for(i=1;i<=n;++i) scanf("%lld%lld",&p[i].x,&p[i].y);
sort(p+1,p+1+n);
stk[++top]=1;
for(i=2;i<=n;++i){
while(top>=2&&((p[stk[top]]-p[stk[top-1]])^(p[i]-p[stk[top]]))<=0){
us[stk[top]]=0;
top--;
}
us[i]=1;
stk[++top]=i;
}
m=top;
for(i=n-1;i>=1;--i){
if(us[i]) continue;
while(top>m&&((p[stk[top]]-p[stk[top-1]])^(p[i]-p[stk[top]]))<=0){
us[stk[top]]=0;
top--;
}
us[i]=1;
stk[++top]=i;
}
LL ans=0;
if(top==3){//一条线！！！不特判会跑死
printf("%lld\n",len2(p[stk[top-1]]-p[stk[top]]));
return 0;
}
for(i=1,j=2;i<=top-1;++i){
while(dis_to_l(i,j)<=dis_to_l(i,(j==top)?1:j+1)){
j=(j==top)?1:j+1;
}
ans=max(ans,dis(i,j));
}
printf("%lld\n",ans);
return 0;
}T8
P3187 最小矩形覆盖
求凸包，枚举边，找左上右三个点
#include<iostream>
#include<stdio.h>
#include<utility>
#include<math.h>
#include<algorithm>
using namespace std;
typedef long double lb;
const int N=100005;
const lb eps=1e-10,inf=1e18,pi=acos(-1);
const lb eps0=0.000005;
#define pii pair<lb,lb>
#define x first
#define y second
pii p[N];
int stk[N],top;
pii operator +(pii a,pii b){return {a.x+b.x,a.y+b.y};}
pii operator -(pii a,pii b){return {a.x-b.x,a.y-b.y};}
lb operator *(pii a,pii b){return a.x*b.x+a.y*b.y;}
lb operator ^(pii a,pii b){return a.x*b.y-a.y*b.x;}
lb len(pii a){return sqrt(a*a);}
lb len2(pii a){return a*a;}
pii sc(lb b,pii a){return {a.x*b,a.y*b};}//数乘
pii si(lb b,pii a){return {a.x/b,a.y/b};}//除法
pii rot(pii a,lb ag){return {a.x*cos(ag)-a.y*sin(ag),a.y*cos(ag)+a.x*sin(ag)};}//旋转
lb shad(pii a,pii b){return (a*b)/len(a);}//投影长
lb chk(lb a){
if(abs(a)<eps0) return 0;
else return a;
}
bool eql(lb a,lb b){
if(abs(b-a)<=eps) return 1;
else return 0;
}
bool btw(int a,int b,int c){//模意义下a在bc之间
if(b<=c) return a>=b&&a<=c;
else return !(a>b&&a<c);
}
int main(){
int n,m,i,j,l,r;
cin>>n;
for(i=1;i<=n;++i) cin>>p[i].x>>p[i].y;
//求凸包
sort(p+1,p+1+n);
stk[++top]=1;
for(i=2;i<=n;++i){//^!!
while(top>=2&&((p[stk[top]]-p[stk[top-1]])^(p[i]-p[stk[top]]))<=0) top--;
stk[++top]=i;
}
m=top;
for(i=n-1;i>=1;--i){
while(top>m&&((p[stk[top]]-p[stk[top-1]])^(p[i]-p[stk[top]]))<=0) top--;
stk[++top]=i;
}
//旋转卡壳
lb h,mn=inf;
pii sqr[4],pl,pr,pq,vl,vr,vi;
for(i=1,j=1,l=1,r=1;i<top;++i){
vi=p[stk[i+1]]-p[stk[i]];//注意^前后顺序！！
while((vi^(p[stk[j%top+1]]-p[stk[i]]))>=(vi^(p[stk[j]]-p[stk[i]]))) j=j%top+1;
//底相同，面积更大，说明高更大
h=(vi^(p[stk[j]]-p[stk[i]]))/len(vi);//高
while((vi*(p[stk[r%top+1]]-p[stk[r]]))>=0) r=r%top+1;//p[stk[]]不要忘stk!!
//右边，投影大于0且增加
pr=p[stk[i]]+sc(((p[stk[r]]-p[stk[i]])*vi)/len2(vi),vi);//右下点
pq=pr+sc(h,rot(si(len(vi),vi),pi/2.0));//右上点
if(btw(l,i,j)) l=j;//先放到左边
while((vi*(p[stk[l%top+1]]-p[stk[l]]))<=0) l=l%top+1;
//左边，投影小于0且减小
pl=p[stk[i]]+sc(((p[stk[l]]-p[stk[i]])*vi)/len2(vi),vi);//左下点
if(len(pr-pl)*h<mn){
mn=len(pr-pl)*h;
sqr[0]=pl;sqr[1]=pr;sqr[2]=pq;
sqr[3]=pl+(pq-pr);
}
}
printf("%.5Lf\n",mn);
mn=sqr[0].y;
int mnid=0;
for(i=1;i<=3;++i){
if(sqr[i].y<mn){
mn=sqr[i].y;mnid=i;
}
else if(eql(sqr[i].y,mn)&&sqr[i].x<sqr[mnid].x){
mnid=i;
}
}
for(i=1,j=mnid;i<=4;++i,j=(j+1)%4){
printf("%.5Lf %.5Lf\n",chk(sqr[j].x),chk(sqr[j].y));//j!!!
}
return 0;
}T9
P4196 凸多边形 / 【模板】半平面交
所有线极角排序，框出所有直线左侧的平面
队列维护所有直线，如果前面的直线交点在该直线右边，把前面的直线弹出，注意队尾队首都弹，先弹队尾，因为队尾对答案的影响更大，最后出来再由队首的线弹队尾的点，因为前面做的时候队尾受的限制比队首少，可能有不行的
最后由交点求多边形面积
#include<iostream>
#include<stdio.h>
#include<utility>
#include<algorithm>
#include<math.h>
using namespace std;
const int N=100005;
typedef double lb;
const lb eps=1e-8;
#define pii pair<lb,lb>
#define x first
#define y second
#define Pii pair<pair<lb,lb>,pair<lb,lb> >
pii p[N],poi[N];
Pii l[N];
int que[N],head,tail;
pii operator +(pii a,pii b){return {a.x+b.x,a.y+b.y};}
pii operator -(pii a,pii b){return {a.x-b.x,a.y-b.y};}
lb operator *(pii a,pii b){return a.x*b.x+a.y*b.y; }
lb operator ^(pii a,pii b){return a.x*b.y-a.y*b.x;}
pii sc(lb b,pii a){//数乘
return {a.x*b,a.y*b};
}
bool chk(pii a,Pii b){//点a在线b左侧
return ((b.y-b.x)^(a-b.x))>0;
}
bool cmp(Pii a,Pii b){
lb t1=atan2((a.x.y-a.y.y),(a.x.x-a.y.x));
lb t2=atan2((b.x.y-b.y.y),(b.x.x-b.y.x));//b!!
if(abs(t2-t1)>eps) return t1<t2;
else{//平行的在左边的放前面
return chk(a.x,b);
}
}
pii get_node(Pii a,Pii b){
// if(a.x==b.x||a.x==b.y) return a.x;
// if(a.y==b.x||a.y==b.y) return a.y;//不然/0了
pii c=a.x;
pii u=a.x-b.x;
pii v=a.y-a.x;
pii w=b.y-b.x;
c=c+sc((w^u)/(v^w),v);//v^w!w!
return c;
}
bool px(Pii a,Pii b){
lb t1=atan2((a.x.y-a.y.y),(a.x.x-a.y.x));
lb t2=atan2((b.x.y-b.y.y),(b.x.x-b.y.x));
return abs(t2-t1)<=eps;
}
int main(){
int n,m,i,j,tot=0;
cin>>n;
while(n--){
cin>>m;
for(i=1+tot;i<=m+tot;++i){
cin>>p[i].x>>p[i].y;
}
for(i=1;i<=m;++i){
j=i%m+1;
l[i+tot]={p[i+tot],p[j+tot]};
}
tot+=m;
}
sort(l+1,l+tot+1,cmp);//极角排序
head=1,tail=0;
for(i=1;i<=tot;++i){//所有点都在直线的左侧
if(px(l[i],l[i-1])) continue;//平行的，后面的在右边，没用
while(tail-head>=1&&chk(poi[tail],l[i])==0) tail--;
while(tail-head>=1&&chk(poi[head+1],l[i])==0) head++;
que[++tail]=i;
if(tail-head>=1){
poi[tail]=get_node(l[que[tail-1]],l[i]);//求直线交点
}
}
while(tail-head>=1&&chk(poi[tail],l[que[head]])==0) tail--;
if(tail-head>=1){
poi[head]=get_node(l[que[head]],l[que[tail]]);
}
lb ans=0;
pii bef=poi[head+1]-poi[head];
for(i=head+2;i<=tail;++i){
pii now=poi[i]-poi[head];
ans=ans+(bef^now);//求多边形面积
bef=now;
}
ans/=2.0;
printf("%.3lf\n",ans);
return 0;
}