UOJ Logo DUSH的博客

博客

简单计算几何

2016-01-10 22:42:44 By DUSH

计算几何基础小结

已解决问题:

  判断点是否在线段上

  判断两线段是否相交

  判断点是否在多边形内

  判断线段、折线、多边形是否在多变形内

  判断上述是否在圆内

  计算上述与线段及直线的交点

  凸包

待解决:

  半平面交

#include<ctime>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#include<cmath>
#include<vector>
#include<stack>
#include<queue>

using namespace std;
//================================ struct =================================
struct V{
    double px,py;
    V operator - (const V c){
        V a;
        a.px = px - c.px,a.py = py - c.py;
        return a;
    }
    V operator + (const V c){
        V a;
        a.px = px + c.px,a.py = py + c.py;
        return a;
    }
    V operator * (const double c){
        V a = *this;
        a.px *= c,a.py *= c;
        return a;
    }
    V operator / (const double c){
        V a = *this;
        a.px /= c,a.py /= c;
        return a;
    }
    bool operator <= (const V c){
        //return px <= c.px && py <= c.py;
        return px*px + py*py <= c.px*c.px + c.py*c.py;
    }
    bool operator == (const V c){
        return px == c.px && py == c.py;
    }
    double MO2(){
        return px * px + py * py;
    }
    double DIS(V p1,V p2,int pa){
        if(pa == 1){
            return abs(py - p1.py);
        }
        if(pa == 2){
            return abs(px - p1.px);
        }
        if(pa == 0){
            double k = (p2.py-p1.py)/(p2.px-p1.px),x0,y0;
            x0 = (k*k * p1.px + k * (py - p1.py) + px) / (k*k + 1);
            y0 = k * (px - p1.px) + p1.py;
            return sqrt((px-x0)*(px-x0) + (py-y0)*(py-y0));
        }
    }
    double CP(V c){
        return px*c.py - py*c.px;
    }//cross product
    double DP(V c){
        return px*c.px + py*c.py;
    }//dot product
    bool min(V a,V b){
        if(b <= a) swap(a,b);
        return a <= *this;
    }
    bool max(V a,V b){
        if(a <= b) swap(a,b);
        return *this <= a;
    }
};//vector

struct L{
    V p1,p2;
    int p0;
/*
    switch(p0){
    case 1: beeline
    case 2: segment
    case 3: ray
    }
*/
    int pa;
/*
    switch(pa){
    case 1: parallel-y;
    case 2: parallel-x;
    case 0: k < oo && k != 0;
    }
*/
    double k;
    void getk(){
        if(p1.px == p2.px) pa = 1;
        else if(p1.py == p2.py) pa = 2;
        else pa = 0,k = (p1.py-p2.py) / (p1.px-p2.px);
    }
    bool IP(V q){
        return q.min(p1,p2) && q.max(p1,p2);
    }//in the ploygon of the line
    bool PA(L c){
        if(pa != c.pa) return false;
        if(pa) return true;
        else if(k == c.k) return true;
        else return false;
    }
    bool CR(L c){
        if((*this).PA(c)) return false;
        V v1 = p2 - p1,v2 = c.p1 - p1,v3 = c.p2 - p1;
        double exp = v1.CP(v2) * v1.CP(v3);
        //if(v1.CP(v2) * v1.CP(v3) < 0) return true;
        switch(p0){
        case 1://B_S
            return exp <= 0;
            break;
        case 2://S_S
            return exp < 0;
            break;
        case 3://R_S
            if(exp > 0) return false;
            else if(!exp){
                V v = v1.CP(v2)? c.p2:c.p1;
                if(p1-v <= p2-v) return false;
                else return true;
            }
            else{
                V v4 = c.p1 - p2,v5 = c.p2 - p2;
                double s1 = abs(v2.CP(v3)),s2 = abs(v4.CP(v5));
                return s1 > s2;
            }
            break;
        }
    }//line cross
};//line

struct BL{
    vector<L> E;
};//broken line

struct PL{
    vector<L> E;
};//ploygon

struct O{
    V o;
    double r;
};//circle

