计算几何随笔

Lifting the Stone

 HDU - 1115 

对于多边形的重心,将多边形三角剖分,求出每一个三角形的重心,和三角形的面积加权以后求平均值,三角形的重心即为三个点的平均值。

循环除6容易爆精度,应该把除数放在循环外,

而且判断double是否为0,应该用==0,<eps是不行的,有负数情况。

复杂度O(n)

#include<bits/stdc++.h>
using namespace std;
const int N=1e6+5000;
typedef double db;
db eps=1e-6;
int n;
struct Point{
  db x,y;
  
  Point(double X=0,double Y=0){x=X;y=Y;}
  Point operator +(Point B){return Point(x+B.x,y+B.y);}
  Point operator -(Point B){return Point(x-B.x,y-B.y);}
  Point operator *(double k){return Point(x*k,y*k);}
  Point operator /(double k){return Point(x/k,y/k);}   

};
typedef Point Vector;

db Cross(Vector A,Vector B){return A.x*B.y-A.y*B.x;}
// Point P[N];
vector<Point>P(N);
db Polygon_area(){
    double ans=0;
    
    for(int i=0;i<n;i++){
      ans+=Cross(P[i],P[(i+1)%n]);
    }

  return ans/2;
}

Point Polygon_center(){
    Point ans(0,0);
    if(Polygon_area()==0)return ans;
    for(int i=0;i<n;i++){
      ans=ans+(P[i]+P[(i+1)%n])*Cross(P[i],P[(i+1)%n]);
    }
    return ans/Polygon_area()/6;
}

