计算几何板子

半平面交

没什么好说的。。计算几何是要多做些题才能对向量熟练运用

#include<bits/stdc++.h>
using namespace std;
const int N=1010;
const double eps=1e-9;
int n,m,cnt,q[N],l,r;double ans;
struct Poi{
    double x,y;
    Poi operator -(const Poi &b){return (Poi){x-b.x,y-b.y};}
    Poi operator +(const Poi &b){return (Poi){x+b.x,y+b.y};}
    double operator ^(const Poi &b){return x*b.y-y*b.x;}
    Poi operator *(const double &b){return (Poi){x*b,y*b};}
}b[N];
struct Line{
    Poi a,b,c;double d;//向量ab左侧的半平面
    void get(){c=b-a;d=atan2(c.y,c.x);}
    bool friend operator <(Line a,Line b){
        return abs(a.d-b.d)>eps?a.d<b.d:(a.c^(b.b-a.a))>eps;//先按极角排序，然后再按相对位置排序
    }
}a[N],c[N];
Poi inter(Line a,Line b){return b.a+b.c*(((b.a-a.a)^a.c)/(a.c^b.c));}//叉积求面积，相似求交点。画图可得。
bool check(Line a,Line b,Line c){return (c.c^(inter(a,b)-c.a))<-eps;}//判a和b的交点是否在c的半平面内
int read(){
    int x=0,c,f=1;
    while(!isdigit(c=getchar()))c=='-'?f=-1:0;
    while(isdigit(c))x=x*10+c-48,c=getchar();
    return x*f;
}
int main(){
    n=read();
    for(int i=1,x;i<=n;i++){
        x=read();
        for(int j=1;j<=x;j++)b[j].x=read(),b[j].y=read();
        b[x+1]=b[1];
        for(int j=1;j<=x;j++)a[++m]=(Line){b[j],b[j+1]};
    }
    for(int i=1;i<=m;i++)a[i].get();
    sort(a+1,a+m+1);
    for(int i=1;i<=m;i++){
        if(abs(a[i].d-a[i-1].d)>eps)++cnt;
        c[cnt]=a[i];//去重
    }
    m=cnt;q[l=r=1]=1;q[++r]=2;
    for(int i=3;i<=m;i++){
        while(l<r&&check(c[q[r-1]],c[q[r]],c[i]))r--;
        while(l<r&&check(c[q[l]],c[q[l+1]],c[i]))l++;
        q[++r]=i;
    }
    while(l<r&&check(c[q[r-1]],c[q[r]],c[q[l]]))r--;
    // while(l<r&&check(c[q[l]],c[q[l+1]],c[q[r]]))l++;
    q[r+1]=q[l];
    for(int i=l;i<=r;i++)b[i]=inter(c[q[i]],c[q[i+1]]);
    b[r+1]=b[l];
    for(int i=l;i<=r;i++)ans+=b[i]^b[i+1];ans/=2;
    if(abs(ans)<eps)ans=0;
    printf("%.3f\n",ans);
    return 0;
}

半平面交还是应该用凸包来理解。因为交出的一定是凸多边形，所以我们按照斜率排序之后维护这个凸包就行了。最初论文中的写法是分别求出上下凸包，然后合并时发现会有无用边
于是需要互相更新一下去掉无用边。
而现在我们见到的算法是经过优化的，我们不再分别求上下凸包，而是整体求一遍。（我们设求凸包的顺序是从 $-\pi$ 开始逆时针求的..）这样我们发现右上和右下的无用边已经被去掉了。然后我们只要在加入边的过程中顺便把左下（队列头）的无用边也顺便去掉，那么剩下的只有左上的无用边，只需最后用队头更新队尾即可。

没见过半平面交可能会是无穷的题目。但是半平面交是空集的还是有可能的。上面这个代码如果是空集的话会输出nan QAQ。怎么判断是不是空集？当已经确定交集一定不是无穷之后，交集为空当且仅当最后队列中的元素小于等于 $2$ 。
（为了防止是无穷的情况，可以先在外面套一个长方形大框）

旋转卡壳

这个是求凸包直径，这真没什么好说的。。就是求完凸包之后双指针扫