//================================ data ===================================
double CP(V a,V b){
    return a.px*b.py - a.py*b.px;
}//cross product

double DP(V a,V b){
    return a.px*b.px + a.py*b.py;
}//dot product

double P_P_DIST(V a,V b){
    return sqrt((a.px-b.px)*(a.px-b.px) + (a.py-b.py)*(a.py-b.py));
}//disdance for points

double P_L_DIST(V a,L b){
    double ret = a.DIS(b.p1,b.p2,b.pa);
    if(b.p0 != 1){
        ret = min(ret,P_P_DIST(a,b.p1));
    }
    if(b.p0 == 2){
        ret = min(ret,P_P_DIST(a,b.p2));
    }
    return ret;
}//least distance for point to line(segment/beeline/ray)

double P_BL_DIST(V p,BL bl){
    double ret = oo;
    for(int i = 0;i < bl.E.size();++ i){
        ret = min(ret,P_L_DIST(p,bl.E[i]));
    }
    return ret;
}//least distance for point to broken line

double P_PL_DIST(V p,PL pl){
    double ret = oo;
    for(int i = 0;i < pl.E.size();++ i){
        ret = min(ret,P_L_DIST(p,pl.E[i]));
    }
    return ret;
}//least distance for point to ploygon

double P_O_DIST(V p,O o){
    return o.r - P_P_DIST(p,o.o);
}//least distance for point to circle

double AREA(PL pl){
    V v0 = pl.E[0].p1;
    double ret = 0;
    for(int i = 0;i < pl.E.size();++ i){
        ret += v0.CP(pl.E[i].p2-pl.E[i].p1);
        v0 = v0 + (pl.E[i].p2-pl.E[i].p1);
    }
    return ret;
}//the measure of area

//================================= judge =================================
bool P_L_ON(V q,L p){
    /*double x1 = p.p1.px,x2 = p.p2.px,y1 = p.p1.py,y2 = p.p2.py,
        x0 = q.px,y0 = q.py;
    if((!p.p0) && (!(x0 <= max(x1,x2) && x0 >= min (x1,x2) && y0 <= max(y1,y2) && y0 >= min(y1,y2))))
      return false;*/
    //if((!p.p0) && (!(q.min(p.p1,p.p2) && q.max(p.p1,p.p2)))) return false;
    if((p.p0 == 2) && (!p.IP(q))) return false;
    int pd = CP(p.p1-q,p.p2-q);
    if(pd != 0) return false;
    pd = DP(p.p1-q,p.p2-q);
    if(pd <= 0) return true;
    else{
        if(p.p0 == 1) return true;
        if(p.p0 == 2) return false;
        if(p.p0 == 3){
            if(p.p1-q <= p.p2-q) return false;
            else return true;
        }
    }
       //if(!CP(q-p.p1,p.p2-p.p1)) return true;
       //else return false;
}//wheather a point is on a line(beeline/segment/ray)

bool S_S_CR(L a,L b){
    return (a.CR(b) && b.CR(a)) || ((P_L_ON(a.p1,b) || P_L_ON(a.p2,b) || P_L_ON(b.p1,a) || P_L_ON(b.p2,a))/* && ((a.p2-a.p1).DP(b.p1-a.p1) * (a.p2-a.p1).DP(b.p2-a.p1) == 0)*/);
}//segments' cross

bool S_L_CR(L a,L b){
    return b.CR(a);
}//segment and beeline's cross

bool S_R_CR(L a,L b){
    return b.CR(a);
}//segment and ray's cross

bool P_PL_IN(V p,PL pl){
    L o;
    o.p1 = p,o.p2.px = p.px-1,o.p2.py = p.py,o.p0 = 3;
    int cnt = 0;
    for(int i = 0;i < pl.E.size();++ i){
        if(P_L_ON(p,pl.E[i])) return true;
        if(pl.E[i].pa != 1){
            if(P_L_ON(pl.E[i].p1,o) || (P_L_ON(pl.E[i].p2,o))){
                V v1 = P_L_ON(pl.E[i].p1,o)? pl.E[i].p1:pl.E[i].p2,
                    v2 = P_L_ON(pl.E[i].p1,o)? pl.E[i].p2:pl.E[i].p1;
                if(v1.py > v2.py) cnt ++;
            }
        }
        else if(S_R_CR(pl.E[i],o)) cnt ++;
    }
    return cnt & 1;
}//point in ploygon