int main(){
  int t;cin>>t;
    
  while(t--){
    
    scanf("%d",&n);

    for(int i=0;i<n;i++)scanf("%lf %lf",&P[i].x,&P[i].y);
  
      Point ans=Polygon_center();
  
      printf("%.2f %.2f
",ans.x,ans.y);

  }

  // system("pause");
  
    return 0;
}
View Code

Surround the Trees

 HDU - 1392 

求凸包周长

将多边形的点按x坐标排序,扫描上凸包和下凸包,把偏向右的点剔除,求凸包一个周长。

复杂度O(nlogn)

#include<bits/stdc++.h>
using namespace std;
const int N=1e3+5000;
typedef double db;
db eps=1e-6;
int sgn(db x){
  if(fabs(x)<eps)return 0;
  else return x<0?-1:1;
}
struct Point{
  db x,y;
  
  Point(double X=0,double Y=0){x=X;y=Y;}
  Point operator +(Point B){return Point(x+B.x,y+B.y);}
  Point operator -(Point B){return Point(x-B.x,y-B.y);}
  Point operator *(double k){return Point(x*k,y*k);}
  Point operator /(double k){return Point(x/k,y/k);}   
  // Point operator <(Point B){
  //   return sgn(x-B.x)<0||(sgn(x-B.x)==0&&sgn(y-B.y)<0);
  // }   
  // Point operator ==(Point B){
  //   return sgn(x-B.x)==0&&sgn(y-B.y)==0;
  // }   
  friend bool operator<(Point A,Point B){
    return sgn(A.x-B.x)<0||(sgn(A.x-B.x)==0&&sgn(A.y-B.y)<0);
  }
  friend bool operator==(Point A,Point B){
    return sgn(A.x-B.x)==0&&sgn(A.y-B.y)==0;
  
  }

};
int n;
bool cmp(Point A,Point B){
  return sgn(A.x-B.x)<0||(sgn(A.x-B.x)==0&&sgn(A.y-B.y)<0);
}

Point P[N],ch[N];
typedef Point Vector;

db Cross(Vector A,Vector B){return A.x*B.y-A.y*B.x;}
db Distance(Point A,Point B){return hypot(A.x-B.x,A.y-B.y);}

int convexHull(){
  sort(P,P+n);
  n=unique(P,P+n)-P;
  
  int cnt=0;
  
  for(int i=0;i<n;i++){
    while(cnt>1&&sgn(Cross(ch[cnt-1]-ch[cnt-2],P[i]-ch[cnt-2]))<=0)cnt--;
    ch[cnt++]=P[i];
  }
  
  int j=cnt;
  for(int i=n-2;i>=0;i--){
    while(cnt>j&&sgn(Cross(ch[cnt-1]-ch[cnt-2],P[i]-ch[cnt-2]))<=0)cnt--;
    ch[cnt++]=P[i];
  }
  if(n>1)cnt--;
  return cnt;

}

int main(){

  while(~scanf("%d",&n),n){
    
    for(int i=0;i<n;i++)scanf("%lf %lf",&P[i].x,&P[i].y);
    int cnt=convexHull();

    double ans=0;
    if(cnt==1)ans=0;
    else if(cnt==2)ans=Distance(ch[0],ch[1]);
    else {
      for(int i=0;i<cnt;i++)ans+=Distance(ch[i],ch[(i+1)%cnt]);
    }
      printf("%.2lf
",ans);

  }

  // system("pause");
  
    return 0;
}
View Code

Run

 HDU - 2297 

一些人跑步,给定若干个人的s=s0+vt,求有多少个人在某一个时刻在第一。

建立坐标系,添加一个沿y轴向下,沿x轴向左的向量,对若干个向量求一个半平面交即可。

复杂度O(log(n))

#include<bits/stdc++.h>
using namespace std;
const int N=1e4+5000;
typedef double db;
const db INF=1e16;
const db eps=1e-6;
int sgn(db x){
  if(fabs(x)<eps)return 0;
  else return x<0?-1:1;
}
struct Point{
  db x,y;
  
  Point(double X=0,double Y=0){x=X;y=Y;}
  Point operator +(Point B){return Point(x+B.x,y+B.y);}
  Point operator -(Point B){return Point(x-B.x,y-B.y);}
  Point operator *(double k){return Point(x*k,y*k);}
  Point operator /(double k){return Point(x/k,y/k);}   

};
typedef Point Vector;

struct Line{
  Point p;
  Vector v;
  db ang;
  Line(){}
  Line(Point p,Vector v):p(p),v(v){ang=atan2(v.y,v.x);}
  bool operator <(Line &L){return ang<L.ang;}

};

db Cross(Vector A,Vector B){return A.x*B.y-A.y*B.x;}
db Distance(Point A,Point B){return hypot(A.x-B.x,A.y-B.y);}

//判断点在向量左边,即半平面里面
bool OnLeft(Line L,Point p){return sgn(Cross(L.v,p-L.p))>0;}
Point Cross_point(Line a,Line b){
  Vector u=a.p-b.p;
  db t=Cross(b.v,u)/Cross(a.v,b.v);
  return a.p+a.v*t;
}

vector<Point> HPI(vector<Line>L){

  int n=L.size();
  sort(L.begin(),L.end());//极角排序
  int first,last;
  vector<Point>p(n);
  vector<Line>q(n);
  q[first=last=0]=L[0];
  
  for(int i=1;i<n;i++){
    // cout<<"que  :"<<endl;
    //覆盖队尾
    while(first<last&&!OnLeft(L[i],p[last-1]))last--;
    //覆盖队头
    while(first<last&&!OnLeft(L[i],p[first]))first++;
    q[++last]=L[i];
    // 极角相同保留左边
    if(fabs(Cross(q[last].v,q[last-1].v))<eps){
      last--;
      if(OnLeft(q[last],L[i].p))q[last]=L[i];
    }
      if(first<last)p[last-1]=Cross_point(q[last-1],q[last]);

  }

  // cout<<"out :"<<endl;
  // 剔除无用的队尾元素
  while(first<last&&!OnLeft(q[first],p[last-1]))last--;

  vector<Point>ans;

  if(last-first<=1)return ans;
  p[last]=Cross_point(q[last],q[first]);
  for(int i=first;i<=last;i++)ans.push_back(p[i]);
    return ans;
}

int main(){
    int t;cin>>t;
    while(t--){

      int n;scanf("%d",&n);
      vector<Line>lines;
      lines.push_back(Line(Point(0,0),Vector(0,-1)));
      lines.push_back(Line(Point(0,INF),Vector(-1,0)));
      // cout<<"test :"<<endl;
      for(int i=1;i<=n;i++){
        db s,v;scanf("%lf %lf",&s,&v);
        lines.push_back(Line(Point(0,s),Vector(1,v)));
      }

      vector<Point>ans=HPI(lines);
      // cout<<"test :"<<endl;
      // for(int i=0;i<ans.size();i++){
      //   printf("test :%.2lf %.2lf 
",ans[i].x,ans[i].y);
      // }
      printf("%d
",ans.size()-2);
    }

  // system("pause");
  
    return 0;
}
View Code

Problem G. Interstellar Travel

 HDU - 6325 

给出若干个点,从1到n,每次走的费用为叉积,求字典序最小的,费用最小的方案。

沿着上凸包走必然最小,把点排序出来,求一个上凸包即可。

【排序开始写的有点问题一直wa】

#include<bits/stdc++.h>
using namespace std;
const int N=2e5+5000;
typedef double  db;
db eps=1e-6;
int sgn(db x){
  if(fabs(x)<eps)return 0;
  else return x<0?-1:1;
}
struct Point{
  db x,y;
  int id;
  Point(double X=0,double Y=0){x=X;y=Y;}
  Point operator +(Point B){return Point(x+B.x,y+B.y);}
  Point operator -(Point B){return Point(x-B.x,y-B.y);}
  Point operator *(double k){return Point(x*k,y*k);}
  Point operator /(double k){return Point(x/k,y/k);}   
  friend bool operator<(Point A,Point B){
    return sgn(A.x-B.x)<0||(sgn(A.x-B.x)==0&&sgn(A.y-B.y)<0);
  }
  friend bool operator==(Point A,Point B){
    return sgn(A.x-B.x)==0&&sgn(A.y-B.y)==0;
  }

};
int n;
bool cmp(Point A,Point B){
  // return sgn(A.x-B.x)<0||(sgn(A.x-B.x)==0&&sgn(A.y-B.y)<0);
  return A==B?A.id<B.id:A<B;
}

Point P[N],ch[N];
typedef Point Vector;

db Cross(Vector A,Vector B){return A.x*B.y-A.y*B.x;}
db Distance(Point A,Point B){return hypot(A.x-B.x,A.y-B.y);}
// vector<int>ans;
int convexHull(){
  sort(P,P+n,cmp);
  n=unique(P,P+n)-P;
  
  int cnt=0;

  for(int i=0;i<n;i++){
    while(cnt>=2){
      int d=sgn(Cross(ch[cnt-2]-ch[cnt-1],ch[cnt-1]-P[i]));
      if(d>0||
        (d==0&&P[i].id<ch[cnt-1].id))  cnt--;
      else break;
    }
    ch[cnt++]=P[i];
  }

  for(int i=0;i<cnt;i++){
    printf("%d%c",ch[i].id," 
"[i==cnt-1]);
  }

  // for(int i=0;i<n;i++){
  //   while(cnt>1&&sgn(Cross(ch[cnt-1]-ch[cnt-2],P[i]-ch[cnt-2]))<=0)cnt--;
  //   ch[cnt++]=P[i];
  // }
  
  // int j=cnt;
  // for(int i=n-2;i>=0;i--){
  //   while(cnt>j&&sgn(Cross(ch[cnt-1]-ch[cnt-2],P[i]-ch[cnt-2]))<0)cnt--;
  //   ch[cnt++]=P[i];
  //   ans.push_back(P[i].id);
  // }
  // reverse(ans.begin(),ans.end());
  // ans.push_back(n);

  // if(n>1)cnt--;
  // return cnt;
}

int main(){

  int t;cin>>t;
  while(t--){
    scanf("%d",&n);
    // ans.clear();    
    for(int i=0;i<n;i++)scanf("%lf %lf",&P[i].x,&P[i].y),P[i].id=i+1;
    // cout<<convexHull()<<endl;
      convexHull();

    // for(int i=0;i<ans.size();i++)printf("%d%c",ans[i]," 
"[i==ans.size()-1]);
    //   printf("%.2lf
",ans);

  }

   // system("pause");
  
    return 0;
}
View Code

Strange fuction

 HDU - 2899 

给出一个函数,求最小值。

对函数求导,再求导,发现函数是个先减后增的函数,直接三分。

#include<bits/stdc++.h>
using namespace std;
const int N=2e5+5000;
typedef double  db;
db eps=1e-6;
// int sgn(db x){
//   if(fabs(x)<eps)return 0;
//   else return x<0?-1:1;
// }

// struct Point{
//   db x,y;
//   Point(double X=0,double Y=0){x=X;y=Y;}
//   Point operator +(Point B){return Point(x+B.x,y+B.y);}
//   Point operator -(Point B){return Point(x-B.x,y-B.y);}
//   Point operator *(double k){return Point(x*k,y*k);}
//   Point operator /(double k){return Point(x/k,y/k);}   
//   friend bool operator<(Point A,Point B){
//     return sgn(A.x-B.x)<0||(sgn(A.x-B.x)==0&&sgn(A.y-B.y)<0);
//   }
//   friend bool operator==(Point A,Point B){
//     return sgn(A.x-B.x)==0&&sgn(A.y-B.y)==0;
//   }

// };
// int n;
// bool cmp(Point A,Point B){
//   return sgn(A.x-B.x)<0||(sgn(A.x-B.x)==0&&sgn(A.y-B.y)<0);
// }

// Point P[N],ch[N];
// typedef Point Vector;
// db Cross(Vector A,Vector B){return A.x*B.y-A.y*B.x;}
// db Distance(Point A,Point B){return hypot(A.x-B.x,A.y-B.y);}
db y;

db func(db x){
  return 6*pow(x,7)+8*pow(x,6)+7*pow(x,3)+5*pow(x,2)-y*x;
}

int main(){

  int t;cin>>t;
  while(t--){
    scanf("%lf",&y);
    db ans=1e18,l=0,r=100;
    while(r>l+eps){
      db k=(r-l)/3;
      db mid1=l+k,mid2=r-k;
      if(func(mid1)>=func(mid2)){
          l=mid1;
          ans=min(ans,func(mid2));
      }
      else r=mid2;
    }
  printf("%.4lf
",func(l));

  }

   // system("pause");
    return 0;
}
View Code

模拟退火做法,

每次随机跳一个点,更优则加入答案,否则设概率为$p=e^{frac{E}{T}}$,协调E,只要p<1,即可,随机跳进去。

神奇的是降温系数r=0.98,其他是奇奇怪怪的东西。

#include<bits/stdc++.h>
using namespace std;
const int N=2e5+5000;
typedef double  db;
db eps=1e-8;
db y;
db func(db x){
  return 6*pow(x,7)+8*pow(x,6)+7*pow(x,3)+5*pow(x,2)-y*x;
}
int dir[]={1,-1};
db solve(){
  db T=100;
  db r=0.98;
  db now=50,next,ans=1e18,nowans=func(50);
  
  while(T>eps){

    next=now+dir[rand()%2]*T/10;
    if(next>100||next<0)continue;
    db nextans=func(next);
    db E=nowans-nextans;
    if(E>eps){
      nowans=nextans;
      now=next;
      ans=min(ans,nowans);
    }

    else if(exp(E/T)*RAND_MAX>rand()){
      nowans=nextans;
      now=next;
      ans=min(ans,nowans);
    }

    T*=r;
  }

  return ans;
}

int main(){
  int t;
  cin>>t;
  while(t--){
    scanf("%lf",&y);
  printf("%.4lf
",solve());

  }
   system("pause");
    return 0;
}
View Code

Segment set

 HDU - 1558 

每次判断两个线段是否相交,加入并查集即可。注意端点相交的情况。

#include<bits/stdc++.h>
using namespace std;
#define pb push_back
const int N=1e3+500;
typedef double db;
const db eps=1e-6;
const db PI=acos(-1.0);
int sgn(db x){
  if(fabs(x)<eps)return 0;
  else return x<0?-1:1;
}
struct Point{
  db x,y;
  
  Point(double X=0,double Y=0){x=X;y=Y;}
  Point operator +(Point B){return Point(x+B.x,y+B.y);}
  Point operator -(Point B){return Point(x-B.x,y-B.y);}
  Point operator *(double k){return Point(x*k,y*k);}
  Point operator /(double k){return Point(x/k,y/k);}   
  friend bool operator<(Point A,Point B){
    return sgn(A.x-B.x)<0||(sgn(A.x-B.x)==0&&sgn(A.y-B.y)<0);
  }
  friend bool operator==(Point A,Point B){
    return sgn(A.x-B.x)==0&&sgn(A.y-B.y)==0;
  
  }

};
typedef Point Vector;

struct Line{
    Point p1,p2;//线上的两个点
    Line(){}
    Line(Point p1,Point p2):p1(p1),p2(p2){}
   // Line(Point x,Point y){p1 = x;p2 = y;}
  //  Point(double x,double y):x(x),y(y){}
//根据一个点和倾斜角 angle 确定直线,0<=angle<pi
    Line(Point p,double angle){
        p1 = p;
        if(sgn(angle - PI/2) == 0){p2 = (p1 + Point(0,1));}
        else{p2 = (p1 + Point(1,tan(angle)));}
    }
//ax+by+c=0
    Line(double a,double b,double c){
        if(sgn(a) == 0){
            p1 = Point(0,-c/b);
            p2 = Point(1,-c/b);
        }
        else if(sgn(b) == 0){
            p1 = Point(-c/a,0);
            p2 = Point(-c/a,1);
        }
        else{
            p1 = Point(0,-c/b);
            p2 = Point(1,(-c-a)/b);
        }
    }
};

typedef Line Segment;   //定义线段,两端点是Point p1,p2

double Cross(Vector A,Vector B){return A.x*B.y-A.y*B.x;}
db Distance(Point A,Point B){return hypot(A.x-B.x,A.y-B.y);}

// bool Cross_segment(Point a,Point b,Point c,Point d){//Line1:ab,  Line2:cd
//   double c1=Cross(b-a,c-a),c2=Cross(b-a,d-a);
//   double d1=Cross(d-c,a-c),d2=Cross(d-c,b-c);
//   return sgn(c1)*sgn(c2)<0 && sgn(d1)*sgn(d2)<0;//注意交点是端点的情况不算在内
// }

bool Cross_segment(Point a,Point b,Point c,Point d)
{
    Point p1=a,p2=b,p3=c,p4=d;
 
    if( min(p1.x,p2.x)>max(p3.x,p4.x) || min(p1.y,p2.y)>max(p3.y,p4.y) || min(p3.x,p4.x)>max(p1.x,p2.x) || min(p3.y,p4.y)>max(p1.y,p2.y) )
        return 0;
 
    double k1,k2,k3,k4;
    k1=(p2.x-p1.x)*(p3.y-p1.y) - (p2.y-p1.y)*(p3.x-p1.x);
    k2=(p2.x-p1.x)*(p4.y-p1.y) - (p2.y-p1.y)*(p4.x-p1.x);
    k3=(p4.x-p3.x)*(p1.y-p3.y) - (p4.y-p3.y)*(p1.x-p3.x);
    k4=(p4.x-p3.x)*(p2.y-p3.y) - (p4.y-p3.y)*(p2.x-p3.x);
    return (k1*k2<=eps && k3*k4<=eps);
}

// Point Cross_point(Line a,Line b){
//   Vector u=a.p-b.p;
//   db t=Cross(b.v,u)/Cross(a.v,b.v);
//   return a.p+a.v*t;
// }
Line lines[N];
int fa[N],size[N];
int find(int x){return fa[x]==x?x:fa[x]=find(fa[x]);}
void build(int x,int y){
  int dx=find(x),dy=find(y);
  if(dx!=dy){
    fa[dx]=dy;
    size[dy]+=size[dx];
  }
}
int main(){
    int t;
    cin>>t;
    while(t--){
      int n;scanf("%d",&n);
      int cnt=1;
      for(int i=1;i<=n;i++)fa[i]=i,size[i]=1;
      while(n--){
        char op[10];scanf("%s",&op);
        if(op[0]=='P'){

          scanf("%lf %lf %lf %lf",&lines[cnt].p1.x,&lines[cnt].p1.y,&lines[cnt].p2.x,&lines[cnt].p2.y);
          for(int i=0;i<cnt;i++){
            if(Cross_segment(lines[i].p1,lines[i].p2,lines[cnt].p1,lines[cnt].p2)){ 
              if(find(cnt)==find(i))continue;
              build(find(i),cnt);
            }
          }

          cnt++;
        }

        else {

          int x;scanf("%d",&x);
          printf("%d
",size[find(x)]);
        }

      }
      if(t)printf("
");
    }


    // system("pause");
    return 0;
}
View Code
原文地址:https://www.cnblogs.com/littlerita/p/14053917.html