#include<bits/stdc++.h>
using namespace std;
const int N=1e5+50;
int n,q[N],m,ans;
struct Poi{
	int x,y;
	bool friend operator <(Poi a,Poi b){return a.x!=b.x?a.x<b.x:a.y<b.y;}
	int operator *(const Poi &b){return x*b.y-y*b.x;}
	Poi operator -(const Poi &b){return (Poi){x-b.x,y-b.y};}
}a[N],b[N];
int cross(Poi a,Poi b,Poi c){return (b-a)*(c-a);}
int dis(Poi a,Poi b){return (b.x-a.x)*(b.x-a.x)+(b.y-a.y)*(b.y-a.y);}
int main(){
	scanf("%d",&n);
	for(int i=1;i<=n;i++)scanf("%d%d",&a[i].x,&a[i].y);
	sort(a+1,a+n+1);q[m=1]=1;
	for(int i=2;i<=n;q[++m]=i++)
		while(m>1&&cross(a[q[m-1]],a[q[m]],a[i])<=0)m--;
	for(int i=n,l=m;i;q[++m]=i--)
		while(m>l&&cross(a[q[m-1]],a[q[m]],a[i])<=0)m--;
	m--;for(int i=0;i<m;i++)b[i]=a[q[i+1]];
	for(int i=0,j=i+1;i<m;i++){
		while((b[(i+1)%m]-b[i])*(b[(j+1)%m]-b[j%m])>0)j++;
		ans=max(ans,dis(b[i],b[j%m]));
	}
	printf("%d\n",ans);
	return 0;
}

最小圆覆盖

只要知道一个定理：
若当前点集的最小覆盖圆为S且新加入的点不被S覆盖，则此点在新圆上。
于是我们只要十分naive地做就好了。具体看注释。还有，外接圆好好求啊233
还有复杂度，最内层循环是 $O(n)$ 的，外层循环有 $\frac{3}{n}$ 的几率进入内层循环。故总复杂度为 $O(n)$ 。
还有就是这个是不会出现三点共线求外接圆的情况的，具体的自己在纸上模拟几种情况即可。

#include<bits/stdc++.h>
#define double long double
using namespace std;
const int N=1e5+50;
const double eps=1e-9;
int n;double r;
struct poi{
    double x,y;
    poi operator +(const poi &b){return (poi){x+b.x,y+b.y};}
    poi operator -(const poi &b){return (poi){x-b.x,y-b.y};}
    double operator *(const poi &b){return x*b.y-y*b.x;}
    poi operator *(const double &b){return (poi){x*b,y*b};}
    poi operator ~(){return (poi){-y,x};}//向量旋转90度
}a[N],O;
double dis(poi a,poi b){
    double xx=a.x-b.x,yy=a.y-b.y;
    return sqrt(xx*xx+yy*yy);
}
bool out(poi x){return dis(O,x)>r+eps;}
poi inter(poi s1,poi t1,poi s2,poi t2){
    poi v1=t1-s1,v2=t2-s2;
    return s2+v2*((s2-s1)*v1/(v1*v2));
}
void circle(poi a,poi b,poi c){
    poi xx=(a+b)*0.5,yy=(b+c)*0.5;
    O=inter(xx+~(a-xx),xx+~(b-xx),yy+~(b-yy),yy+~(c-yy));//中垂线交点。
    r=dis(O,a);
}
int main(){
    scanf("%d",&n);
    for(int i=1;i<=n;i++)scanf("%Lf%Lf",&a[i].x,&a[i].y);
    random_shuffle(a+1,a+n+1);O=a[1];r=0;
    for(int i=1;i<=n;i++)if(out(a[i])){//外层循环为求前i个点的最小覆盖圆。
        O=a[i];r=0;//因为这个点一定在圆上。
        for(int j=1;j<i;j++)if(out(a[j])){//第二层循环为求前j个点和第i个点的最小覆盖圆。
            O=(a[i]+a[j])*0.5;r=dis(a[i],a[j])/2;//因为这两个点都一定在圆上。
            for(int k=1;k<j;k++)if(out(a[k])){//这层...
                circle(a[i],a[j],a[k]);//这三个点都一定在圆上。
            }
        }
    }//完事233
    printf("%.10Lf\n%.10Lf %.10Lf",r,O.x,O.y);
    return 0;
}

辛普森积分

让你求一个要多奇怪有多奇怪的函数的积分。
首先我们知道二次函数的积分满足这个：

\int_l^rf(x)dx=\frac{m}{3}(f(l)+f(r)+4f(m))

其中 $m=\frac{l+r}{2}$ 。推导可以见洛谷题解，就是拆一波再合一波。
这东西只用求 $O(1)$ 次原函数即可。所以可以考虑用二次函数拟合原函数。
然后我们分几段呢？分得多会T，分得少会WA。那么我们可以让程序自己判断是否够精度，不够就继续往下分。其实原理非常简单啦！