bool S_PL_IN(L p,PL pl){
    if(!(P_PL_IN(p.p1,pl) && P_PL_IN(p.p2,pl))) return false;
    vector<V> v;
    for(int i = 0;i < pl.E.size();++ i){
        if(P_L_ON(p.p1,pl.E[i]) || P_L_ON(p.p2,pl.E[i])){
            if(P_L_ON(p.p1,pl.E[i])) v.push_back(p.p1);
            if(P_L_ON(p.p2,pl.E[i])) v.push_back(p.p2);
        }
        else if(P_L_ON(pl.E[i].p1,p) || P_L_ON(pl.E[i].p2,p)){
            if(P_L_ON(pl.E[i].p1,p)) v.push_back(pl.E[i].p1);
            if(P_L_ON(pl.E[i].p2,p)) v.push_back(pl.E[i].p2);
        }
        else if(p.CR(pl.E[i])) return false;
    }
    for(int i = 0;i < v.size() - 1;++ i){
        if(!P_PL_IN((v[i] + v[i+1])/2,pl)) return false;
    }
    return true;
}//segment in ploygon

bool BL_PL_IN(BL bl,PL pl){
    for(int i = 0;i < bl.E.size();++ i){
        if(!S_PL_IN(bl.E[i],pl)) return false;
    }
    return true;
}//broken line in ploygon

bool PL_PL_IN(PL pl1,PL pl2){
    for(int i = 0;i < pl1.E.size();++ i){
        if(!S_PL_IN(pl1.E[i],pl2)) return false;
    }
    return true;
}//ploygon in ploygon

bool O_PL_IN(O o,PL pl){
    double dis = oo;
    for(int i = 0;i < pl.E.size();++ i){
        dis = min(dis,P_L_DIST(o.o,pl.E[i]));
    }
    return o.r <= dis;
}//circle in ploygon

bool P_O_IN(V p,O o){
    return (p-o.o).MO2() <= o.r * o.r;
}//point in circle

bool S_O_IN(L s,O o){
    return P_O_IN(s.p1,o) && P_O_IN(s.p2,o);
}//segment in circle

bool BL_O_IN(BL bl,O o){
    for(int i = 0;i < bl.E.size();++ i){
        if(!P_O_IN(bl.E[i].p1,o)) return false;
    }
    return P_O_IN(bl.E[bl.E.size()-1].p2,o);
}//broken line in circle

bool PL_O_IN(PL pl,O o){
    for(int i = 0;i < pl.E.size();++ i){
        if(!P_O_IN(pl.E[i].p1,o)) return false;
    }
    return true;
}//ploygon in circle

bool O_O_IN(O o1,O o2){
    if(o1.r > o2.r) return false;
    else return o2.r-o1.r <= P_P_DIST(o1.o,o2.o);
}//circle in circle

//================================= node ==================================
void P_O_ND(V p,O o,double &x,double &y){
    if(p == o.o){
        x = y = oo;
        return;
    }
    L a;
    double lth = P_O_DIST(p,o);
    a.p1 = o.o,a.p2 = p;
    a.getk();
    if(a.pa == 1){
        if(p.py > o.o.py) x = p.px,y = o.o.py + o.r;
        else x = p.px,y = o.o.py - o.r;
    }
    else if(a.pa == 2){
        if(p.px > o.o.px) y = p.py,x = o.o.px + o.r;
        else y = p.py,x = o.o.px - o.r;
    }
    else{
        double xp = o.r / sqrt(a.k*a.k + 1),yp = abs(a.k * xp);
        x = o.o.px,y = o.o.py;
        if(p.px > o.o.px) x += xp;
        else x -= xp;
        if(p.py > o.o.py) y += yp;
        else y -= yp;
    }
}//closest node of point and circle

