2020寒假专题训练：计算几何

Day1 20.01.27

本日任务：

「HNOI2019」鱼

【题目链接】http://59.61.75.5:8018/problem/4382

【题解】

显然要从6个点中选择一个点枚举。

由于其中4个点跟D有关，我们考虑枚举D。

A,B,C为一部分，E,F为另一部分，根据要求，在两个半平面内，考虑分别统计。

先考虑半平面如何划分。注意到题目给出的角度限制都与A有关，于是可以考虑枚举A。

角度问题显然可以转化为极角序，则DA与DE,DF的夹角 $ heta$ 应在$[frac{pi}{2},frac{3pi}{2}]$内。

显然该范围内任意两个到D距离相等的点都可以选择，即相当于有多个团，每个大小为 $x$ 团的贡献为 $frac{x(x-1)}{2}$。

用双指针维护E,F所在的一段区间，枚举A的时候推一下就处理好了。

接下来考虑B,C。根据垂直平分线的判定，AD垂直平分BC。

这部分需要判定两个内容，一个是垂直，另一个是AD与BC中点有交点。注意是线段相交。

垂直、中点预处理一下然后直接找就可以了。垂直的判定可以用点积。

对于每个A，答案为B,C的方案数与E,F的方案数相乘。

综上，效率 $O(n^2log n)$。

注意一下，本题略卡精度，能用整数尽量使用整数运算。

我的代码里只有极角序的atan2用到double类型，其余用的都是long long。

