各种数学知识3
大概有一些计算几何基础
1.向量
因为文化课没学好，还是写写点/叉积吧(丢人
向量点积：对应项乘起来/模长相乘然后乘夹角余弦值
向量叉积：$x1y2-x2y1$，表示向量$1$与其逆时针方向上的向量$2$表示的平行四边形的面积
向量旋转：将$x,y$逆时针旋转$θ$，用复数转：$(x+yi)(sinθ+cosθ)$
2.三点定圆
大力推式子，这里给一段代码
1 void Getcir(a xx,a yy,a zz) 2 { 3 double x1=xx.x,x2=yy.x,x3=zz.x; 4 double y1=xx.y,y2=yy.y,y3=zz.y; 5 double a=x2-x1,b=x3-x1,c=y2-y1,d=y3-y1,bas=2*(b*c-a*d); 6 double sq1=x1*x1+y1*y1,sq2=x2*x2+y2*y2,sq3=x3*x3+y3*y3; 7 nde.x=(c*(sq3-sq1)-d*(sq2-sq1))/bas; 8 nde.y=(b*(sq2-sq1)-a*(sq3-sq1))/bas; 9 rdi=calc(nde,xx); 10 }
3.判断点是否在多边形内
射线法：从该点引一条射线，如果与多边形有偶数个交点则在多边形外，否则在多边形内
实际实现不需要引射线，直接在搞一个一定在多边形外的点弄出一条线段变成线段判交即可
扫描线法：按照询问从上往下扫，每次只取出当前判断的点的下面的部分来判断
优点是当点很多的时候复杂度上界是O(横坐标数目*多边形边数的)，而射线法是O(点数*多边形边数的)；缺点是细节多而不太好写
注意两个方法都要先扫一遍多边形判掉过多边形顶点的情况、
HDU 1756是个模板题
1 #include<cmath> 2 #include<cstdio> 3 #include<cstring> 4 #include<cstdlib> 5 #include<algorithm> 6 #define db double 7 using namespace std; 8 const int N=105; 9 const db eps=1e-8; 10 int n,m; db t1,t2,mx,my,xx[N],yy[N]; 11 bool Equ(db a,db b) 12 { 13 return fabs(a-b)<=eps; 14 } 15 bool Exact(db x1,db y1,db x2,db y2) 16 { 17 db slp=(y1-y2)/(x1-x2),cut=y1-slp*x1; 18 for(int i=1;i<=n;i++) 19 if(Equ(yy[i],slp*xx[i]+cut)) return true; 20 return false; 21 } 22 db Cha(db x0,db y0,db x1,db y1,db x2,db y2) 23 { 24 return (x1-x0)*(y2-y0)-(x2-x0)*(y1-y0); 25 } 26 bool Across(db sx1,db sy1,db ex1,db ey1,db sx2,db sy2,db ex2,db ey2) 27 { 28 return Cha(sx1,sy1,ex1,ey1,sx2,sy2)*Cha(sx1,sy1,ex1,ey1,ex2,ey2)<-eps&& 29 Cha(sx2,sy2,ex2,ey2,sx1,sy1)*Cha(sx2,sy2,ex2,ey2,ex1,ey1)<-eps; 30 } 31 bool Inpoly(db x,db y) 32 { 33 bool inp=false; 34 db mxx=mx+rand(),mxy=my+rand(); 35 while(Exact(x,y,mxx,mxy)) 36 mxx=mx+rand(),mxy=my+rand(); 37 for(int i=1;i<=n;i++) 38 inp^=Across(x,y,mxx,mxy,xx[i],yy[i],xx[i+1],yy[i+1]); 39 return inp; 40 } 41 int main() 42 { 43 srand(20020516); 44 while(scanf("%d",&n)!=EOF) 45 { 46 mx=my=0; 47 for(int i=1;i<=n;i++) 48 { 49 scanf("%lf%lf",&xx[i],&yy[i]); 50 mx=max(mx,xx[i]),my=max(my,yy[i]); 51 } 52 xx[n+1]=xx[1],yy[n+1]=yy[1]; 53 scanf("%d",&m); 54 for(int i=1;i<=m;i++) 55 { 56 scanf("%lf%lf",&t1,&t2); 57 Inpoly(t1,t2)?puts("Yes"):puts("No"); 58 } 59 } 60 return 0; 61 }
4.几何中的欧拉公式
V-E+F=2(三个变量分别为顶点数，边数，面数)
5.凸包(凸猫?)相关
凸的(整个图形在任意一边的延长线一侧)，周长最小(的包住所有点的图形)，斜率单调(字面意思)
求凸包
常见的一种方法：排序增量法
两种排法：极角排序/坐标排序(猫老师推荐后者，精度比较高，不过极角排序可以叉积搞就不丢精度了)
上下凸包分开搞，用一个栈维护，依照斜率单调性，依次加入然后每次按斜率弹栈
1 #include<cmath> 2 #include<cstdio> 3 #include<cstring> 4 #include<algorithm> 5 using namespace std; 6 const int N=10005; 7 struct a{double x,y;}pts[N]; 8 int n,top,stk[N]; double ans; 9 bool cmp(a n1,a n2) 10 { 11 return n1.y==n2.y?n1.x<n2.x:n1.y<n2.y; 12 } 13 double calc(int a,int b) 14 { 15 double difx=pts[a].x-pts[b].x,dify=pts[a].y-pts[b].y; 16 return sqrt(difx*difx+dify*dify); 17 } 18 bool check(int a,int b,int c) 19 { 20 return (pts[a].x-pts[b].x)*(pts[b].y-pts[c].y)<(pts[a].y-pts[b].y)*(pts[b].x-pts[c].x); 21 } 22 int main() 23 { 24 scanf("%d",&n); 25 for(int i=1;i<=n;i++) 26 scanf("%lf%lf",&pts[i].x,&pts[i].y); 27 sort(pts+1,pts+1+n,cmp); 28 stk[1]=1,stk[2]=2,top=2; 29 for(int i=3;i<=n;i++) 30 { 31 while(top>=2&&check(i,stk[top],stk[top-1])) top--; 32 stk[++top]=i; 33 } 34 for(int i=1;i<top;i++) ans+=calc(stk[i],stk[i+1]); 35 stk[1]=1,stk[2]=2,top=2; 36 for(int i=3;i<=n;i++) 37 { 38 while(top>=2&&!check(i,stk[top],stk[top-1])) top--; 39 stk[++top]=i; 40 } 41 for(int i=1;i<top;i++) ans+=calc(stk[i],stk[i+1]); 42 printf("%.2f",ans); 43 return 0; 44 }
6.闵可夫斯基和
即两个点集$n^2$加起来之后得到的集合
7.最小圆覆盖
随机增量法：把给出顺序打乱，然后$n^3$枚举点来做，但是期望是$O(n)$的，不知道为啥
1 #include<cmath> 2 #include<cstdio> 3 #include<cstring> 4 #include<algorithm> 5 using namespace std; 6 const int N=100005; 7 struct a{double x,y;}pts[N],nde; 8 int n,m,idx[N]; double rdi; 9 void Swap(int a,int b) 10 { 11 int tmp=idx[a]; 12 idx[a]=idx[b],idx[b]=tmp; 13 } 14 double calc(a aa,a bb) 15 { 16 double difx=aa.x-bb.x,dify=aa.y-bb.y; 17 return sqrt(difx*difx+dify*dify); 18 } 19 void Getcir(a xx,a yy,a zz) 20 { 21 double x1=xx.x,x2=yy.x,x3=zz.x; 22 double y1=xx.y,y2=yy.y,y3=zz.y; 23 double a=x2-x1,b=x3-x1,c=y2-y1,d=y3-y1,bas=2*(b*c-a*d); 24 double sq1=x1*x1+y1*y1,sq2=x2*x2+y2*y2,sq3=x3*x3+y3*y3; 25 nde.x=(c*(sq3-sq1)-d*(sq2-sq1))/bas; 26 nde.y=(b*(sq2-sq1)-a*(sq3-sq1))/bas; 27 rdi=calc(nde,xx); 28 } 29 int main() 30 { 31 srand(20020513); 32 scanf("%d",&n),m=10*n; 33 for(int i=1;i<=n;i++) idx[i]=i; 34 for(int i=1;i<=m;i++) Swap(rand()%n+1,rand()%n+1); 35 for(int i=1;i<=n;i++) 36 scanf("%lf%lf",&pts[idx[i]].x,&pts[idx[i]].y); 37 nde=pts[1]; 38 for(int i=1;i<=n;i++) 39 if(calc(nde,pts[i])>rdi) 40 { 41 nde=pts[i],rdi=0; 42 for(int j=1;j<=i-1;j++) 43 if(calc(nde,pts[j])>rdi) 44 { 45 nde.x=(pts[i].x+pts[j].x)/2; 46 nde.y=(pts[i].y+pts[j].y)/2; 47 rdi=calc(nde,pts[j]); 48 for(int k=1;k<=j-1;k++) 49 if(calc(nde,pts[k])>rdi) 50 Getcir(pts[i],pts[j],pts[k]); 51 } 52 } 53 printf("%.10f\n%.10f %.10f",rdi,nde.x,nde.y); 54 return 0; 55 }
8.平面最近点对
分治法，没用归并$O(n\log^2n)$，归并的话可以再少一个$log$，具体看这里
1 #include<cmath> 2 #include<cstdio> 3 #include<cstring> 4 #include<algorithm> 5 using namespace std; 6 const int N=200005; 7 const double inf=1e9; 8 struct a{double x,y;}pts[N]; 9 int n,mem[N]; double ans; 10 bool cmp(a n1,a n2) 11 { 12 return n1.y==n2.y?n1.x<n2.x:n1.y<n2.y; 13 } 14 bool com(int n1,int n2) 15 { 16 return pts[n1].x<pts[n2].x; 17 } 18 double calc(int a,int b) 19 { 20 double difx=pts[a].x-pts[b].x,dify=pts[a].y-pts[b].y; 21 return sqrt(difx*difx+dify*dify); 22 } 23 double DC(int l,int r) 24 { 25 if(l==r) return inf; 26 if(l==r-1) return calc(l,r); 27 int mid=(l+r)/2,pt=0; 28 double dis=min(DC(l,mid),DC(mid+1,r)); 29 for(int i=l;i<=r;i++) 30 if(fabs(pts[mid].y-pts[i].y)<=dis) 31 mem[++pt]=i; 32 sort(mem+1,mem+1+pt,com); 33 for(int i=1;i<=pt;i++) 34 for(int j=i+1;j<=pt;j++) 35 { 36 if(pts[mem[j]].x-pts[mem[i]].x>=dis) 37 break; dis=min(dis,calc(mem[i],mem[j])); 38 } 39 return dis; 40 } 41 int main() 42 { 43 scanf("%d",&n); 44 for(int i=1;i<=n;i++) 45 scanf("%lf%lf",&pts[i].x,&pts[i].y); 46 sort(pts+1,pts+1+n,cmp); 47 printf("%.4f",DC(1,n)); 48 return 0; 49 }
9.自适应辛普森法
函数的拟合：用图像相似的初等函数“代替”一些奇怪的函数
辛普森法：下面的这个东西可以完美拟合一次和二次函数，对于其他的就有误差了
$S(l,r)=∫_l^rf(x)dx=\frac{f(l)+4f(\frac{l+r}{2})+f(r)}{6}(f(x)$
那什么是“自适应”呢？
每次判断一下如果误差较大就继续递归，否则返回，这就叫“自适应”
1 #include<cmath> 2 #include<cstdio> 3 #include<cstring> 4 #include<algorithm> 5 using namespace std; 6 const double eps=1e-8; 7 double a,b,c,d,ll,rr; 8 double Func(double x) 9 { 10 return (c*x+d)/(a*x+b); 11 } 12 double Simpson(double l,double r) 13 { 14 double mid=(l+r)/2; 15 return (Func(l)+Func(r)+4*Func(mid))*(r-l)/6; 16 } 17 double Auto_Simpson(double l,double r,double f) 18 { 19 double mid=(l+r)/2,f1=Simpson(l,mid),f2=Simpson(mid,r); 20 double diff=(f1+f2-f)/16; 21 if(fabs(diff)<=eps) 22 return f1+f2+diff; 23 else 24 return Auto_Simpson(l,mid,f1)+Auto_Simpson(mid,r,f2); 25 } 26 int main() 27 { 28 scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&ll,&rr); 29 printf("%.6f",Auto_Simpson(ll,rr,Simpson(ll,rr))); 30 return 0; 31 }
10.伪·动态凸包
离线后倒着加，用set维护凸包
1 #include<set> 2 #include<cmath> 3 #include<cstdio> 4 #include<cstring> 5 #include<algorithm> 6 using namespace std; 7 const int N=200005; 8 int T,n,m,cx,cy,t1,t2,p1,p2; 9 int del[N],qry[N]; double ans,outp[N]; 10 struct a 11 { 12 int xp,yp; 13 }era[N],pts[N]; 14 bool operator < (a x,a y) 15 { 16 return x.xp==y.xp?x.yp<y.yp:x.xp<y.xp; 17 } 18 set<a> st; 19 double Dist(a x,a y) 20 { 21 int x1=x.xp,x2=y.xp,y1=x.yp,y2=y.yp; 22 return sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)); 23 } 24 int Cha(a x,a y) 25 { 26 return x.xp*y.yp-x.yp*y.xp; 27 } 28 a Diff(a x,a y) 29 { 30 x.xp-=y.xp,x.yp-=y.yp; 31 return x; 32 } 33 void Add(int xx,int yy) 34 { 35 a p=(a){xx,yy}; 36 set<a> ::iterator lt,rt,tt; 37 lt=rt=st.lower_bound(p),lt--; 38 a lp=*lt,rp=*rt; 39 if(Cha(Diff(rp,lp),Diff(p,lp))<0) return; 40 ans-=Dist(lp,rp),st.insert(p); 41 while(true) 42 { 43 tt=rt++; 44 if(rt==st.end()||Cha(Diff(*rt,p),Diff(*tt,p))>0) break; 45 ans-=Dist(*tt,*rt),st.erase(tt); 46 } 47 while(lt!=st.begin()) 48 { 49 tt=lt--; 50 if(Cha(Diff(*tt,p),Diff(*lt,p))>0) break; 51 ans-=Dist(*tt,*lt),st.erase(tt); 52 } 53 st.insert(p),lt=rt=st.find(p); 54 lt--,rt++,ans+=Dist(*lt,p)+Dist(*rt,p); 55 } 56 int main() 57 { 58 scanf("%d%d%d",&m,&cx,&cy); 59 st.insert((a){0,0}); 60 st.insert((a){m,0}); 61 st.insert((a){cx,cy}); 62 ans+=Dist((a){0,0},(a){cx,cy}); 63 ans+=Dist((a){m,0},(a){cx,cy}); 64 scanf("%d",&n); 65 for(int i=1;i<=n;i++) 66 scanf("%d%d",&pts[i].xp,&pts[i].yp); 67 scanf("%d",&T); 68 while(T--) 69 { 70 scanf("%d",&t1); 71 if(t1==1) scanf("%d",&t2),era[++p1]=pts[t2],del[t2]=true; 72 else qry[++p2]=p1; 73 } 74 for(int i=1;i<=n;i++) 75 if(!del[i]) Add(pts[i].xp,pts[i].yp); 76 for(int i=p2;i;i--) 77 { 78 while(p1>qry[i]) 79 Add(era[p1].xp,era[p1].yp),p1--; 80 outp[i]=ans; 81 } 82 // for(set<a> ::iterator it=st.begin();it!=st.end();it++) 83 // printf("%d %d\n",it->xp,it->yp); 84 for(int i=1;i<=p2;i++) 85 printf("%.2f\n",outp[i]); 86 return 0; 87 }
浙公网安备 33010602011771号