void sameL_S_S_ND(L a,L b,double &x,double &y){
    V h = a.p1,v1,v2,v3,v4;
    if(a.pa == 1) h.px += 100;
    else if(a.pa == 2) h.py += 100;
    else a.getk(),h.px += 100,h.py += 100 * a.k;
    if(a.p1.MO2() > a.p2.MO2()) swap(a.p1,a.p2);
    if(b.p1.MO2() > b.p2.MO2()) swap(b.p1,b.p2);
    if(a.p2.MO2() > b.p2.MO2()) swap(a,b);
    v1 = a.p1 - h,v2 = a.p2 - h,v3 = b.p1 - h,v4 = b.p2 - h;
    double s1 = abs(CP(v1,v2)),s2 = abs(CP(v3,v4)),s3 = abs(CP(v1,v4));
    if(s1 + s2 < s3) x = y = - oo;//no node
    else if(s1 + s2 > s3) x = y = oo;//oo nodes
    else x = a.p2.px,y = a.p2.py;
}//node of segments which are on the same beeline

void S_SorB_ND(L a,L b,double &x,double &y){
    if(!((b.p0 == 1 && S_L_CR(a,b))||(b.p0 == 2 && S_S_CR(a,b)))){
        x = y = - oo;//no node
        return;
    }
    if(a.pa == 1){
        if(b.pa == 1){
            if(b.p0 == 1) x = y = oo;//oo nodes
            else sameL_S_S_ND(a,b,x,y);
        }
        else{
            x = a.p1.px;
            b.getk();
            y = b.k * (x - b.p1.px) + b.p1.py;
        }
    }
    else if(b.pa == 1){
        x = b.p1.px;
        a.getk();
        y = a.k * (x - a.p1.px) + a.p1.py;
    }
    else if(a.pa == 2){
        if(b.pa == 2){
            if(b.p0 == 1) x = y = oo;//oo nodes
            else sameL_S_S_ND(a,b,x,y);
        }
        else{
            y = a.p1.py;
            b.getk();
            x = 1.0/a.k * (y - b.p1.py) + b.p1.px;
        }
    }
    else if(b.pa == 2){
        y = b.p1.py;
        a.getk();
        x = 1.0/b.k * (y - a.p1.py) + a.p1.px;
    }
    else{
        a.getk(),b.getk();
        if(a.k == b.k) sameL_S_S_ND(a,b,x,y);
        else{
            x = (a.k*a.p1.px - b.k*b.p1.px + b.p1.py - a.p1.px) / (a.k - b.k);
            y = a.k * (x - a.p1.px) + a.p1.py;
        }
    }
}//node of segment and segment(or beeline)

void SorB_BL_ND(L l,BL bl,V v[]){
    int j = 1;
    for(int i = 0;i < bl.E.size();++ i){
        S_SorB_ND(l,bl.E[i],v[j].px,v[j].py);
        if(v[j].px != oo && v[j].px != -oo)
            j ++;
    }    
}//node(s) of segment(or beeline) and broken line

void SorB_PL_ND(L l,PL pl,V v[]){
    int j = 1;
    for(int i = 0;i < pl.E.size();++ i){
        S_SorB_ND(l,pl.E[i],v[j].px,v[j].py);
        if(v[j].px != oo && v[j].px != -oo)
            j ++;
    }
}//node(s) of segment(or beeline) and ploygon