#include<bits/stdc++.h>
const double eps=1e-10;
const double pi=acos(-1);
long long ans;
struct Point
{
    long long x,y;
    Point(){}
    Point(long long _x,long long _y){x=_x;y=_y;}
    inline friend Point operator + ( const Point &u,const Point &v ) { return Point(u.x+v.x,u.y+v.y); }
    inline friend Point operator - ( const Point &u,const Point &v ) { return Point(u.x-v.x,u.y-v.y); }
}A[4010];
inline Point rotate ( const Point &u,const Point &v )
{
    Point p=v-u;
    std::swap(p.x,p.y);p.x=-p.x;
    return u+p;
}
struct Vector
{
    Point p;double k;
    Vector(){}
    Vector(Point _p,double _k){p=_p;k=_k;}
    inline friend bool operator < ( const Vector &u,const Vector &v ) { return u.k<v.k; }
}B[4010];
inline long long get_abs ( long long x ) { return x<0 ? -x : x ; }
struct Line
{
    long long a,b,c,x;
    Line(){}
    Line(long long _a,long long _b,long long _c,long long _x){a=_a;b=_b;c=_c;x=_x;}
    inline friend bool operator < ( const Line &u,const Line &v )
    {
        if ( u.a!=v.a ) return u.a<v.a;
        if ( u.b!=v.b ) return u.b<v.b;
        if ( u.c!=v.c ) return u.c<v.c;
        return u.x<v.x;
    }
};
inline Line calc ( Point u,Point v )
{
    long long dx=v.x-u.x,dy=v.y-u.y;
    long long g=std::__gcd(get_abs(dx),get_abs(dy));
    dx/=g;dy/=g;std::swap(dx,dy);dx=-dx;
    if ( dx==0 ) dy=get_abs(dy);
    if ( dx<0 ) dx=-dx,dy=-dy;
    return Line(dx,dy,u.x*dx+u.y*dy,0);
}
std::vector<Line> X,Y;
long long dis[4010];
int n,num[4010];
std::set<long long> s;
std::map<long long,int> map;
signed main()
{
    scanf("%d",&n);
    for ( int i=1;i<=n;i++ ) scanf("%lld%lld",&A[i].x,&A[i].y),A[i].x<<=1,A[i].y<<=1;
    for ( int i=1;i<=n;i++ ) for ( int j=1;j<=n;j++ ) if ( i!=j )
    {
        Point mid=Point((A[i].x+A[j].x)>>1,(A[i].y+A[j].y)>>1);
        Line tmp=calc(rotate(mid,A[i]),rotate(mid,A[j]));
        X.push_back(Line(tmp.a,tmp.b,tmp.c,mid.x));
        Y.push_back(Line(tmp.a,tmp.b,tmp.c,mid.y));
    }
    std::sort(X.begin(),X.end());
    std::sort(Y.begin(),Y.end());
    for ( int i=1;i<=n;i++ )
    {
        int tot=0,cnt=0;s.clear();map.clear();
        for ( int j=1;j<=n;j++ ) if ( i!=j ) B[++tot]=Vector(A[j],atan2(A[j].y-A[i].y,A[j].x-A[i].x));
        std::sort(B+1,B+n);
        for ( int j=1;j<n;j++ ) B[j+n-1]=Vector(B[j].p,B[j].k+2*pi);
        for ( int j=1;j<n;j++ ) s.insert(dis[j]=(A[i].x-B[j].p.x)*(A[i].x-B[j].p.x)+(A[i].y-B[j].p.y)*(A[i].y-B[j].p.y));
        for ( long long x:s ) map[x]=++cnt;
        for ( int j=1;j<n;j++ ) dis[j+n-1]=dis[j]=map[dis[j]];
        for ( int j=1;j<=cnt;j++ ) num[j]=0;
        long long ans1=0,ans2=0;
        for ( int j=1,l=0,r=0;j<n;j++ )
        {
            while ( r<2*n-1 and B[j].k+1.5*pi-B[r+1].k>eps ) r++,ans1+=num[dis[r]]++;
            while ( l<2*n-1 and B[j].k+0.5*pi-B[l+1].k>-eps ) l++,ans1-=--num[dis[l]];
            Line t=calc(A[i],B[j].p),tl=t,tr=t;
            if ( A[i].x!=B[j].p.x )
                tl.x=std::min(A[i].x,B[j].p.x),tr.x=std::max(A[i].x,B[j].p.x),
                ans2=std::lower_bound(X.begin(),X.end(),tr)-std::upper_bound(X.begin(),X.end(),tl);
            else
                tl.x=std::min(A[i].y,B[j].p.y),tr.x=std::max(A[i].y,B[j].p.y),
                ans2=std::lower_bound(Y.begin(),Y.end(),tr)-std::upper_bound(Y.begin(),Y.end(),tl);
            ans+=ans1*std::max(ans2,0LL);
        }
    }
    return !printf("%lld
",ans<<1);
}

---

---

Day2 20.01.28

 本日任务：

算法一：最小圆覆盖

算法二：自适应simpson积分法

最小圆覆盖

最小圆覆盖通常采用的算法为随机增量法，除此之外，模拟退火法也是可以应用在最小圆覆盖中的。

例题：「AHOI2012」信号塔

【题目链接】http://59.61.75.5:8018/problem/2167

【题解】

方法一：随机增量法

考虑将每个点依次加入。

不难证明，新加入的一个点，若不在原来的点的最小圆覆盖内，则一定在加入后的最小圆覆盖上。此时确定了一个点。

接下来贪心地选取。假设加入的点为第 $i$ 个点，那么先确定一个以 $1,i$ 为直径的圆。

枚举前 $i-1$ 个点，若点 $j$ 不在圆内，则 $j$ 必定在前 $j$ 个点和 $i$ 的最小圆覆盖上。此时确定了两个点。

同理，枚举前 $j-1$ 个点，若不在圆内，则令其在圆上。此时确定了三个点。

根据三点确定一个圆，此时我们已经得到前 $i$ 个点的最小圆覆盖。递推即可。

但注意到一个问题，此时需要三重循环，效率是 $O(n^3)$ 的，考虑优化。

由于三点确定一个圆，则一个点在圆上的概率为 $frac{3}{n}$，于是考虑将给出的点随机排序，这样效率的期望就是 $O(n)$ 的了。

#include<bits/stdc++.h>
#define clr(a) memset(a,0,sizeof(a))
#define print(x) std::cerr<<#x<<'='<<x<<'
'
const double eps=1e-6;
int n;
struct point { double x,y; }p[500000+10];
struct circle { point p;double r; }ans;
inline double dist ( point p1,point p2 ) { return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)); }
inline circle Circle ( point p1,point p2,point p3 )
{
    circle tmp;
    double a=dist(p1,p2),b=dist(p1,p3),c=dist(p2,p3);
    double a1=p2.x-p1.x,b1=p2.y-p1.y,c1=(a1*a1+b1*b1)/2;
    double a2=p3.x-p1.x,b2=p3.y-p1.y,c2=(a2*a2+b2*b2)/2;
    double d=a1*b2-a2*b1;
    tmp.p=(point){p1.x+(c1*b2-c2*b1)/d,p1.y+(a1*c2-a2*c1)/d};
    tmp.r=dist(tmp.p,p1);
    return tmp;

}
int main()
{
    scanf("%d",&n);
    for ( int i=1;i<=n;i++ ) scanf("%lf%lf",&p[i].x,&p[i].y);
    std::random_shuffle(p+1,p+n+1);ans.p=p[1];ans.r=0;
    for ( int i=2;i<=n;i++ ) if ( dist(p[i],ans.p)>ans.r+eps )
    {
        ans.p=p[i];ans.r=0;
        for ( int j=1;j<i;j++ ) if ( dist(p[j],ans.p)>ans.r+eps )
        {
            ans.p=(point){(p[i].x+p[j].x)/2,(p[i].y+p[j].y)/2};
            ans.r=dist(p[i],p[j])/2;
            for ( int k=1;k<j;k++ ) if ( dist(p[k],ans.p)>ans.r+eps ) ans=Circle(p[i],p[j],p[k]);
        }
    }
    return !printf("%.2lf %.2lf %.2lf
",ans.p.x,ans.p.y,ans.r);
}

方法二：模拟退火法

其实模拟退火的算法并不适用于较大数据范围的题，但由于网上有一大部分题解都有提到该算法，因此在此讲解。

由于此算法的不稳定性，暂时没有本题的代码，这里只给出实现方式。

首先，要明白模拟退火算法在做什么。

假设我们要求的值是一个函数。模拟退火算法的关键就在于如何找到函数的最值。

对于具有单调性的函数，可以采用二分或三分的方法解决，但对于一个普通的函数，就没有好的解决方法了。

我们依然考虑随即增量，即向左或向右偏移一点。若偏移后的函数值比原函数值大，则向偏移的位置移动。

一直移动下去，总会到达某处，使得左右两边的函数值都小于此处的函数值。此时有两种情况：

1. 该函数为单峰函数，此处即为最大值。

2. 此处为该函数的局部最大值。

到这里就是爬山算法，可以解决单峰函数，也就是情况1的问题。

接下来考虑如何处理情况2。我们唯一的方式就是跳出这个点，选择一个可能更劣的点。

现在考虑跳出的情况。显然我们并不知道我们所处的位置是不是峰值，不能直接跳出，因此考虑随机跳出。

模拟退火算法的核心就在于随机值如何确定。这里采用的是模仿金属退火过程的温度变化进行随机判定，具体实现可以去网上查找其他题解。

现在回到本题。要求的值是一个关于位置的函数，函数值为当前位置到所有点的距离的最大值。进行模拟退火算法即可。

但注意到对精确度的要求，所以需要一些奇技淫巧小技巧：

1. 可以逐位考虑，在前一位确定的情况下提高精度，进行新一轮的退火。

2. 显然不在凸包上的点不可能在最小圆覆盖上，所以可以先求一遍凸包，然后大致定在凸包的中心处开始退火。

再次重申：由于数据范围较大，本题并不适合用模拟退火法实现，这里只是提供一个思路而已。

自适应simpson积分法

 这又是一个计算几何中的冷门算法。。。上次考查此算法为NOI2005的月下柠檬树，在下方的例题中有讲到。

现在

考虑这么一个问题：求任意曲线与 $x$ 轴，直线 $x=l$，直线 $x=r$ 包围的图形的面积。

曲边梯形的面积问题，通常用定积分可以解决。但任意曲线就出现了一个问题：如何求解定积分？

显然对于任意函数，求解定积分不是那么容易，我们考虑用近似的方式解决问题。

假设当前处理 $[l,r]$ 部分的曲边梯形面积。由于是曲边梯形，考虑用二次函数近似。

设二次函数：$f(x)=ax^2+bx+c(x in [l,r])$。

积分得：
$$egin{aligned} int_{l}^{r}f(x)dx
&=int_{l}^{r}(ax^2+bx+c)dx
\ &=frac{1}{3}ax^3+frac{1}{2}bx^2+cx+CBig|_l^r
\&=(frac{1}{3}ar^3+frac{1}{2}br^2+cr+C)-(frac{1}{3}al^3+frac{1}{2}bl^2+cl+C)
\&=frac{1}{6}(r-l)(2a(r^2+lr+l^2)+3b(l+r)+6c)
\&=frac{1}{6}(r-l)(2al^2+2ar^2+2alr+3bl+3br+6c)
\&=frac{1}{6}(r-l)((al^2+bl+c)+(ar^2+br+c)+(al^2+2alr+ar^2+2bl+2br+4c))
\&=frac{1}{6}(r-l)(f(l)+f(r)+4(a(frac{l+r}{2})^2+b(frac{l+r}{2})+c))
\&=frac{1}{6}(r-l)(f(l)+f(r)+4f(frac{l+r}{2}))
end{aligned}$$
于是我们就解决了 $[l,r]$ 部分的曲边梯形面积。

下面考虑近似的边界。我们在处理 $[l,r]$ 时，可以先用上面给出的公式算出 $[l,r],[l,frac{l+r}{2}],[frac{l+r}{2},r]$ 三部分曲边梯形的近似面积。若左右两边的曲边梯形近似面积和总的曲边梯形近似面积差值小于 $eps$，则返回，否则两边递归求解。这就是“自适应”的方法。一般代码实现会涉及到一些常数以及 $eps$ 的精度在传递过程中需不断提高的问题，具体参考代码。

#include<bits/stdc++.h>
double a,b,c,d;
inline double f ( double x ) { return (c*x+d)/(a*x+b); }
inline double simpson ( double a,double b )
{
    double c=(a+b)*0.5;
    return (f(a)+f(b)+4.0*f(c))*(b-a)/6.0;
}
inline double ars ( double a,double b,double eps )
{
    double c=(a+b)*0.5,mid=simpson(a,b),l=simpson(a,c),r=simpson(c,b);
    if ( fabs(l+r-mid)<=15*eps ) return l+r+(l+r-mid)/15.0;
    return ars(a,c,eps*0.5)+ars(c,b,eps*0.5); 
}
int main()
{
    double l,r;
    scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&l,&r);
    printf("%.6lf
",ars(l,r,1e-6));
    return 0;
}

simpson的例题较少，一般有两种考查方式：

1. 结合数学分析进行考查，如luogu P4526：https://www.luogu.com.cn/problem/P4526

本题主要考察指数函数的性质与simpson模板的应用，代码如下：

#include<bits/stdc++.h>
long double a;
inline long double f ( long double x ) { return pow(x,a/x-x); }
inline long double simpson ( long double a,long double b )
{
    long double c=(a+b)*0.5;
    return (f(a)+f(b)+4.0*f(c))*(b-a)/6.0;
}
inline long double ars ( long double a,long double b,long double eps )
{
    long double c=(a+b)*0.5,mid=simpson(a,b),l=simpson(a,c),r=simpson(c,b);
    if ( fabs(l+r-mid)<=15*eps ) return l+r+(l+r-mid)/15.0;
    return ars(a,c,eps*0.5)+ars(c,b,eps*0.5); 
}
int main()
{
    scanf("%LF",&a);
    if ( a<0 ) puts("orz");
    else printf("%.5LF
",ars(1e-7,15,1e-7));
    return 0;
}

2. 结合计算几何其他内容（如相交圆）进行考查，如NOI2005月下柠檬树：https://www.luogu.com.cn/problem/P4207

动手画画应该就明白了吧。。。这里只给出网上的一张经典图片：

接下来的部分交给数学。什么相似解几之类的式子推一推就完事了。

#include<bits/stdc++.h>
const double eps=1e-6;
typedef std::pair<double,double> P;
#define fst first
#define snd second
#define G(x,y) std::make_pair(x,y)
struct Line
{
    P s,t;
    double k,b;
    Line(){}
    Line(P _s,P _t){s=_s;t=_t;k=(s.snd-t.snd)/(s.fst-t.fst);b=s.snd-s.fst*k;}
    inline double f ( double x ) { return k*x+b; }
}L[10000];
int n,m;double H[10000],alpha;P C[10000];
double sina,cosa,tana;
inline void add ( P a,P b )
{
    double sina=(a.snd-b.snd)/(b.fst-a.fst),cosa=sqrt(1-pow(sina,2));
    L[++m]=Line(G(a.fst+a.snd*sina,a.snd*cosa),G(b.fst+b.snd*sina,b.snd*cosa));
}
inline double f ( double x )
{
    double res=0;
    for ( int i=1;i<=m;i++ ) if ( L[i].s.fst<=x and x<=L[i].t.fst ) res=std::max(res,L[i].f(x));
    for ( int i=1;i<=n;i++ ) if ( C[i].fst-C[i].snd<=x and x<=C[i].fst+C[i].snd ) res=std::max(res,sqrt(pow(C[i].snd,2)-pow(x-C[i].fst,2)));
    return res;
}
inline double simpson ( double a,double b )
{
    long double c=(a+b)*0.5;
    return (f(a)+f(b)+4.0*f(c))*(b-a)/6.0;
}
inline double ars ( double a,double b,double eps )
{
    long double c=(a+b)*0.5,mid=simpson(a,b),l=simpson(a,c),r=simpson(c,b);
    if ( fabs(l+r-mid)<=15*eps ) return l+r+(l+r-mid)/15.0;
    return ars(a,c,eps*0.5)+ars(c,b,eps*0.5); 
}
inline int cmp ( double x ) { return (fabs(x)<eps) ? 0 : ( x>0 ? 1 : -1 ) ; }
int main()
{
    scanf("%d%lf",&n,&alpha);
    for ( int i=1;i<=n+1;i++ ) scanf("%lf",&H[i]),H[i]+=H[i-1];
    for ( int i=1;i<=n;i++ ) scanf("%lf",&C[i].snd);
    alpha=tan(alpha);
    P p=G(H[n+1]/alpha,0);
    double Lmin=2e9,Rmax=0;
    Rmax=std::max(Rmax,p.fst);
    C[n].fst=H[n]/alpha;
    double x=C[n].fst;
    double r=C[n].snd;
    Lmin=std::min(Lmin,x-r);
    Rmax=std::max(Rmax,x+r);
    if ( x+r<p.fst )
    {
        double l=pow(r,2)/(p.fst-x);
        L[++m]=Line(G(x+l,sqrt(pow(r,2)-pow(l,2))),p);
    }
    for ( int i=n-1;i;i-- )
    {
        C[i].fst=H[i]/alpha;x=C[i].fst;r=C[i].snd;
        Lmin=std::min(Lmin,x-r);Rmax=std::max(Rmax,x+r);
        if ( cmp(C[i+1].fst-x-fabs(C[i+1].snd-r))>0 ) add(C[i],C[i+1]);
    }
    return !printf("%.2lf
",2*ars(Lmin,Rmax,eps));
}

hdu上还有一道题是求椭圆在 $[l,r]$ 范围的面积，直接套板子即可。

---

---

 Day3 20.01.29

本日任务：

「SDOI2016」平凡的骰子

【题目链接】http://59.61.75.5:8018/problem/2688

【题解】

Step 1：求重心

根据重心的定义式有：

$$overrightarrow{x_c}=frac{sum m_poverrightarrow{x_p}}{sum m_p}$$

对于三角形的重心很好计算，同理不难得到四面体的重心公式。

对于多边形的重心，通常采用三角形分割的方式；同理，对于多面体的重心，通常采用四面体分割的方式。

我们一般在多面体内随便找一个点（常用所有点坐标的平均数），然后在每个面进行三角形分割，将三角形与找的点连成四面体。

下面给出四面体的体积公式：

$$V=frac{1}{6}lvert overrightarrow{a} imes overrightarrow{b} cdot overrightarrow{c}  vert$$

其中 $overrightarrow{a},overrightarrow{b},overrightarrow{c}$ 为从四面体同一个顶点出发的三条边的向量。

至此就可以求出多面体的重心。

Step 2：算面积

根据题目给出的提示，我们只需进行三角形分割及二面角计算即可。

不难发现，面积即为内角和减去 $(n-2) imes pi$，其中 $n$ 为点数。

于是对于每个二面角，先求出平面上两个向量的叉积求出平面面的法向量，然后直接点积，求反余弦值即可。

效率 $O(nf)$，轻松过掉本题。

PS：本题数据范围有误，应为 $f leq 100$。

#include<bits/stdc++.h>
const double pi=acos(-1);
int n,f,cnt[200];
struct Vector
{
    double x,y,z;
    Vector(){x=y=z=0;}
    Vector(double _x,double _y,double _z){x=_x;y=_y;z=_z;}
    inline friend Vector operator + ( Vector u,Vector v ) { return Vector(u.x+v.x,u.y+v.y,u.z+v.z); }
    inline friend Vector operator - ( Vector u,Vector v ) { return Vector(u.x-v.x,u.y-v.y,u.z-v.z); }
    inline friend Vector operator * ( Vector u,double v ) { return Vector(u.x*v,u.y*v,u.z*v); }
    inline friend Vector operator / ( Vector u,double v ) { return Vector(u.x/v,u.y/v,u.z/v); }
    inline friend double operator * ( Vector u,Vector v ) { return u.x*v.x+u.y*v.y+u.z*v.z; }
    inline friend Vector operator ^ ( Vector u,Vector v ) { return Vector(u.y*v.z-u.z*v.y,u.z*v.x-u.x*v.z,u.x*v.y-u.y*v.x); }
    inline double Len ( void ) { return sqrt(x*x+y*y+z*z); }
    inline Vector unit ( void ) { return *this/Len(); }
}p[200],g[200][200];
inline double V ( Vector u,Vector v,Vector w ) { return ((u^v)*w)/6; }
inline int get ( int x,int v ) { return (x-1)%v+1; }
inline double calc ( Vector G,int t )
{
    double ans=0;
    for ( int i=1;i<=cnt[t];i++ )
    {
        Vector u=g[t][i]-G,v=g[t][get(i+1,cnt[t])]-G,w=g[t][get(i+2,cnt[t])]-G;
        ans+=acos(((u^v).unit())*((w^v).unit()));
    }
    return ans-pi*(cnt[t]-2);
}
signed main()
{
    scanf("%d%d",&n,&f);
    for ( int i=1;i<=n;i++ ) scanf("%lf%lf%lf",&p[i].x,&p[i].y,&p[i].z);
    for ( int i=1;i<=f;i++ )
    {
        scanf("%d",&cnt[i]);
        for ( int j=1,x;j<=cnt[i];j++ ) scanf("%d",&x),g[i][j]=p[x];
    }
    Vector c,d;
    double sum=0;
    for ( int i=1;i<=n;i++ ) c=c+p[i];
    c=c/n;
    for ( int i=1;i<=f;i++ ) for ( int j=2;j<cnt[i];j++ )
    {
        Vector u=g[i][1]-c,v=g[i][j]-c,w=g[i][j+1]-c;
        double res=V(u,v,w);d=d+((u+v+w)/4*res);sum+=res;
    }
    Vector G=c+(d/sum);
    for ( int i=1;i<=f;i++ ) printf("%.7lf
",calc(G,i)/(4*pi));
    return 0;
}

---

---

Day4 20.01.30

本日任务：

算法一：扫描线

算法二：半平面交

扫描线

扫描线是处理二维平面问题的一个重要工具，在数据结构中也有重要应用，常与线段树或平衡树相结合。

例题1：矩形面积并

【题目链接】https://www.luogu.com.cn/problem/T110664

【题解】

先放张图（没错就是样例的图）：

  解决面积并一类问题，在数据范围较小的时候有一个很漂亮的做法：容斥。由于是求并的面积，可以容斥计算。

但这题的范围高达 $10^5$，显然需要一种优秀的做法。

考虑我们小学怎么算这种图形的面积：割补法。我们考虑割补。由于是矩形，因此边界都在格子上，比较好算。

我们考虑竖着切割。继续放图：

  切割完，我们需要求的就变成每个小矩形的面积。

考虑先处理在哪里切割。切割的地方显然应该是矩形左右两边。因此我们需要把每个矩形的两边离散化后存下来（本题数据较小可以不用离散化）。

每个小矩形的面积应为相邻两条直线的距离乘上直线上被覆盖的长度。显然可以用数据结构维护被覆盖的区域。

因此将每个矩形拆解为在 $x_1$ 处插入区间 $[y_1,y_2]$ 与在 $x_2$ 处删除区间 $[y_1,y_2]$ 两个操作，从左往右扫一遍，用线段树维护即可。

以上即为扫描线。除了计算几何，一些偏序问题也可以转化到二维平面上用扫描线的思想解决。

#include<bits/stdc++.h>
inline int read ( void )
{
    int x=0;char ch;bool f=true;
    while ( !isdigit(ch=getchar()) ) if ( ch=='-' ) f=false;
    for ( x=ch^48;isdigit(ch=getchar()); ) x=(x<<1)+(x<<3)+(ch^48);
    return f ? x : -x ;
}
const int maxn=200000+10;
int n;long long ans;
std::vector<std::pair<int,int>> I[maxn],D[maxn];
#define ls (k<<1)
#define rs (k<<1|1)
#define mid ((l+r)>>1)
struct tree { int res,sum; } t[maxn<<2];
inline void pushup ( int k,int l,int r )
{
    if ( t[k].res ) t[k].sum=r-l+1;
    else if ( l==r ) t[k].sum=0;
    else t[k].sum=t[ls].sum+t[rs].sum;
}
inline void update ( int k,int l,int r,int ql,int qr,int f )
{
    if ( ql<=l and r<=qr ) { t[k].res+=f;pushup(k,l,r);return; }
    if ( ql<=mid ) update(ls,l,mid,ql,qr,f);
    if ( qr>mid ) update(rs,mid+1,r,ql,qr,f);
    pushup(k,l,r);
}
#undef ls
#undef rs
#undef mid
signed main()
{
//    freopen(".in","r",stdin);
//    freopen(".out","w",stdout);
    n=read();
    for ( int i=1;i<=n;i++ )
    {
        int X1=read()+100000,X2=read()+100000,Y1=read(),Y2=read();
        I[X1].push_back(std::make_pair(Y1,Y2-1));
        D[X2].push_back(std::make_pair(Y1,Y2-1));
    }
    for ( int lst=0,x=0;x<=200000;x++ ) if ( (int)I[x].size() or (int)D[x].size() )
    {
        ans+=1LL*(x-lst)*t[1].sum;
        for ( auto p:I[x] ) update(1,-100000,100000,p.first,p.second,1);
        for ( auto p:D[x] ) update(1,-100000,100000,p.first,p.second,-1);
        lst=x;
    }
    return !printf("%lld
",ans);
}

例题2：「DTOJ4428」团

【题目链接】http://59.61.75.5:8018/problem/4428

【题解】

这是一个很经典的问题。不难发现题目要求的即为一个转了 $45$ 度的正方形内点数的最大值。

考虑旋转坐标系，转化为正常的正方形。旋转时注意处理 $R$ 的变化。

现在来解决这个问题。根据上一题的经验，采用扫描线算法。

由于本题维护的是一个正方形，因此扫描线应该维护的是一段区间。

另外一个方向同理用线段树维护最值即可。注意应将点的权值转化为区间加，答案统计线段树的每个叶子节点会被多少个点对应的区间覆盖到。

#include<bits/stdc++.h>
const int Lmin=-200000000,Rmax=100000000,maxn=300000+10;
std::pair<int,int> p[maxn];
struct tree { int max,tag,ls,rs; } t[maxn*100];
int n,R,ans,tot,Root;
#define mid ((l+r)>>1)
#define ls t[k].ls
#define rs t[k].rs
inline void pushdown ( int k )
{
    if ( !ls ) ls=++tot;
    if ( !rs ) rs=++tot;
    t[ls].max+=t[k].tag;t[ls].tag+=t[k].tag;
    t[rs].max+=t[k].tag;t[rs].tag+=t[k].tag;
    t[k].tag=0;
}
inline void update ( int &k,int l,int r,int ql,int qr,int val )
{
    if ( !k ) k=++tot;
    if ( ql<=l and r<=qr ) { t[k].max+=val;t[k].tag+=val;return; }
    if ( t[k].tag ) pushdown(k);
    if ( ql<=mid ) update(ls,l,mid,ql,qr,val);
    if ( qr>mid ) update(rs,mid+1,r,ql,qr,val);
    t[k].max=std::max(t[ls].max,t[rs].max);
}
signed main()
{
    scanf("%d%d",&n,&R);
    for ( int i=1,x,y;i<=n;i++ ) scanf("%d%d",&x,&y),p[i].first=x+y,p[i].second=x-y;
    std::sort(p+1,p+n+1);
    for ( int l=1,r=0;l<=n;l++ )
    {
        while ( r<n and p[r+1].first-p[l].first<=R ) r++,update(Root,Lmin,Rmax,p[r].second-R,p[r].second,1);
        ans=std::max(ans,t[Root].max);update(Root,Lmin,Rmax,p[l].second-R,p[l].second,-1);
    }
    return !printf("%d
",ans);
}

半平面交

例题1：「CQOI2006」凸多边形

【题目链接】http://59.61.75.5:8018/problem/1687

【题解】

先不急着讨论本题。先讲讲半平面交。

定义：

半平面：一条直线将一个平面分为两个半平面。这里我们统一规定，直线为有向直线（点向式，即用一个点和一个向量表示直线），我们需要的为直线的左侧（沿直线方向的左侧）。

多边形的核：多边形中的一块区域，该区域中的任意一个点都能直接看到多边形的任意一条边（视线与边相切也可以），即每条边上任意一点与其连线与其它边均无交点。

不难发现，将多边形的边逆时针遍历一遍，每条边所在直线（注意是有向的）左侧的半平面的交即为多边形的核。

现在回到本题。由于每个多边形都是凸多边形，凸多边形的核即其本身。因此交集即为核的交集。核的交集即半平面交。

将每个多边形拆成直线，执行半平面交算法即可。

下面介绍半平面交算法：

Step 0：补全边界

由于处理问题，可以先补全边界，防止出现无限大的半平面交（虽然数据并不会，但有些题目需要）。

Step 1：排序

将每条直线按方向（不是倾斜角）的极角序排序，同方向（注意不是平行）的直线取最左侧那条即可。

Step 2：加入直线

考虑每次新加入一条直线，如何维护已有直线及交点。

 加入红色线时，与之前最后一条直线产生交点。

  加入红色线时，与倒数第二条线有交点。此时应在加入前先判断，删除最后几条不必要的直线。

  加入新直线时，与首尾直线均有交点，因此两边都需要维护判定。

 与第二条直线有交点，从而需要删去第一条直线。

具体的判定都是判交点在直线的左右侧。可以用叉积的正负性判断。

综上所述，我们需要维护直线与交点，同时支持以下操作：

1. 删除最前端

2. 删除最后端

3. 从最后端加入

于是就可以用双端队列维护。具体实现参考代码：

#include<algorithm>
#include<cstdio>
#include<cmath>
const double eps=1e-6,lim=1e3;
inline int dcmp ( double x ) { return fabs(x)<eps ? 0 : ( x<0 ? -1 : 1 ) ; }
struct Vector
{
    double x,y;
    Vector(double _x=0,double _y=0){x=_x;y=_y;}
    inline friend Vector operator + ( Vector u,Vector v ) { return Vector(u.x+v.x,u.y+v.y); }
    inline friend Vector operator - ( Vector u,Vector v ) { return Vector(u.x-v.x,u.y-v.y); }
    inline friend Vector operator * ( Vector u,double v ) { return Vector(u.x*v,u.y*v); }
    inline friend Vector operator / ( Vector u,double v ) { return Vector(u.x/v,u.y/v); }
    inline friend double operator * ( Vector u,Vector v ) { return u.x*v.x+u.y*v.y; }
    inline friend double operator ^ ( Vector u,Vector v ) { return u.x*v.y-u.y*v.x; }
    inline double calc_angle ( void ) { return atan2(y,x); }
}Q2[30000],p[30000];
struct Line
{
    Vector s,t;double A;
    Line(Vector _s=Vector(),Vector _t=Vector()){s=_s;t=_t-_s;A=t.calc_angle();}
    inline friend bool operator < ( Line u,Line v )
    {
        if ( dcmp(u.A-v.A)!=0 ) return dcmp(u.A-v.A)==-1;
        return dcmp(u.t^(v.t+v.s-u.s))==-1;
    }
}L[30000],Q[30000];
inline bool PX ( Line u,Line v ) { return dcmp(u.t^v.t)==0; }
inline Vector JD ( Line u,Line v ) { return u.s+u.t*((v.t^(u.s-v.s))/(u.t^v.t)); }
inline bool OnRight ( Line u,Vector v ) { return dcmp(u.t^(v-u.s))<0; }
int n,m,tot;
signed main()
{
    scanf("%d",&tot);
    for ( int i=1;i<=tot;i++ )
    {
        int cnt,X[100]={0},Y[100]={0};scanf("%d",&cnt);
        for ( int i=1;i<=cnt;i++ ) scanf("%d%d",&X[i],&Y[i]);
        for ( int i=2;i<=cnt;i++ ) L[++n]=Line(Vector(X[i-1],Y[i-1]),Vector(X[i],Y[i]));
        L[++n]=Line(Vector(X[cnt],Y[cnt]),Vector(X[1],Y[1]));
    }
    Vector p1=Vector(-lim,-lim),p2=Vector(lim,-lim),p3=Vector(lim,lim),p4=Vector(-lim,lim);
    L[++n]=Line(p1,p2);L[++n]=Line(p2,p3);L[++n]=Line(p3,p4);L[++n]=Line(p4,p1);
    std::sort(L+1,L+n+1);
    int l=0,r=0;Q[0]=L[1];
    for ( int i=2;i<=n;i++ ) if ( dcmp(L[i].A-L[i-1].A) )
    {
        if ( l<r and ( PX(Q[l],Q[l+1]) or PX(Q[r],Q[r-1]) ) ) return !puts("0.000");
        while ( l<r and OnRight(L[i],Q2[r-1]) ) r--;
        while ( l<r and OnRight(L[i],Q2[l]) ) l++;
        Q[++r]=L[i];
        if ( l<r ) Q2[r-1]=JD(Q[r],Q[r-1]);
    }
    while ( l<r and OnRight(Q[l],Q2[r-1]) ) r--;
    while ( l<r and OnRight(Q[r],Q2[l]) ) l++;
    if ( r-l<=1 ) return !puts("0.000");
    Q2[r]=JD(Q[l],Q[r]);
    for ( int i=l;i<=r;i++ ) p[++m]=Q2[i];
    double ans=p[m]^p[1];
    for ( int i=2;i<=m;i++ ) ans+=p[i-1]^p[i];
    return !printf("%.3lf
",ans*0.5);
}

例题2：「SCOI2015」小凸想跑步

【题目链接】http://59.61.75.5:8018/problem/2348

【题解】

考虑选取点 $P$ 使得面积最小的要求。

直接用向量叉积算面积，可以推出来一个式子，形式为 $Ax+By+C<0$。

看到这个式子应该非常熟悉，想起高中解几学的直线的一般式。显然直线的一侧即为半平面。

因此每个条件当成一个半平面，直接半平面交即可。

注意：选取点需要在多边形内，因此将多边形每条边也作为半平面的限制。

#include<bits/stdc++.h>
const double eps=1e-6;
inline int dcmp ( double x ) { return fabs(x)<eps ? 0 : ( x<0 ? -1 : 1 ) ; }
struct Vector
{
    double x,y;
    Vector(double _x=0,double _y=0){x=_x;y=_y;}
    inline friend Vector operator + ( Vector u,Vector v ) { return Vector(u.x+v.x,u.y+v.y); }
    inline friend Vector operator - ( Vector u,Vector v ) { return Vector(u.x-v.x,u.y-v.y); }
    inline friend Vector operator * ( Vector u,double v ) { return Vector(u.x*v,u.y*v); }
    inline friend Vector operator / ( Vector u,double v ) { return Vector(u.x/v,u.y/v); }
    inline friend double operator * ( Vector u,Vector v ) { return u.x*v.x+u.y*v.y; }
    inline friend double operator ^ ( Vector u,Vector v ) { return u.x*v.y-u.y*v.x; }
    inline double calc_angle ( void ) { return atan2(y,x); }
}Q2[200000],p[200000];
struct Line
{
    Vector s,t;double A;
    Line(Vector _s=Vector(),Vector _t=Vector()){s=_s;t=_t;A=t.calc_angle();}
    inline friend bool operator < ( Line u,Line v )
    {
        if ( dcmp(u.A-v.A)!=0 ) return dcmp(u.A-v.A)==-1;
        return dcmp(u.t^(v.t+v.s-u.s))==-1;
    }
}L[200000],Q[200000];
inline bool PX ( Line u,Line v ) { return dcmp(u.t^v.t)==0; }
inline Vector JD ( Line u,Line v ) { return u.s+u.t*((v.t^(u.s-v.s))/(u.t^v.t)); }
inline bool OnRight ( Line u,Vector v ) { return dcmp(u.t^(v-u.s))<0; }
int n,cnt,m;long long X[200000],Y[200000];
signed main()
{
    scanf("%d",&n);double sum=0,ans=0;
    for ( int i=1;i<=n;i++ ) scanf("%lld%lld",&X[i],&Y[i]),p[i]=Vector(X[i],Y[i]);
    X[n+1]=X[1];Y[n+1]=Y[1];p[n+1]=p[1];
    for ( int i=1;i<=n;i++ ) L[++cnt]=Line(p[i],p[i+1]-p[i]),sum+=p[i]^p[i+1];
    for ( int i=2;i<=n;i++ )
    {
        double a=X[2]-X[1]+X[i]-X[i+1];
        double b=Y[2]-Y[1]+Y[i]-Y[i+1];
        double c=X[2]*Y[1]+X[i]*Y[i+1]-X[1]*Y[2]-X[i+1]*Y[i];
        if ( dcmp(a)!=0 ) L[++cnt]=Line(Vector(0,c/a),Vector(-a,-b));
        else if ( dcmp(b)!=0 ) L[++cnt]=Line(Vector(-c/b,0),Vector(-a,-b));
    }
    std::sort(L+1,L+cnt+1);
    int l=0,r=0;Q[0]=L[1];
    for ( int i=2;i<=cnt;i++ ) if ( dcmp(L[i].A-L[i-1].A) )
    {
        if ( l<r and ( PX(Q[l],Q[l+1]) or PX(Q[r],Q[r-1]) ) ) { ans=0;goto E; }
        while ( l<r and OnRight(L[i],Q2[r-1]) ) r--;
        while ( l<r and OnRight(L[i],Q2[l]) ) l++;
        Q[++r]=L[i];
        if ( l<r ) Q2[r-1]=JD(Q[r],Q[r-1]);
    }
    while ( l<r and OnRight(Q[l],Q2[r-1]) ) r--;
    while ( l<r and OnRight(Q[r],Q2[l]) ) l++;
    if ( r-l<=1 ) { ans=0;goto E; }
    Q2[r]=JD(Q[l],Q[r]);
    for ( int i=l;i<=r;i++ ) p[++m]=Q2[i];
    ans=p[m]^p[1];
    for ( int i=2;i<=m;i++ ) ans+=p[i-1]^p[i];
    E:printf("%.4lf
",ans/sum);
    return 0;
}

PS：由于答案在多边形内，一定不是无限大区域，可以不用加边界。

---

---

 Day5 20.01.31

本日任务：

算法一：凸包（复习）

算法二：旋转卡壳

算法三：平面图转对偶图

凸包

凸包属于计算几何的基本内容，在此不进行讲解，只选择几道题目来研究。

例题1：「DTOJ4212」旅行规划

【题目链接】http://59.61.75.5:8018/problem/4212

【题解】

此类不太好用线段树之类的数据结构维护的题一般都可以考虑分块。

分块后发现题意即为维护一个凸包。直接维护即可。

#include<bits/stdc++.h>
inline int read(void)
{
    int x = 0; char ch; bool f = true;
    while ( !isdigit(ch = getchar()) )if ( ch == '-' )f = false;
    for ( x = ch ^ 48; isdigit(ch = getchar());)x = ( x << 1 ) + ( x << 3 ) + ( ch ^ 48 );
    return f ? x : -x;
}
const long long inf = 10000000000000000LL;
const int maxn = 100000 + 10;
const int maxs = 500 + 10;
int n, s, cnt, bl[maxn], L[maxs], R[maxs], st[maxs], p[maxs][maxs], tot[maxs];
long long a[maxn], fst[maxs], d[maxs], add[maxs];
inline double calc(int x, int y) { return 1.0 *( a[y] - a[x] ) / ( y - x ); }
inline void build(int x)
{
    int tp = 1; st[tp] = L[x];
    for ( int i = L[x] + 1; i <= R[x]; i++ )
    {
        while ( tp >= 2 and calc(st[tp - 1], st[tp]) < calc(st[tp - 1], i) ) tp--;
        st[++tp] = i;
    }
    st[0] = 0; st[tp + 1] = n + 1; tot[x] = tp;
    for ( int i = 0; i <= tp + 1; i++ ) p[x][i] = st[i];
}
inline long long solve(int x)
{
    if ( x == 0 || x == n + 1 ) return -inf;
    return a[x] + fst[bl[x]] + d[bl[x]] * ( x - L[bl[x]] ) + add[bl[x]];
}
inline long long query(int x)
{
    int l = 1, r = tot[x];
    while ( l <= r )
    {
        int mid = ( l + r ) >> 1;
        long long t1 = solve(p[x][mid - 1]), t2 = solve(p[x][mid]), t3 = solve(p[x][mid + 1]);
        if ( t1 < t2 and t2 < t3 ) l = mid + 1;
        else if ( t1 > t2 and t2 > t3 ) r = mid - 1;
        else return t2;
    }
    return -1;
}
signed main()
{
    n = read();
    for ( int i = 1; i <= n; i++ ) a[i] = a[i - 1] + read();
    a[0] = a[n + 1] = -inf;
    s = (int) floor(sqrt(n) + 0.5); cnt = ( n - 1 ) / s + 1;
    for ( int i = 1; i <= n; i++ ) bl[i] = ( i - 1 ) / s + 1;
    for ( int i = 1; i <= cnt; i++ ) L[i] = ( i - 1 ) * s + 1, R[i] = std::min(n, i * s);
    for ( int i = 1; i <= cnt; i++ ) build(i);
    for ( int m = read(); m; m-- )
    {
        int op = read(), x = read(), y = read(), l = bl[x], r = bl[y];
        if ( op )
        {
            long long ans = -inf;
            if ( l == r )
            {
                for ( int i = x; i <= y; i++ ) ans = std::max(ans, solve(i));
            }
            else
            {
                for ( int i = l + 1; i <= r - 1; i++ ) ans = std::max(ans, query(i));
                for ( int i = x; i <= R[l]; i++ ) ans = std::max(ans, solve(i));
                for ( int i = L[r]; i <= y; i++ ) ans = std::max(ans, solve(i));
            }
            printf("%lld
", ans);
        }
        else
        {
            long long k = read();
            if ( l == r )
            {
                long long res = fst[l];
                for ( int i = L[l]; i <= R[l]; i++ ) a[i] += res, res += d[l], a[i] += add[l];
                fst[l] = d[l] = add[l] = 0;
                res = k;
                for ( int i = x; i <= y; i++ ) a[i] += res, res += k;
                build(l);
                res = k * ( y - x + 1 );
                for ( int i = y + 1; i <= R[l]; i++ ) a[i] += res;
                build(r);
                for ( int i = l + 1; i <= cnt; i++ ) add[i] += res;
            }
            else
            {
                long long res = k * ( L[l + 1] - x + 1 );
                for ( int i = l + 1; i <= r - 1; i++ ) fst[i] += res, d[i] += k, res += s * k;
                res = fst[l];
                for ( int i = L[l]; i <= R[l]; i++ ) a[i] += res, res += d[l], a[i] += add[l];
                fst[l] = d[l] = add[l] = 0;
                res = k;
                for ( int i = x; i <= std::min(y, R[l]); i++ ) a[i] += res, res += k;
                build(l);
                res = fst[r];
                for ( int i = L[r]; i <= R[r]; i++ ) a[i] += res, res += d[r], a[i] += add[r];
                fst[r] = d[r] = add[r] = 0;
                res = k * ( L[r] - x + 1 );
                for ( int i = L[r]; i <= y; i++ ) a[i] += res, res += k;
                res = k * ( y - x + 1 );
                for ( int i = y + 1; i <= R[r]; i++ ) a[i] += res;
                build(r);
                for ( int i = r + 1; i <= cnt; i++ ) add[i] += res;
            }
        }
    }
    return 0;
}

（代码是曾经有段时间用 VS 写的，格式化较严重）

例题2：「DTOJ4515」Generator 3

【题目链接】http://59.61.75.5:8018/problem/4515

【题解】

看到这么大的数据就知道肯定跟造数据方式有关。。。

考虑观察数据的生成。如果有做过类似题，不难发现会存在循环节。循环的方式是 $ ho$ 型的。$x,y$ 均按此方式循环，循环节小于 $10^5$。

因此 $ ho$ 前面一部分暴力处理，考虑解决两个环的问题。

思考凸包的做法。维护上凸壳、下凸壳。显然 $x$ 相同时，只有 $y$ 的最大最小值有用。在环上倍增维护即可。

#include<bits/stdc++.h>
inline int read ( void )
{
    int x=0;char ch;bool f=true;
    while ( !isdigit(ch=getchar()) ) if ( ch=='-' ) f=false;
    for ( x=ch^48;isdigit(ch=getchar()); ) x=(x<<1)+(x<<3)+(ch^48);
    return f ? x : -x ;
}
const int maxn=200000+10;
int X0,Y0,ax,ay,bx,by,px,py,vis[maxn],sx,sy,tx,ty,st,tp,s[maxn<<2],min[maxn<<1][35],max[maxn<<1][35],tot;
std::vector<int> C[maxn],Ax,Ay;
std::pair<int,int> p[maxn<<2];
long long n,ans;
inline long long cp ( std::pair<int,int> p1,std::pair<int,int> p2,std::pair<int,int> p3 )
{
    return 1LL*(p3.second-p1.second)*(p2.first-p1.first)-1LL*(p2.second-p1.second)*(p3.first-p1.first);
}
inline int getx ( int x )
{
    if ( x<sx ) return Ax[x];
    else return Ax[(x-sx)%tx+sx];
}
inline int gety ( int y )
{
    if ( y<sy ) return Ay[y];
    else return Ay[(y-sy)%ty+sy];
}
inline void init ( int id )
{
    int m=2*(int)C[id].size();
    for ( int i=0;i<m;i++ ) min[i][0]=max[i][0]=C[id][i%(int)C[id].size()];
    for ( int j=1;(1<<j)<=m;j++ ) for ( int i=0;i+(1<<j)<=m;i++ )
        min[i][j]=std::min(min[i][j-1],min[i+(1<<(j-1))][j-1]),
        max[i][j]=std::max(max[i][j-1],max[i+(1<<(j-1))][j-1]);
}
inline int rmqmin ( int l,int r )
{
    int k=log2(r-l+1);
    return std::min(min[l][k],min[r-(1<<k)+1][k]);
}
inline int rmqmax ( int l,int r )
{
    int k=log2(r-l+1);
    return std::max(max[l][k],max[r-(1<<k)+1][k]);
}
signed main()
{
    X0=read();Y0=read();ax=read();ay=read();bx=read();by=read();px=read();py=read();std::cin>>n;
    memset(vis,-1,sizeof(vis));
    vis[X0]=0;Ax.push_back(X0);
    for ( int i=1,res=X0;i<=(px<<1);i++ )
    {
        res=(1LL*res*ax+bx)%px;
        if ( ~vis[res] ) { sx=vis[res];tx=i-vis[res];break; }
        vis[res]=i;Ax.push_back(res);
    }
    memset(vis,-1,sizeof(vis));
    vis[Y0]=0;Ay.push_back(Y0);
    for ( int i=1,res=Y0;i<=(py<<1);i++ )
    {
        res=(1LL*res*ay+by)%py;
        if ( ~vis[res] ) { sy=vis[res];ty=i-vis[res];break; }
        vis[res]=i;Ay.push_back(res);
    }
    memset(vis,-1,sizeof(vis));
    int g=std::__gcd(tx,ty);st=std::max(sx,sy);
    for ( int i=0;i<g;i++ )
    {
        int tmp=st+i;
        for ( int j=0,k;j<tx/g;j++ ) k=getx(tmp),C[i].push_back(k),vis[k]=(int)C[i].size()-1,tmp=(tmp-st+ty)%tx+st;
    }
    for ( int i=0;i<g;i++ )
    {
        init(i);
        for ( int j=st+i;j<ty+st and j<n;j+=g )
            if ( n>=1LL*ty*(int)C[i].size()+j )
                p[++tot]=std::make_pair(rmqmin(0,(int)C[i].size()-1),gety(j)),
                p[++tot]=std::make_pair(rmqmax(0,(int)C[i].size()-1),gety(j));
            else
            {
                int k=vis[getx(j)],ed=(n-j-1)/ty;
                p[++tot]=std::make_pair(rmqmin(k,k+ed),gety(j)),
                p[++tot]=std::make_pair(rmqmax(k,k+ed),gety(j));
            }
    }
    for ( int i=0;i<st;i++ ) p[++tot]=std::make_pair(getx(i),gety(i));
    std::sort(p+1,p+tot+1);
    s[++tp]=1;
    for ( int i=2;i<=tot;i++ )
    {
        while ( tp>1 and cp(p[s[tp-1]],p[i],p[s[tp]])<=0 ) tp--;
        s[++tp]=i;
    }
    for ( int i=tp;i>=3;i-- ) ans+=cp(p[s[1]],p[s[i]],p[s[i-1]]);
    tp=0;
    s[++tp]=1;
    for ( int i=2;i<=tot;i++ )
    {
        while ( tp>1 and cp(p[s[tp-1]],p[i],p[s[tp]])>=0 ) tp--;
        s[++tp]=i;
    }
    for ( int i=3;i<=tp;i++ ) ans+=cp(p[s[1]],p[s[i-1]],p[s[i]]);
    return !printf("%lld
",ans);
}

例题3： 「FJWC2015」信用卡凸包

【题目链接】http://59.61.75.5:8018/problem/2197

【题解】

根据小学奥数计算此类曲边图形的周长方法，将圆弧部分与线段部分拆开考虑。

显然圆心角之和为 $2pi$，因此圆弧部分总长为 $2pi r$。

考虑将线段部分拼接在一起。向内平移至圆心处即可。

即每张卡拆成四个角的圆心，算凸包长度即可。

#include<bits/stdc++.h>
const double pi=acos(-1);
struct Graph
{
    struct Vector
    {
        double x,y;
        Vector(double _x=0,double _y=0){x=_x;y=_y;}
        inline friend Vector operator - ( const Vector &u,const Vector &v ) { return Vector(u.x-v.x,u.y-v.y); }
        inline friend double operator ^ ( const Vector &u,const Vector &v ) { return u.x*v.y-u.y*v.x; }
        inline double length ( void ) { return sqrt(x*x+y*y); }
    }st[1000000];int tp;
    std::vector<Vector> E;
    inline void add ( double x,double y ) { E.push_back(Vector(x,y)); }
    inline double solve ( void )
    {
        double ans=0;
        std::sort(E.begin(),E.end(),[&](const Vector &u,const Vector &v){return fabs(u.x-v.x)<1e-8 ? u.y<v.y : u.x<v.x;});
        for ( auto p:E )
        {
            while ( tp>1 and ((p-st[tp-1])^(st[tp]-st[tp-1]))<=0 ) tp--;
            st[++tp]=p;
        }
        for ( int i=2;i<=tp;i++ ) ans+=(st[i]-st[i-1]).length();
        tp=0;
        for ( auto p:E )
        {
            while ( tp>1 and ((p-st[tp-1])^(st[tp]-st[tp-1]))>=0 ) tp--;
            st[++tp]=p;
        }
        for ( int i=2;i<=tp;i++ ) ans+=(st[i]-st[i-1]).length();
        return ans;
    }
}G;
int n;double a,b,r,R,alpha;
signed main()
{
    scanf("%d%lf%lf%lf",&n,&a,&b,&r);a-=2*r;b-=2*r;
    R=sqrt(a*a+b*b)*0.5;alpha=atan(a/b);
    for ( int i=1;i<=n;i++ )
    {
        double x,y,theta;scanf("%lf%lf%lf",&x,&y,&theta);
        G.add(x+R*cos(theta+alpha),y+R*sin(theta+alpha));
        G.add(x+R*cos(theta-alpha),y+R*sin(theta-alpha));
        G.add(x-R*cos(theta+alpha),y-R*sin(theta+alpha));
        G.add(x-R*cos(theta-alpha),y-R*sin(theta-alpha));
    }
    return !printf("%.2lf
",G.solve()+2*pi*r);
}

旋转卡壳

例题：Beauty Contest

【题目链接】http://59.61.75.5:8018/problem/2165

【题解】

先给个网上的参考读音。。。

 

 再给个输入法的现状。$2^{3_{3^{3_3{3^{3_{3^{3_{3}}}}}}}}$。

思考有关距离的一般定义，都有个垂直。所以我们需要一个垂直。

 现在我们有个垂直。于是距离就很显然了：在凸包用一组平行线转一圈，平行线间距离最大值即为凸包的距离。

（到此结束上代码）

显然平行线不会超过 $O(n)$ 组（自己看图）。因此效率是正确的。

考虑如何实现转一圈的过程。用点来转显然非常麻烦，还需要做垂直。

我们考虑将平行线旋转一下：

现在看起来好多了！现在开始观察性质。

我们看看其它点到这条边的距离有什么性质：逆时针呈单峰函数！

于是就有这么一个做法：逆时针按顺序枚举边，点用推指针的方式找到（显然点是单调的），判定可以用面积法。效率 $O(n)$。

#include<bits/stdc++.h>
struct Vector
{
    int x,y;
    Vector(int _x=0,int _y=0){x=_x;y=_y;}
    inline friend Vector operator + ( Vector u,Vector v ) { return Vector(u.x+v.x,u.y+v.y); }
    inline friend Vector operator - ( Vector u,Vector v ) { return Vector(u.x-v.x,u.y-v.y); }
    inline friend long long operator ^ ( Vector u,Vector v ) { return 1LL*u.x*v.y-1LL*u.y*v.x; }
}G[50010],p[50010];
int n,tp,st[50010],m=-1;
std::set<int> s;
inline long long dis ( Vector u,Vector v ) { return 1LL*(u.x-v.x)*(u.x-v.x)+1LL*(u.y-v.y)*(u.y-v.y); }
signed main()
{
    scanf("%d",&n);
    for ( int i=1;i<=n;i++ ) scanf("%d%d",&p[i].x,&p[i].y);
    std::sort(p+1,p+n+1,[&](const Vector &p1,const Vector &p2){return p1.x==p2.x ? p1.y<p2.y : p1.x<p2.x;});
    tp=0;st[++tp]=1;
    for ( int i=2;i<=n;i++ )
    {
        while ( tp>1 and ((p[i]-p[st[tp-1]])^(p[st[tp]]-p[st[tp-1]]))<=0 ) tp--;
        st[++tp]=i;
    }
    for ( int i=1;i<=tp;i++ ) G[++m]=p[st[i]],s.insert(st[i]);
    tp=0;st[++tp]=1;
    for ( int i=2;i<=n;i++ )
    {
        while ( tp>1 and ((p[i]-p[st[tp-1]])^(p[st[tp]]-p[st[tp-1]]))>=0 ) tp--;
        st[++tp]=i;
    }
    for ( int i=tp;i;i-- ) if ( !s.count(st[i]) ) G[++m]=p[st[i]];
    std::reverse(G,G+m+1);G[++m]=G[0];
    long long ans=dis(G[0],G[1]);
    for ( int i=0,j=2;i<m;i++ )
    {
        while ( ((G[i+1]-G[i])^(G[j]-G[i]))<((G[i+1]-G[i])^(G[j+1]-G[i])) ) j=(j+1)%m;
        ans=std::max(ans,std::max(dis(G[i],G[j]),dis(G[i+1],G[j])));
    }
    return !printf("%lld
",ans);
}

PS：其实当坐标绝对值范围较小时，凸包上的点往往个数很少，因此以下算法活跃涌现：随机化、暴力枚举……（好像都能过）

平面图转对偶图

 例题1：「ZJOI2012」旅游

【题目链接】http://59.61.75.5:8018/problem/1232

【题解】

这题好像跟计算几何没多大关系。。

直接平面图转对偶图，显然结果是棵树，找树的直径即可。

#include<bits/stdc++.h>
const int maxn=200000+10;
std::vector<int> e[maxn];
int dis[maxn],n,ans,p;
std::map<std::pair<int,int>,int> map;
inline void Add ( int u,int v ) { e[u].push_back(v);e[v].push_back(u); }
inline void dfs ( int u,int fr )
{
    if ( (dis[u]=dis[fr]+1)>ans ) ans=dis[u],p=u;
    for ( int v:e[u] ) if ( v!=fr ) dfs(v,u);
}
signed main()
{
    scanf("%d",&n);n-=2;
    for ( int i=1;i<=n;i++ )
    {
        int p,q,r;scanf("%d%d%d",&p,&q,&r);
        if ( p>q ) std::swap(p,q);
        if ( q>r ) std::swap(q,r);
        if ( p>q ) std::swap(p,q);
        if ( p+1!=q )
        {
            if ( map.count(std::make_pair(p,q)) ) Add(map[std::make_pair(p,q)],i);
            map[std::make_pair(p,q)]=i;
        }
        if ( q+1!=r )
        {
            if ( map.count(std::make_pair(q,r)) ) Add(map[std::make_pair(q,r)],i);
            map[std::make_pair(q,r)]=i;
        }
        if ( p!=1 or r!=n+2 )
        {
            if ( map.count(std::make_pair(p,r)) ) Add(map[std::make_pair(p,r)],i);
            map[std::make_pair(p,r)]=i;
        }
    }
    dfs(1,0);dfs(p,0);
    return !printf("%d
",ans);
}

例题2：「HNOI2016」矿区

【题目链接】http://59.61.75.5:8018/problem/4547

【题解】

这题是非常标准的平面图转对偶图。

考虑转的过程。我们把每条边分离成两条：

原本是这样：

分离后是这样：

按照以往做法，每个多边形的边按逆时针排序。

这时我们注意到：

每条原来的边被拆成两条有向边，分别对应两个多边形。

把多边形当成新图的点，把边转换成连接两个多边形的新边即可。

考虑具体实现的过程：

首先需要找出每个多边形：对每个点的出边按极角序排序，从每个点出发沿最小转角绕一圈即可。

其次建图。

最后维护信息：

我们把无限大的部分当成一个点，任意找出对偶图的一个生成树。

不难证明，答案一定是在树上的一个连通块上。

先在树上对每条边定个方向，按给出点的顺序遍历每条边。若正向则答案加上子树答案和，若反向则减去。即可维护。

关于分数输出，可以上下同时 $ imes 4$，然后进行处理。

#include<bits/stdc++.h>
inline int read ( void )
{
    int x=0;char ch;bool f=true;
    while ( !isdigit(ch=getchar()) ) if ( ch=='-' ) f=false;
    for ( x=ch^48;isdigit(ch=getchar()); ) x=(x<<1)+(x<<3)+(ch^48);
    return f ? x : -x ;
}
const int maxn=200000+10,maxm=1200000+10;
const double eps=1e-10;
inline int dcmp ( double x ) { return fabs(x)<eps ? 0 : ( x<0 ? -1 : 1 ) ; }
int cnt=1,tot,root,nxt[maxm],bl[maxm],f[maxm],vis[maxm],tree[maxm],z[maxm];
long long s[maxm],S[maxm];
struct Point
{
    int x,y;
    Point(int _x=0,int _y=0){x=_x;y=_y;}
    inline friend Point operator + ( Point u,Point v ) { return Point(u.x+v.x,u.y+v.y); }
    inline friend Point operator - ( Point u,Point v ) { return Point(u.x-v.x,u.y-v.y); }
    inline friend long long operator ^ ( Point u,Point v ) { return 1LL*u.x*v.y-1LL*u.y*v.x; }
}p[maxn];
struct Edge { int id,u,v;double A; } e[maxm];
inline bool operator < ( Edge e1,Edge e2 ) { return dcmp(e1.A-e2.A)==0 ? e1.v<e2.v : dcmp(e1.A-e2.A)==-1 ; }
std::vector<Edge> h[maxn],E[maxm];
inline void addedge ( int u,int v )
{
    ++cnt;e[cnt]=(Edge){cnt,u,v,atan2(p[v].y-p[u].y,p[v].x-p[u].x)};h[u].push_back(e[cnt]);
    ++cnt;e[cnt]=(Edge){cnt,v,u,atan2(p[u].y-p[v].y,p[u].x-p[v].x)};h[v].push_back(e[cnt]);
}
inline void dfs ( int u,int fr )
{
    f[u]=fr;S[u]=1LL*s[u]*s[u];s[u]<<=1;vis[u]=1;
    for ( Edge e:E[u] ) if ( !vis[e.v] ) tree[e.id]=tree[e.id^1]=1,dfs(e.v,u),s[u]+=s[e.v],S[u]+=S[e.v];
}
inline long long gcd ( long long a,long long b ) { return b ? gcd(b,a%b) : a ; }
signed main()
{
    int n=read(),m=read(),q=read();
    for ( int i=1;i<=n;i++ ) p[i].x=read(),p[i].y=read();
    while ( m-- ) addedge(read(),read());
    for ( int i=1;i<=n;i++ ) std::sort(h[i].begin(),h[i].end());
    for ( int i=2;i<=cnt;i++ )
    {
        int v=e[i].v;
        auto it=std::lower_bound(h[v].begin(),h[v].end(),e[i^1]);
        if ( it==h[v].begin() ) it=h[v].end();
        nxt[i]=(*(--it)).id;
    }
    for ( int i=2;i<=cnt;i++ ) if ( !bl[i] )
    {
        bl[i]=bl[nxt[i]]=++tot;
        for ( int j=nxt[i];e[j].v!=e[i].u;j=nxt[j],bl[j]=tot ) s[tot]+=(p[e[j].u]-p[e[i].u])^(p[e[j].v]-p[e[i].u]);
        if ( s[tot]<=0 ) root=tot;
    }
    for ( int i=2;i<=cnt;i++ ) E[bl[i]].push_back((Edge){i,bl[i],bl[i^1],0});
    dfs(root,0);
    for ( long long P=0,Q=0;q--; )
    {
        int k=(read()+P)%n+1;
        for ( int i=1;i<=k;i++ ) z[i]=(read()+P)%n+1;
        z[k+1]=z[1];P=Q=0;
        for ( int i=1;i<=k;i++ )
        {
            int x=z[i],y=z[i+1];
            auto it=lower_bound(h[x].begin(),h[x].end(),(Edge){0,x,y,atan2(p[y].y-p[x].y,p[y].x-p[x].x)});
            int j=(*it).id;
            if ( !tree[j] ) continue;
            if ( f[bl[j]]==bl[j^1] ) P+=S[bl[j]],Q+=s[bl[j]];
            else P-=S[bl[j^1]],Q-=s[bl[j^1]];
        }
        long long g=gcd(P,Q);P/=g;Q/=g;
        printf("%lld %lld
",P,Q);
    }
    return 0;
}

---

---

Day6 20.02.01

本日任务：

算法一：单位圆覆盖

算法二：仿射变换与矩阵

单位圆覆盖

【题目链接】http://59.61.75.5:8018/problem/1698

【题解】

方法一：

显然一定存在一个最优解使得至少两个点在圆上。

$O(n^2)$ 枚举两个点，找到圆心，$O(n)$ 判断每个点是否在圆内即可。

效率 $O(n^3)$。

#include<bits/stdc++.h>
const int eps=1e-10;
struct Vector
{
    double x,y;
    Vector(double _x=0,double _y=0){x=_x;y=_y;}
    inline friend Vector operator + ( Vector u,Vector v ) { return Vector(u.x+v.x,u.y+v.y); }
    inline friend Vector operator - ( Vector u,Vector v ) { return Vector(u.x-v.x,u.y-v.y); }
    inline friend Vector operator * ( Vector u,double v ) { return Vector(u.x*v,u.y*v); }
    inline friend Vector operator / ( Vector u,double v ) { return Vector(u.x/v,u.y/v); }
    inline double Len ( void ) { return sqrt(x*x+y*y); }
    inline Vector unit ( void ) { return *this/Len(); }
}p[1000],C;
int n;
inline int calc ( void )
{
    int res=0;
    for ( int i=1;i<=n;i++ ) if ( (C-p[i]).Len()<=1 ) res++;
    return res;
}
signed main()
{
    while ( ~scanf("%d",&n) and n )
    {
        int ans=1;
        for ( int i=1;i<=n;i++ )
        {
            double x,y;
            scanf("%lf%lf",&x,&y);
            p[i]=Vector(x,y);
        }
        for ( int i=1;i<n;i++ ) for ( int j=i+1;j<=n;j++ ) if ( (p[i]-p[j]).Len()<=2 )
        {
            Vector mid=(p[i]+p[j])/2;
            Vector d=p[j]-p[i];
            Vector a=Vector(d.y,-d.x).unit();
            double dis=sqrt(1-(p[i]-mid).Len()*(p[i]-mid).Len());
            a=a*dis;
            C=mid+a;ans=std::max(ans,calc());
            C=mid-a;ans=std::max(ans,calc());
        }
        printf("%d
",ans);
    }
    return 0;
}

方法二：

考虑先确定其中一个点。

显然圆心应选在以确定点为圆心，单位长为半径的圆上。因此考虑每段弧的贡献。

把另一点转化为相交的两个端点的弧，则以弧上任一点为圆心的圆必定包含另一点。

因此每个点转化为一段区间，求最大区间覆盖数，排序后扫一遍即可。效率 $O(n^2 log n)$。

#include<bits/stdc++.h>
const double pi=acos(-1);
int n;
struct Vector
{
    double x,y;
    Vector(double _x=0,double _y=0){x=_x;y=_y;}
    inline friend Vector operator + ( Vector u,Vector v ) { return Vector(u.x+v.x,u.y+v.y); }
    inline friend Vector operator - ( Vector u,Vector v ) { return Vector(u.x-v.x,u.y-v.y); }
    inline double Len ( void ) { return sqrt(x*x+y*y); }
    inline double calc ( void ) { return atan2(y,x); }
} p[1000];
struct Angle
{
    double A;bool f;
    Angle(double _A=0,bool _f=false){A=_A;f=_f;}
    inline bool friend operator < (  const Angle &A1,const Angle &A2 ) { return A1.A<A2.A; }
} q[1000];
signed main()
{
    while ( ~scanf("%d",&n) and n )
    {
        int ans=1;
        for ( int i=1;i<=n;i++ )
        {
            double x,y;
            scanf("%lf%lf",&x,&y);
            p[i]=Vector(x,y);
        }
        for ( int i=1;i<=n;i++ )
        {
            int tot=0;
            for ( int j=1;j<=n;j++ ) if ( i!=j and (p[i]-p[j]).Len()<=2 )
            {
                double res=acos((p[i]-p[j]).Len()/2);
                double ang=(p[j]-p[i]).calc();
                q[++tot]=Angle(ang-res,true);
                q[++tot]=Angle(ang+res,false);
            }
            if ( tot<=ans ) continue;
            std::sort(q+1,q+tot+1);
            int cnt=1;
            for ( int j=1;j<=tot;j++ )
            {
                if ( q[j].f ) cnt++;
                else cnt--;
                ans=std::max(ans,cnt);
            }
        }
        printf("%d
",ans);
    }
}

仿射变换与矩阵

从本专题起已经进入计算几何的冷门内容，下述两个算法均参考刘汝佳著作《算法竞赛入门经典——训练指南》。

仿射变换的定义：https://baike.baidu.com/item/%E4%BB%BF%E5%B0%84%E5%8F%98%E6%8D%A2/4289056?fr=aladdin

简单来说，仿射变换就是线性变换与平移的组合。

线性变换的定义：https://baike.baidu.com/item/%E7%BA%BF%E6%80%A7%E5%8F%98%E6%8D%A2/5904192

简单来说，线性变换就是旋转、放缩等基本变换。

高中选修数学中有提及关于此部分的内容，提及了用矩阵来表示变换（其实本质是因为向量本身就是矩阵的一种形式）。

现在给出具体的做法：

我们用矩阵：$$egin{bmatrix}x\yend{bmatrix}$$ 表示点 $(x,y)$。

那么有（下列变换一律采用左乘的形式）：

变换一：旋转 $alpha$ 度（角度一律按照高中数学课本定义）对应的矩阵为：$$egin{bmatrix}cos{alpha}&-sin{alpha}\sin{alpha}&cos{alpha}end{bmatrix}$$

变换二：放缩至 $k$ 倍对应的矩阵为：$$egin{bmatrix}k&0\0&kend{bmatrix}$$

其余线性变换均可用上述变换复合得成，例如反射即为放缩至 $-1$ 倍。

但注意到：平移无法被表示。于是我们引入第三维参数 $1$，进行平移变换，有：

用矩阵：$$egin{bmatrix}x\y\1end{bmatrix}$$ 表示点 $(x,y)$。

变换一：旋转 $alpha$ 度对应的矩阵为：$$egin{bmatrix}cos{alpha}&-sin{alpha}&0\sin{alpha}&cos{alpha}&0\0&0&1end{bmatrix}$$

变换二：放缩至 $k$ 倍对应的矩阵为：$$egin{bmatrix}k&0&0\0&k&0\0&0&1end{bmatrix}$$

变换三：平移 $(dx,dy)$ 对应的矩阵为：$$egin{bmatrix}1&0&dx\0&1&dy\0&0&1end{bmatrix}$$

于是我们就可以用矩阵乘法的形式表示所有变换了。

同时在这种做法下，我们也能真正地区分点与向量的区别，即第三维参数为 $0$ 时为向量，为 $1$ 时为点。

这样同时也满足了点与向量的关系：

（1）向量+向量=向量

（2）点+向量=点

（3）点+点无意义

上述做法被称为点的齐次坐标，若有兴趣的读者可自行查阅其余资料了解。

例题：UVa12303

【题目链接】https://vjudge.net/problem/UVA-12303

【题解】

依题意，直接用仿射矩阵记录每次操作即可。

对于平面，应用三点式而非点法式储存。使用点法式直接变换会导致答案错误，因为法向量无法按此规则变换，且平面外一点变换后与平面的相对位置关系改变。

#include<bits/stdc++.h>
const double eps=1e-8;
const double pi=acos(-1.0);
inline int dcmp ( double x ) { return fabs(x)<eps ? 0 : ( x<0 ? -1 : 1 ) ; }
struct matrix
{
    int n,m;double a[10][10];
    matrix(){}
    inline double * operator [] ( const int x ) { return a[x]; }
    inline void init ( int _n,int _m,bool flag )
    {
        n=_n;m=_m;
        for ( int i=1;i<=n;i++ ) for ( int j=1;j<=m;j++ ) a[i][j]=0;
        if ( flag ) for ( int i=1;i<=n;i++ ) a[i][i]=1;
    }
    inline friend matrix operator * ( matrix A,matrix B )
    {
        matrix C;C.init(A.n,B.m,false);
        for ( int i=1;i<=A.n;i++ ) for ( int k=1;k<=B.m;k++ ) for ( int j=1;j<=A.m;j++ ) C[i][k]+=A[i][j]*B[j][k];
        return C; 
    }
};
struct Point
{
    double x,y,z;
    Point(double _x=0,double _y=0,double _z=0){x=_x;y=_y;z=_z;}
    inline friend Point operator + ( Point u,Point v ) { return Point(u.x+v.x,u.y+v.y,u.z+v.z); }
    inline friend Point operator - ( Point u,Point v ) { return Point(u.x-v.x,u.y-v.y,u.z-v.z); }
    inline friend Point operator * ( Point u,double v ) { return Point(u.x*v,u.y*v,u.z*v); }
    inline friend Point operator / ( Point u,double v ) { return Point(u.x/v,u.y/v,u.z/v); }
    inline friend double operator * ( Point u,Point v ) { return u.x*v.x+u.y*v.y+u.z*v.z; }
    inline friend Point operator ^ ( Point u,Point v ) { return Point(u.y*v.z-u.z*v.y,u.z*v.x-u.x*v.z,u.x*v.y-u.y*v.x); }
    inline friend bool operator == ( Point u,Point v ) { return dcmp(u.x-v.x)==0 and dcmp(u.y-v.y)==0 and dcmp(u.z-v.z)==0; }
    inline double Len ( void ) { return sqrt(x*x+y*y+z*z); }
    inline Point unit ( void ) { return *this/Len(); }
} p[100000];
inline Point modify ( Point u,matrix A )
{
    matrix B;B.init(4,1,false);
    B[1][1]=u.x;B[2][1]=u.y;B[3][1]=u.z;B[4][1]=1;
    A=A*B;
    return Point(A[1][1],A[2][1],A[3][1]);
}
struct Plane { double a,b,c,d; } T[100000];
signed main()
{
    int n,m,Q;scanf("%d%d%d",&n,&m,&Q);
    for ( int i=1;i<=n;i++ )
    {
        double a,b,c;
        scanf("%lf%lf%lf",&a,&b,&c);
        p[i]=Point(a,b,c);
    }
    for ( int i=1;i<=m;i++ ) scanf("%lf%lf%lf%lf",&T[i].a,&T[i].b,&T[i].c,&T[i].d);
    matrix Ans;Ans.init(4,4,true);
    while ( Q-- )
    {
        char op[20]={''};scanf(" %s",op);
        matrix trans;trans.init(4,4,true);
        if ( op[0]=='T' ) scanf("%lf%lf%lf",&trans[1][4],&trans[2][4],&trans[3][4]);
        else if ( op[0]=='R' )
        {
            double a,b,c,d;scanf("%lf%lf%lf%lf",&a,&b,&c,&d);
            Point P=Point(a,b,c).unit();
            a=P.x;b=P.y;c=P.z;d=d/180*pi;
            double Cos=cos(d),Sin=sin(d);
            trans[1][1]=(1-Cos)*a*a+Cos;trans[1][2]=(1-Cos)*a*b-Sin*c;trans[1][3]=(1-Cos)*a*c+Sin*b;
            trans[2][1]=(1-Cos)*b*a+Sin*c;trans[2][2]=(1-Cos)*b*b+Cos;trans[2][3]=(1-Cos)*b*c-Sin*a;
            trans[3][1]=(1-Cos)*c*a-Sin*b;trans[3][2]=(1-Cos)*c*b+Sin*a;trans[3][3]=(1-Cos)*c*c+Cos;
        }
        else if ( op[0]=='S' ) scanf("%lf%lf%lf",&trans[1][1],&trans[2][2],&trans[3][3]);
        Ans=trans*Ans;
    }
    for ( int i=1;i<=n;i++ ) p[i]=modify(p[i],Ans),printf("%.2lf %.2lf %.2lf
",p[i].x,p[i].y,p[i].z);
    for ( int i=1;i<=m;i++ )
    {
        double a=T[i].a,b=T[i].b,c=T[i].c,d=-T[i].d;
        double max=std::max(fabs(a),std::max(fabs(b),fabs(c)));
        Point P[4];
        if ( fabs(a)==max ) P[1]=Point(d/a,0,0),P[2]=Point((d-b)/a,1,0),P[3]=Point((d-c)/a,0,1);
        else if ( fabs(b)==max ) P[1]=Point(0,d/b,0),P[2]=Point(1,(d-a)/b,0),P[3]=Point(0,(d-c)/b,1);
        else P[1]=Point(0,0,d/c),P[2]=Point(1,0,(d-a)/c),P[3]=Point(0,1,(d-b)/c);
        for ( int i=1;i<=3;i++ ) P[i]=modify(P[i],Ans);
        Point u=((P[3]-P[1])^(P[2]-P[1])).unit();
        printf("%.2lf %.2lf %.2lf %.2lf
",u.x,u.y,u.z,-(u*P[1]));
    }
    return 0;
}

---

---

 Day7 20.02.02

 本日任务：算法：运动规划（一）

运动规划（一）

*在本篇 $blog$ 中，模拟退火与运动规划的两部分内容已经略微超越平常的算法竞赛，与实际应用最为紧密。

运动规划是来自于机器人领域的一个问题，一般综合考查图论与计算几何，代码复杂，容易出错。

实际上，最短路模型是运动规划里最为简单、最为基础的内容。下面的例题将会带来最短路与计算几何的结合内容。

例题1：UVa11921

【题目链接】https://vjudge.net/problem/UVA-11921

【题解】

题意：给定平面上两点与 $n$ 个障碍（障碍为圆形），求两点间的最短路径。

这是个最短路问题，直接上 $SPFA$ 就完了。

这是一道典型的运动规划问题。考虑将其转化为一般最短路模型求解。

显然随机游走的答案是不优的。考虑实际应当怎么走。

显然最优路线一定是沿直线走或沿障碍物走。

因此考虑在障碍物、起点 $S$、终点 $T$ 间连边。有以下四类：

① $S-T$ 连边

② $S/T$ 向障碍物连切线

③ 障碍物与障碍物连公切线

④ 障碍物上不同点连边

分别处理即可。考虑如何处理。

第①类，只需判断线段是否与障碍物有交点。

第②类，作切线，转化为连线段，同①。

第③类，作切线，同②。

第④类，考虑障碍物上哪些点会被走到。第②③类的点记录一下，按极角序排个序，相邻两个连边就可以了。

然后跑你喜欢的最短路算法 $SPFA$ 或 $dijkstra$ 即可。

顺路提醒一下，代码强烈建议封装。期望长度： $6k$。

#include<bits/stdc++.h>
const double eps=1e-8;
const double pi=acos(-1);
inline int dcmp ( double x ) { return fabs(x)<eps ? 0 : ( x<0 ? -1 : 1 ) ; }
struct Vector
{
    double x,y;
    Vector(double _x=0,double _y=0){x=_x;y=_y;}
    inline friend Vector operator + ( Vector u,Vector v ) { return Vector(u.x+v.x,u.y+v.y); }
    inline friend Vector operator - ( Vector u,Vector v ) { return Vector(u.x-v.x,u.y-v.y); }
    inline friend Vector operator * ( Vector u,double v ) { return Vector(u.x*v,u.y*v); }
    inline friend Vector operator / ( Vector u,double v ) { return Vector(u.x/v,u.y/v); }
    inline friend double operator * ( Vector u,Vector v ) { return u.x*v.x+u.y*v.y; }
    inline friend double operator ^ ( Vector u,Vector v ) { return u.x*v.y-u.y*v.x; }
    inline friend bool operator == ( Vector u,Vector v ) { return dcmp(u.x-v.x)==0 and dcmp(u.y-v.y)==0; }
    inline friend bool operator < ( Vector u,Vector v )
    {
        if ( dcmp(u.x-v.x)!=0 ) return dcmp(u.x-v.x)==-1;
        else return dcmp(u.y-v.y)==-1;
    }
    inline Vector rotate ( double rad ) { return Vector(x*cos(rad)-y*sin(rad),x*sin(rad)+y*cos(rad)); }
    inline double Len ( void ) { return sqrt(x*x+y*y); }
    inline double angle ( void ) { return atan2(y,x); }
    inline Vector unit ( void ) { return *this/Len(); }
}S,T;
struct Segment
{
    Vector s,t;
    Segment(Vector _s=Vector(),Vector _t=Vector()){s=_s;t=_t;}
    inline double Len ( void ) { return (t-s).Len(); }
    inline double dist_line ( Vector u ) { return fabs((u-s)^(u-t))/Len(); }
    inline double dist_segment ( Vector u )
    {
        if ( s==t ) return (u-s).Len();
        Vector A=u-s,B=u-t;
        if ( dcmp((t-s)*A)==-1 ) return A.Len();
        else if ( dcmp((t-s)*B)==1 ) return B.Len();
        else return fabs(A^B)/Len();
    }
};
struct Line
{
    Vector p,v;
    Line(Vector _p=Vector(),Vector _v=Vector()){p=_p;v=_v;}
    inline Vector OnLine ( Vector u ) { return p+v*((v*(u-p))/(v*v)); } 
};
struct Circle
{
    Vector O;
    double r;
    Circle(Vector _O=Vector(),double _r=0){O=_O;r=_r;}
    inline Vector point ( double rad ) { return O+Vector(cos(rad),sin(rad))*r; }
    inline int On_Cricle ( Vector p ) { return dcmp((p-O).Len()-r)==0; }
    inline int Tangent_Point ( Vector p,Vector *res )
    {
        double d=(O-p).Len();
        if ( dcmp(d-r)==-1 ) return 0;
        else if ( dcmp(d-r)==0 )
        {
            res[0]=(O-p).rotate(pi/2);
            return 1;
        }
        else
        {
            double angle=asin(r/d);
            res[0]=(O-p).rotate(-angle);
            res[1]=(O-p).rotate(angle);
            return 2;
        }
    }
}C[100000];
inline int Tangent_Circle ( Circle u,Circle v,Segment *res )
{
    if ( u.r<v.r ) std::swap(u,v);
    int cnt=0;
    double d=(u.O-v.O).Len(),del=fabs(u.r-v.r),sum=u.r+v.r;
    if ( del>d ) return 0; 
    double rad=(v.O-u.O).angle(),theta=acos(del/d);
    res[cnt++]=Segment(u.point(rad+theta),v.point(rad+theta));
    res[cnt++]=Segment(u.point(rad-theta),v.point(rad-theta));
    if ( dcmp(d-sum)==0 ) res[cnt++]=Segment(u.point(rad),v.point(rad+pi));
    else
    {
        double angle=acos(sum/d);
        res[cnt++]=Segment(u.point(rad+angle),v.point(rad+angle+pi));
        res[cnt++]=Segment(u.point(rad-angle),v.point(rad-angle+pi));
    }
    return cnt;
}
struct edge { int v,nxt;double w; } e[100000];
int e_cnt,h[100000],n;
inline void addedge ( int u,int v,double w )
{
    e[++e_cnt].nxt=h[u];e[h[u]=e_cnt].v=v;e[e_cnt].w=w;
    e[++e_cnt].nxt=h[v];e[h[v]=e_cnt].v=u;e[e_cnt].w=w;
}
inline double Get_Angle ( double rad )
{
    while ( rad<0 ) rad+=2*pi;
    while ( rad>2*pi ) rad-=2*pi;
    return std::min(rad,2*pi-rad);
}
std::vector<Vector> Points;
std::map<Vector,int> map;
inline int Get_Id ( Vector p )
{
    if ( !map.count(p) ) Points.push_back(p),map[p]=(int)Points.size()-1;
    return map[p];
}
const double inf=1e12;
double dis[100000];
bool vis[100000];
inline int On_Circles ( Vector p )
{
    int res=0;
    for ( int i=1;i<=n;i++ ) if ( C[i].On_Cricle(p) ) res++;
    return res;
}
inline void Add_Segment ( Segment segment,int cnt )
{
    for ( int i=1;i<=n;i++ )
    {
        double d=segment.dist_segment(C[i].O);
        if ( dcmp(d-C[i].r)==-1 or ( dcmp(d-C[i].r)==0 and cnt--==0 ) ) return;
    }
    int u=Get_Id(segment.s),v=Get_Id(segment.t);
    if ( u>=0 and v>=0 and u!=v ) addedge(u,v,segment.Len());
}
inline void Add_Tangent ( Vector p )
{
    for ( int i=1;i<=n;i++ )
    {
        Vector res[2];int cnt=C[i].Tangent_Point(p,res);
        for ( int j=0;j<cnt;j++ ) Add_Segment(Segment(p,Line(p,res[j]).OnLine(C[i].O)),1);
    }
}
struct Double
{
    double val;
    Double(double _val=0){val=_val;}
    inline friend bool operator == ( Double x,Double y ) { return dcmp(x.val-y.val)==0; }
    inline friend bool operator < ( Double x,Double y ) { return dcmp(x.val-y.val)==-1; }
};
inline double solve ( void )
{
    int s=Get_Id(S),t=Get_Id(T);
    Add_Segment(Segment(S,T),0);
    Add_Tangent(S);Add_Tangent(T);
    for ( int i=1;i<n;i++ ) for ( int j=i+1;j<=n;j++ )
    {
        Segment res[4];int cnt=Tangent_Circle(C[i],C[j],res);
        for ( int k=0;k<cnt;k++ ) Add_Segment(res[k],2);
    }
    for ( int i=1;i<=n;i++ )
    {
        std::vector<double> Angle;
        for ( Vector p:Points ) if ( C[i].On_Cricle(p) ) Angle.push_back((p-C[i].O).angle());
        std::sort(Angle.begin(),Angle.end());
        Angle.resize(unique(Angle.begin(),Angle.end())-Angle.begin());
        if ( (int)Angle.size()<2 ) continue;
        Angle.push_back(Angle[0]);
        for ( int j=1;j<(int)Angle.size();j++ )
        {
            double dist=C[i].r*Get_Angle(Angle[j]-Angle[j-1]);
            int u=Get_Id(C[i].point(Angle[j])),v=Get_Id(C[i].point(Angle[j-1]));
            if ( u>=0 and v>=0 and u!=v ) addedge(u,v,dist);
        }
    }
    for ( int i=0;i<(int)Points.size();i++ ) dis[i]=inf,vis[i]=false;
    std::queue<int> q;q.push(s);dis[s]=0;vis[s]=true;
    while ( !q.empty() )
    {
        int u=q.front();q.pop();vis[u]=false;
        if ( On_Circles(Points[u])>=2 ) continue;
        for ( int i=h[u];i;i=e[i].nxt ) if ( dis[e[i].v]>dis[u]+e[i].w )
        {
            dis[e[i].v]=dis[u]+e[i].w;
            if ( !vis[e[i].v] ) q.push(e[i].v),vis[e[i].v]=true;
        }
    }
    return dis[t];
}
signed main()
{
    int Q;scanf("%d",&Q);
    for ( int Case=1;Case<=Q;Case++ )
    {
        scanf("%lf%lf%lf%lf%d",&S.x,&S.y,&T.x,&T.y,&n);
        for ( int i=1;i<=n;i++ ) scanf("%lf%lf%lf",&C[i].O.x,&C[i].O.y,&C[i].r);
        printf("Case %d: %.8lf
",Case,solve());
        for ( int i=0;i<(int)Points.size();i++ ) h[i]=0;
        Points.clear();map.clear();e_cnt=0;
    }
    return 0;
}

---

---

 Day8 20.02.03

本日任务：算法：运动规划（二）

运动规划（一）

例题2：「JOISC 2018 Day1」栅栏

【题目链接】http://59.61.75.5:8018/problem/4198

【题解】

还是按照运动规划的一般思路分析，考虑转化为最短路模型。

考虑连边。显然有三种连边：

①栅栏与栅栏连边

②栅栏与牧场连边

③牧场与牧场连边

第①类，需要分类讨论栅栏间的最短距离：垂线段、端点连线。

第②类，放上官方题解的一张图：

这就是九大类情况。

此外还有三种形式的连边：

 

 

分类讨论即可。

第三类即直接连边。

考虑答案的计算。

显然要求一个多边形，能够包裹住牧场。

也就是说，牧场应在多边形内部。

判定多边形内部常用射线法，因此任意作一条射线 $OM$。

连边时判定与 $OM$ 的相交情况，最后跑最短路时令交点数量为奇数即可。

具体细节参考代码。

#include<bits/stdc++.h>
const double eps=1e-8;
const double inf=1e9;
const double pi=acos(-1);
inline int dcmp ( double x ) { return fabs(x)<eps ? 0 : ( x<0 ? -1 : 1 ) ; }
int n,S;double dis[1000][1000][2],Dis[1000][1000][2];
struct Vector
{
    double x,y;
    Vector(double _x=0,double _y=0){x=_x;y=_y;}
    inline friend Vector operator + ( Vector u,Vector v ) { return Vector(u.x+v.x,u.y+v.y); }
    inline friend Vector operator - ( Vector u,Vector v ) { return Vector(u.x-v.x,u.y-v.y); }
    inline friend Vector operator * ( Vector u,double v ) { return Vector(u.x*v,u.y*v); }
    inline friend Vector operator / ( Vector u,double v ) { return Vector(u.x/v,u.y/v); }
    inline friend double operator * ( Vector u,Vector v ) { return u.x*v.x+u.y*v.y; }
    inline friend double operator ^ ( Vector u,Vector v ) { return u.x*v.y-u.y*v.x; }
    inline friend bool operator == ( Vector u,Vector v ) { return dcmp(u.x-v.x)==0 and dcmp(u.y-v.y)==0; }
    inline friend bool operator < ( Vector u,Vector v )
    {
        if ( dcmp(u.x-v.x)!=0 ) return dcmp(u.x-v.x)==-1;
        else return dcmp(u.y-v.y)==-1;
    }
    inline Vector rotate ( double rad ) { return Vector(x*cos(rad)-y*sin(rad),x*sin(rad)+y*cos(rad)); }
    inline double length ( void ) { return sqrt(x*x+y*y); }
    inline double angle ( void ) { return atan2(y,x); }
    inline Vector unit ( void ) { return *this/length(); }
};
struct Segment
{
    Vector s,t;
    Segment(Vector _s=Vector(),Vector _t=Vector()){s=_s;t=_t;}
    inline double length ( void ) { return (t-s).length(); }
    inline Vector dist_segment ( Vector u )
    {
        double val=((t-s)*(u-s))/length()/length();
        if ( val<0 ) return s;
        if ( val>1 ) return t;
        else return s+(t-s)*val;
    }
}Line_M,L[200];
inline bool ifsect ( Segment u,Segment v )
{
    return dcmp(((v.s-u.s)^(u.t-u.s))*((v.t-u.s)^(u.t-u.s)))<0 and dcmp(((u.s-v.s)^(v.t-v.s))*((u.t-v.s)^(v.t-v.s)))<=0;
}
inline bool ifsectbound ( Segment p )
{
    if ( ifsect(p,Segment(Vector(-S,-S),Vector(-S,S))) ) return true;
    if ( ifsect(p,Segment(Vector(-S,-S),Vector(S,-S))) ) return true;
    if ( ifsect(p,Segment(Vector(-S,S),Vector(S,S))) ) return true;
    if ( ifsect(p,Segment(Vector(S,-S),Vector(S,S))) ) return true;
    if ( ifsect(p,Segment(Vector(-S,-S),Vector(S,S))) ) return true;
    if ( ifsect(p,Segment(Vector(S,-S),Vector(-S,S))) ) return true;
    return false;
}
inline void Solve_Fences ( Segment u,Segment v,int x ,int y )
{
    Vector pos;bool w=false;
    dis[x][y][0]=dis[x][y][1]=inf;
    
    pos=u.dist_segment(v.s);
    w=ifsect(Line_M,Segment(v.s,pos))^ifsect(Line_M,Segment(u.s,pos));
    if ( !ifsectbound(Segment(pos,v.s)) ) dis[x][y][w]=std::min(dis[x][y][w],(v.s-pos).length());
    
    pos=u.dist_segment(v.t);
    w=ifsect(Line_M ,v)^ifsect(Line_M,Segment(pos,v.t))^ifsect(Line_M,Segment(u.s,pos));
    if ( !ifsectbound(Segment(pos,v.t)) ) dis[x][y][w]=std::min(dis[x][y][w],(v.t-pos).length());

    pos=v.dist_segment(u.s);
    w=ifsect(Line_M,Segment(u.s,pos))^ifsect(Line_M,Segment(v.s,pos));
    if ( !ifsectbound(Segment(pos,u.s)) ) dis[x][y][w]=std::min(dis[x][y][w],(u.s-pos).length());

    pos=v.dist_segment(u.t);
    w=ifsect(Line_M,u)^ifsect(Line_M,Segment(u.t,pos))^ifsect(Line_M,Segment(v.s,pos));
    if ( !ifsectbound(Segment(pos,u.t)) ) dis[x][y][w]=std::min(dis[x][y][w],(u.t-pos).length());
}

int main()
{
    srand(time(NULL));scanf("%d%d",&n,&S);
    double Angle=2*pi*rand()/RAND_MAX;
    Line_M=Vector(500*cos(Angle),500*sin(Angle));
    for ( int i=1,X1,Y1,X2,Y2;i<=n;i++ ) scanf("%d%d%d%d",&X1,&Y1,&X2,&Y2),L[i]=Segment(Vector(X1,Y1),Vector(X2,Y2));
    L[++n]=Segment(Vector(S,S),Vector(S+eps,S));
    L[++n]=Segment(Vector(S,-S),Vector(S+eps,-S));
    L[++n]=Segment(Vector(-S,S),Vector(-S-eps,S));
    L[++n]=Segment(Vector(-S,-S),Vector(-S-eps,-S));
    for ( int i=1;i<=n;i++ )
    {
        Dis[i][i][1]=dis[i][i][1]=inf;
        for ( int j=i+1;j<=n;j++ )
        {
            Solve_Fences(L[i],L[j],i,j);
            for ( int k=0;k<2;k++ ) Dis[i][j][k]=Dis[j][i][k]=dis[j][i][k]=dis[i][j][k];
        }
    }
    double Ans=inf;
    for ( int i=1;i<=n;i++ )
    {
        for ( int j=1;j<i;j++ ) for ( int k=j+1;k<i;k++ ) for ( int p=0;p<2;p++ ) for ( int q=0;q<2;q++ )
            Ans=std::min(Ans,Dis[j][k][p^q^1]+dis[i][j][p]+dis[i][k][q]);
        for ( int j=1;j<=n;j++ ) for ( int k=1;k<=n;k++ ) for ( int p=0;p<2;p++ ) for ( int q=0;q<2;q++ )
            Dis[j][k][p^q]=std::min(Dis[j][k][p^q],Dis[i][j][p]+Dis[i][k][q]);
    }
    return !printf("%.10lf
",Ans);
}