#include<bits/stdc++.h>
using namespace std;
const double eps=1e-7;
double a,b,c,d,l,r;
double f(double x){return (c*x+d)/(a*x+b);}//要积分的原函数
double simpson(double l,double r){return (r-l)*(f(l)+f(r)+4*f((l+r)/2))/6;}//
double calc(double l,double r,double eps,double ans){
    double mid=(l+r)/2;
    double ll=simpson(l,mid),rr=simpson(mid,r);
    if(abs(ll+rr-ans)<eps*15)return ll+rr;
    return calc(l,mid,eps/2,ll)+calc(mid,r,eps/2,rr);
}
int main(){
    scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&l,&r);
    printf("%.6f",calc(l,r,eps,simpson(l,r)));
    return 0;
}

至于*15是为什么？

大概是这样可以还原真实的误差（听说是实验结果？）这篇文章叫AdaptiveQuadProof。

然后就有一个特别经典的例题：月下柠檬树
题意是有一棵树，它是好多台体和顶上的锥体堆起来的，给出月光的角度，求这棵树投影的面积。
首先因为是平行投影，所以相当于是把所有圆平移了，然后因为圆与圆之间是平滑地相接的，所以会有公切线（相当于梯形）。
面积是关于x轴对称的，只要求上方的面积就好了。求面积又让人想到积分，即使这函数是分了无数段233，不过它还是连续的，这就可以积分了。
于是求出所有的公切线（用几何性质导角什么的），然后一个点的函数值就是定义域包含它的所有函数的最大值，直接暴力扫取max即可，然后套板子，其实特别短。

#include<bits/stdc++.h>
#define y1 orzjackpei
using namespace std;
const int N=510;
const double eps=1e-6;
int n;double alpha,x[N],R[N],l[N],r[N],y1[N],y2[N],mx,mn;
double f(double X){
    double Y=0;
    for(int i=1;i<n;i++)if(X>x[i]-R[i]+eps&&X<x[i]+R[i]-eps)
        Y=max(Y,sqrt(R[i]*R[i]-(X-x[i])*(X-x[i])));
    for(int i=1;i<n;i++)if(X>l[i]-eps&&X<r[i]+eps)
        Y=max(Y,y1[i]+(y2[i]-y1[i])*(X-l[i])/(r[i]-l[i]));
    return Y;
}
double simpson(double l,double r){return (r-l)*(f(l)+f(r)+4*f((l+r)/2))/6;}
double calc(double l,double r,double eps,double A){
    double mid=(l+r)/2,ll=simpson(l,mid),rr=simpson(mid,r);
    if(abs(ll+rr-A)<eps*15)return ll+rr;
    return calc(l,mid,eps/2,ll)+calc(mid,r,eps/2,rr);
}
int main(){
    scanf("%d%lf",&n,&alpha);alpha=1/tan(alpha);n++;
    for(int i=1;i<=n;i++)scanf("%lf",&x[i]),x[i]+=x[i-1];
	for(int i=1;i<=n;i++)x[i]*=alpha;
    for(int i=1;i<n;i++)scanf("%lf",&R[i]),mx=max(mx,x[i]+R[i]),mn=min(mn,x[i]-R[i]);
    mx=max(mx,x[n]);
    for(int i=1;i<n;i++){
        double x1=x[i],x2=x[i+1],r1=R[i],r2=R[i+1];
        if(x2-x1<abs(r1-r2)+eps){l[i]=-1e9,r[i]=-1e9;continue;}
        double s=abs(r1-r2)/(x2-x1),c=sqrt(1-s*s);
        if(r2<r1)l[i]=x1+r1*s,r[i]=x2+r2*s;
        else l[i]=x1-r1*s,r[i]=x2-r2*s;
        y1[i]=r1*c;y2[i]=r2*c;
    }
    printf("%.2f",calc(mn,mx,eps,simpson(mn,mx))*2);
    return 0;
}

有一个要注意的，就是当函数有大段为0（或者说为一个相同值）的时候需要把所有连续的不为0的段分别积分。如求圆面积并时这么一个数据：

2
0 0 1
100 100 1

直接求会使程序认为这个函数就是一个常为0的函数，从而积出0。
还有，目前已知除了求圆并以外别的图形都有更高效的算法，且辛普森积分确实会有精度问题，不得已再用吧qwq

Prutekoi's blog

半平面交

旋转卡壳

最小圆覆盖

辛普森积分

线性规划，单纯形和对偶问题