void SorB_O_ND(L l,O o,V v[]){
    if(l.pa == 2 && P_O_IN(l.p1,o) && P_O_IN(l.p2,o)) return;
    //L l0 = l;l0.pa = 1;
    double dis = P_L_DIST(o.o,l);
    V v1,v2;
    if(dis > o.r) return;
    if(l.p0 == 1){
        double h = sqrt(o.r * o.r - dis * dis);
        if(h == 0){
            v[1].px = l.p1.px,v[1].py = o.o.py;
        }
        else{
            v1.px = v2.px = l.p1.px;
            v1.py = v2.py = o.o.py;
            v1.py += h,v2.py -= h;
            if(l.pa == 1){
                v[1] = v1,v[2] = v2;
            }
            else{
                if((P_O_IN(l.p1,o) && P_O_IN(l.p2,o))||(!P_O_IN(l.p1,o) && !P_O_IN(l.p2,o))) return;
                if(P_L_ON(v1,l) && P_L_ON(v2,l)){
                    v[1] = v1,v[2] = v2;
                }
                else if(P_L_ON(v1,l)){
                    v[1] = v1;
                }
                else if(P_L_ON(v2,l)){
                    v[1] = v2;
                }
            }
        }
    }
    else if(l.p0 == 2){
        double w = sqrt(o.r * o.r - dis * dis);
        if(w == 0){
            v[1].py = l.p1.py,v[1].px = o.o.px;
        }
        else{
            v1.py = v2.py = l.p1.py;
            v1.px = v2.px = o.o.px;
            v1.px += w,v2.px -= w;
            if(l.pa == 1){
                v[1] = v1,v[2] = v2;
            }
            else{
                if((P_O_IN(l.p1,o) && P_O_IN(l.p2,o))||(!P_O_IN(l.p1,o) && !P_O_IN(l.p2,o))) return;
                if(P_L_ON(v1,l) && P_L_ON(v2,l)){
                    v[1] = v1,v[2] = v2;
                }
                else if(P_L_ON(v1,l)){
                    v[1] = v1;
                }
                else if(P_L_ON(v2,l)){
                    v[1] = v2;
                }
            }
        }
    }
    else{
        l.getk();
        double A = l.k * l.k + 1,B = 2 * (-l.k * l.p1.px + l.p1.py - o.o.py - o.o.px),C = (-l.k * l.p1.px + l.p1.py - o.o.py) * (-l.k * l.p1.px + l.p1.py - o.o.py) - o.r * o.r,D;// = 4 * ((-l.k * l.p1.px + l.p1.py - o.o.py - o.o.px) * (-l.k * l.p1.px + l.p1.py - o.o.py - o.o.px) - (l.k * l.k + 1) * ((-l.k * l.p1.px + l.p1.py - o.o.py) * (-l.k * l.p1.px + l.p1.py - o.o.py) - o.r * o.r));
        D = B * B - 4 * A * C;
        v1.px = (-B + sqrt(D))/(2*A),v1.py = l.k * (v1.px - l.p1.px) + l.p1.py;
        v2.px = (-B - sqrt(D))/(2*A),v2.py = l.k * (v2.px - l.p1.px) + l.p1.py;
        if(D == 0){
            v[1] = v1;
        }
        else{
            if(l.pa == 1){
                v[1] = v1,v[2] = v2;
            }
            else{
                if((P_O_IN(l.p1,o) && P_O_IN(l.p2,o))||(!P_O_IN(l.p1,o) && !P_O_IN(l.p2,o))) return;
                if(P_L_ON(v1,l) && P_L_ON(v2,l)){
                    v[1] = v1,v[2] = v2;
                }
                else if(P_L_ON(v1,l)){
                    v[1] = v1;
                }
                else if(P_L_ON(v2,l)){
                    v[1] = v2;
                }
            }
        }
    }
}//node(s) of segment(or beeline) and circle

//============================ convex closure =============================
bool cmp(V a,V b){
    return CP(a,b) > 0;
}

stack<V,vector<V> > GCC(V src[],int n){
    V p0 = src[1];
    int o0 = 1;
    for(int i = 2;i <= n;++ i){
        if(p0.py != src[i].py){
            p0 = p0.py < src[i].py? p0:(o0 = i,src[i]);
        }
        else p0 = p0.px < src[i].py? p0:(o0 = i,src[i]);
    }
    V tmp[maxn];
    int j = 0;
    for(int i = 1;i <= n;++ i){
        if(i != o0){
            tmp[++j] = src[i];
        }
    }
    sort(tmp+1,tmp+n,cmp);
    stack<V,vector<V> > S;
    S.push(p0),S.push(tmp[1]),S.push(tmp[2]);
    V ut = tmp[1];//ut: undertop
    for(int i = 3;i < n;++ i){
        while(CP(S.top() - ut,tmp[i] - S.top()) <= 0){
            S.pop(),S.pop();
            V tp = ut;ut = S.top();
            S.push(tp);
        }
        S.push(tmp[i]);
    }
    return S;
}//get convex closure
//================================= main ==================================
const double OO = 1e-8;
const double oo = 1e9;
const int maxn = 10000;



void init(){
}

void work(){
}

int main(){
    init();
    work();
    return 0;
}

求批评指正

字符串哈希C++实现

2015-10-06 20:56:05 By DUSH
最近遇到了好几道字符串比较的题目,感觉直接比较还是挺耗时间的,于是想到了哈希,但是又觉得运算的时候可能会把int数据类型爆掉,所以在讨论哈希之前还是总结一下快速幂。

1、Quickpow

    这是个模板,我总结了一下网络上的各种写法,比这短的反正我是不会的:
int quickpow(int a,int b,int c) { // return (a^b)%c ---------- 把a的类型声明改成long long也许会更好
    int ret = 1;
    a %= c;
    while(b){
        if(b&1) ret = (ret*a) % c;
        b >>= 1;
        a = (a*a) % c;
    }
    return ret;
}
    因为这里用到了分治的思想,所以用递归代码是很容易实现的,但是我并不推荐使用,很容易就TLE了:
int quickpow_2(int a,int b,int c) {
    int ret = 1;
    a %= c;
    ret = (quickpow_2(a,b/2,c) * quickpow_2(a,b/2,c)) %c
    if(b&1){
        ret = (ret * a) % c;
    }
    return ret;
}

2、Hash

    重头戏当然就是哈希函数,这里用到了一个大数和一个大质数用来取模,双关键字比较:
void str_hash(char *s,int &d1,int &d2){ //d1,d2为用于比较的关键字 ---------- sys为进制数
    int i,len = strlen(s);
    for(i = 0;i < len;++ i) {
        d1 = (d1 + (s[i]-'A')*quickpow(sys,i,mod1)) % mod1;
        d2 = (d2 + (s[i]-'A')*quickpow(sys,i,mod2)) % mod2;
    }
}
    利用进制之间的转化公式,很容易得到我们想要的十进制数,也就是代表该字符串的两个数字。

    至于为什么是两个嘛……你会发现只用一个的话,在数据范围比较大的时候容易出现同一个十进制数要表示两个并不相等的字符串,虽然也可以用在哈希表上面挂链表的方式解决,但是对我个人而言是太麻烦了,毕竟咱oier还是不要写太复杂的代码……

3、All

    经过上述讨论,就可以写出完整代码了!!!(Ctrl + C党的福利):
#include<cstdio>
#include<algorithm>
#include<cstring>

using namespace std;
typedef long long LL;
const int sys = 26; //假定字符串只含 A ~ Z,即26进制
const int maxn = 3000; //假定字符串长度小于等于100
const int mod1 = 100000007; //一个常用的取模大数
const int mod2 = 299999321; //一个大质数

char a[maxn];
int ha[maxn][2],n;

int quickpow(LL a,int b,int c) { // return (a^b)%c
    int ret = 1;
    a %= c;
    while(b){
        if(b&1) ret = (ret*a) % c;
        b >>= 1;
        a = (a*a) % c;
    }
    return ret;
}

void str_hash(char *s,int &d1,int &d2){ //d1,d2为用于比较的关键字
    int i,len = strlen(s);
    for(i = 0;i < len;++ i) {
        d1 = (d1 + (s[i]-'A')*quickpow(sys,i,mod1)) % mod1;
        d2 = (d2 + (s[i]-'A')*quickpow(sys,i,mod2)) % mod2;
    }
}

int main(){
    scanf("%d",&n);
    for(int i = 1;i <= n;++ i){
        scanf("%s",a);
        str_hash(a,ha[i][0],ha[i][1]); //读入并处理字符串
//注意,这里并没有保存每个字符串,但是实际运用的时候是需要保存的,毕竟还要输出符合条件的字符串
    }
    for(int i = 1;i <= n;++ i){
        printf("%d %d\n",ha[i][0],ha[i][1]); //输出哈希表
    }
    return 0;
}
DUSH Avatar