上传一份板子,QAQ

上传一份板子,QAQ

前言

可能这是最后的一份板子了。里面能说比较有用的可能就是一些数学上的东西吧,不过很多常见的推导过程都没有写,群论到最后也没学,建议后面的同学(pxa)还是去看一下这样很多数论方面的东西才能找到第一步,FFT什么的建议搭配一般型的母函数构造进行食用,洲阁筛还没学,不过应用据说还不如min_25,复杂度不如杜教筛。几何是抄的任板的2D几何,后面还放了份偷来的吉爷爷的板子(里面有3D的部分)。图论什么的还是建议多去想想一些奇怪的建图方式,很多时候图建出来就做完了。
其实暂时还没想好后面再干啥 。有没有做项目的童鞋带带啊QAQ 不过应该会参加leetcode上的周赛什么的。

数据结构
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;//快读
inline long long read()
{
	long long kk=0,f=1;
	char cc=getchar();
	while(cc<'0'||cc>'9'){if(cc=='-')f=-1;cc=getchar();}
	while(cc>='0'&&cc<='9'){kk=(kk<<1)+(kk<<3)+cc-'0';cc=getchar();}
	return kk*f;
}
Trie与异或

给出数组a求max(a[i]^(a[j]+a[k]))其中i,j,k不相等。

LL trie[33*maxn][2],tot=1,n,a[maxn];
int cnt[33*maxn];
void insert(LL now)
{
	int p=1;
	for(int i=33;i>=0;--i)
	{
		int ch=(now>>i)&1;
		if(!trie[p][ch])trie[p][ch]=++tot;
		p=trie[p][ch];cnt[p]++;
	}
}
void delet(LL now)
{
	int p=1;
	for(int i=33;i>=0;--i)
	{
		int ch=(now>>i)&1;
		p=trie[p][ch];cnt[p]--;
	}
}
LL search(LL now)
{
	int p=1;LL asd=0;
	for(int i=33;i>=0;--i)
	{
		int ch=(now>>i)&1;
		if(trie[p][ch^1]&&cnt[trie[p][ch^1]])
		{
			p=trie[p][ch^1];
			asd|=(1ll<<i);
		}
		else {p=trie[p][ch];}
	}
	return asd;
}
void solve()
{
	sfl(n);
	memset(trie,0,sizeof(trie));
	memset(cnt,0,sizeof(cnt));
	tot=1;
	for(int i=1;i<=n;++i){sfl(a[i]);insert(a[i]);}
	LL ans=0;
	for(int i=1;i<=n;++i)
	{
		delet(a[i]);
		for(int j=i+1;j<=n;++j)
		{
			delet(a[j]);
			ans=max(ans,search(a[i]+a[j]));
			insert(a[j]);
		}
		insert(a[i]);
	}
	printf("%lld\n",ans);
}
ST表
const LL N=2000222;//st表
int a[N];int n,m;
int f[N][30];
void pre_st()
{
	for(int i=1;i<=n;++i)f[i][0]=a[i];
	int t=log2(n)+1;
	for(int j=1;j<t;++j)
	for(int i=1;i<=(n-(1<<j)+1);++i)
	{
		f[i][j]=max(f[i][j-1],f[i+(1<<(j-1))][j-1]);
	}
}
int query_st(int l,int r)
{
	int len=(log2(r-l+1));
	return max(f[l][len],f[r-(1<<(len))+1][len]);
}
线段树
const int N=100022;//线段树
LL ans[N];int n,m;
struct segment_tree
{
	int l,r;LL dat,add;
	#define l(x) tree[x].l
	#define r(x) tree[x].r
	#define dat(x) tree[x].dat
	#define add(x) tree[x].add
}tree[4*N];
void update(int p){dat(p)=dat(p*2)+dat(p*2+1);}
void build(int p,int l,int r)
{
	l(p)=l;r(p)=r;
	if(l==r){dat(p)=ans[l];return;}
	int mid=l+(r-l)/2;
	build(p*2,l,mid);build(p*2+1,mid+1,r);
	update(p);
}
void spread(int p)
{
	if(add(p))
	{
		dat(p*2)+=(r(p*2)-l(p*2)+1)*add(p);
		dat(p*2+1)+=(r(p*2+1)-l(p*2+1)+1)*add(p);
		add(p*2)+=add(p);
		add(p*2+1)+=add(p);
		add(p)=0;
	}
}
void change(int p,int l,int r,int v)
{
	if(l<=l(p)&&r(p)<=r){dat(p)+=(LL)v*(r(p)-l(p)+1);add(p)+=v;return;}
	spread(p);
	int mid=l(p)+(r(p)-l(p))/2;
	if(l<=mid)change(p*2,l,r,v);
	if(r>mid)change(p*2+1,l,r,v);
	update(p);
}
LL query(int p,int l,int r)
{
	if(l<=l(p)&&r(p)<=r)return dat(p);
	spread(p);
	int mid=l(p)+(r(p)-l(p))/2;
	LL val=0;
	if(l<=mid)
	{
		val+=query(p*2,l,r);
	}
	if(r>mid)
	{
		val+=query(p*2+1,l,r);
	}
	return val;
}
最大权值和矩阵、线段树维护最大子段和
int n,m;
struct Point
{
	int x,y;LL val;
	bool operator<(const Point b)const
	{
		if(y==b.y)return x<b.x;
		return y>b.y;
	}
}node[maxn];
int hashx[maxn],hashy[maxn];
struct Segmentree
{
	LL dat,ldat,rdat,sum;int l,r;
	#define l(p) tre[p].l
	#define r(p) tre[p].r
	#define dat(p) tre[p].dat
	#define sum(p) tre[p].sum
	#define ldat(p) tre[p].ldat
	#define rdat(p) tre[p].rdat
}tre[maxn<<2];
void pushup(int p)
{
	sum(p)=sum(p<<1)+sum(p<<1|1);
	ldat(p)=max(ldat(p<<1),sum(p<<1)+ldat(p<<1|1));
	rdat(p)=max(rdat(p<<1|1),sum(p<<1|1)+rdat(p<<1));
	dat(p)=max(rdat(p<<1)+ldat(p<<1|1),max(dat(p<<1|1),dat(p<<1)));
	dat(p)=max(dat(p),max(sum(p),max(ldat(p),rdat(p))));
}
void build(int p,int l,int r)
{
	l(p)=l;r(p)=r;dat(p)=ldat(p)=rdat(p)=sum(p)=0;
	if(l==r){return;}
	int mid=(l(p)+r(p))/2;
	build(p<<1,l,mid);build(p<<1|1,mid+1,r);
	pushup(p);
}
void addtre(int p,int x,LL val)
{
	if(l(p)==r(p)&&l(p)==x)
	{
		sum(p)=sum(p)+val;
		dat(p)=ldat(p)=rdat(p)=sum(p);
		return;
	}
	int mid=(l(p)+r(p))/2;
	if(x<=mid)addtre(p<<1,x,val);
	else addtre(p<<1|1,x,val);
	pushup(p);
}
int main()
{
	int _;sfi(_);
	while(_--)
	{
		sfi(n);
		for(int i=1;i<=n;++i)
		{
			sfi(node[i].x);sfi(node[i].y);sfl(node[i].val);
			hashx[i]=node[i].x;hashy[i]=node[i].y;
		}
		sort(hashx+1,hashx+1+n);
		sort(hashy+1,hashy+1+n);
		for(int i=1;i<=n;++i)
		{
			node[i].x=lower_bound(hashx+1,hashx+1+n,node[i].x)-hashx;
			node[i].y=lower_bound(hashy+1,hashy+1+n,node[i].y)-hashy;
		}
		sort(node+1,node+1+n);
		LL ans=-1e18;
		for(int i=1;i<=n;++i)
		{
			if(i!=1&&node[i].y==node[i-1].y)continue;
			build(1,1,n);
			for(int j=i;j<=n;)
			{
				int k=j;
				while(k<=n&&node[j].y==node[k].y)
				{
					addtre(1,node[k].x,node[k].val);
					k++;
				}
				j=k;
				ans=max(ans,dat(1));
			}
		}
		printf("%lld\n",ans);
	}
}
吉老师线段树

modify 将区间大于x的数变成x
query 询问区间和
单次复杂度 O(log^2(n))

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

/*
    modify 将区间大于x的数变成x
    query 询问区间和
    单次复杂度 O(log^2(n))
*/

const ll INF=0xc0c0c0c0c0c0c0c0ll;
const int MAXN=200005;
ll seg[MAXN<<2],m1[MAXN<<2],m2[MAXN<<2],cnt[MAXN<<2],tag[MAXN<<2],a[MAXN];
int n,q;

void pushdown(int rt)
{
    if(!tag[rt]) return;
    ll y=m1[rt];
    if(y<m1[rt<<1])
    {
        tag[rt<<1]=1;
        seg[rt<<1]-=(m1[rt<<1]-y)*cnt[rt<<1];
        m1[rt<<1]=y;
    }
    if(y<m1[rt<<1|1])
    {
        tag[rt<<1|1]=1;
        seg[rt<<1|1]-=(m1[rt<<1|1]-y)*cnt[rt<<1|1];
        m1[rt<<1|1]=y;
    }
    tag[rt]=0;
}

void pushup(int rt)
{
    seg[rt]=seg[rt<<1]+seg[rt<<1|1];
    if(m1[rt<<1]==m1[rt<<1|1])
    {
        m1[rt]=m1[rt<<1];
        cnt[rt]=cnt[rt<<1]+cnt[rt<<1|1];
        m2[rt]=max(m2[rt<<1],m2[rt<<1|1]);
    }
    else if(m1[rt<<1]>m1[rt<<1|1])
    {
        m1[rt]=m1[rt<<1];
        cnt[rt]=cnt[rt<<1];
        m2[rt]=max(m2[rt<<1],m1[rt<<1|1]);
    }
    else
    {
        m1[rt]=m1[rt<<1|1];
        cnt[rt]=cnt[rt<<1|1];
        m2[rt]=max(m2[rt<<1|1],m1[rt<<1]);
    }
}

void build(int rt,int l,int r)
{
    tag[rt]=0;
    if(l==r)
    {
        seg[rt]=m1[rt]=a[l];
        cnt[rt]=1;
        m2[rt]=INF;
        return;
    }
    int m=l+r>>1;
    if(l<=m) build(rt<<1,l,m);
    if(m<r) build(rt<<1|1,m+1,r);
    pushup(rt);
}

void modify(int rt,int l,int r,int L,int R,ll y)
{
    if(y>=m1[rt]) return;
    if(L<=l&&r<=R&&y>m2[rt])
    {
        tag[rt]=1;
        seg[rt]-=(m1[rt]-y)*cnt[rt];
        m1[rt]=y;
        return;
    }
    pushdown(rt);
    int m=l+r>>1;
    if(L<=m) modify(rt<<1,l,m,L,R,y);
    if(m<R) modify(rt<<1|1,m+1,r,L,R,y);
    pushup(rt);
}

ll query(int rt,int l,int r,int L,int R)
{
    if(L<=l&&r<=R) return seg[rt];
    int m=l+r>>1;
    pushdown(rt);
    ll ret=0;
    if(L<=m) ret+=query(rt<<1,l,m,L,R);
    if(m<R) ret+=query(rt<<1|1,m+1,r,L,R);
    pushup(rt);
    return ret;
}
扫描线之求两个以上的矩形交的面积
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstdlib>
#include<cstring>
using namespace std;
typedef long long LL;
const int maxn=2022;
#define sfi(a) scanf("%d",&a)
#define sfl(a) scanf("%lld",&a)
#define sff(a) scanf("%lf",&a)
struct Line
{
	double x;
	double rupy,rdowny;
	int upy,downy;
	int flag;
	bool operator<(const Line b)const
	{
		return x<b.x;
	}
}line[maxn<<1];
struct Segmentree
{
	int cov=0;
	int l,r;
	double dat1,dat2;
	#define l(p) tre[p].l
	#define r(p) tre[p].r
	#define cov(p) tre[p].cov
	#define dat1(p) tre[p].dat1
	#define dat2(p) tre[p].dat2
}tre[maxn<<3];
LL ix[maxn<<1];int n,m;
double ry[maxn<<1];

void pushup(int p)
{
	if(cov(p))
	{
		dat1(p)=ry[r(p)+1]-ry[l(p)];
		if(cov(p)>=2)
		{
			dat2(p)=ry[r(p)+1]-ry[l(p)];
		}
		else if(cov(p)==1)
		{
			dat2(p)=dat1(p<<1)+dat1(p<<1|1);
		}
	}
	else 
	{
		dat1(p)=dat1(p<<1)+dat1(p<<1|1);
		dat2(p)=dat2(p<<1)+dat2(p<<1|1);
	}
}
void build(int p,int l,int r)
{
	l(p)=l;r(p)=r;
	if(l==r){return;}
	int mid=(l(p)+r(p))/2;
	build(p<<1,l,mid);build(p<<1|1,mid+1,r);
}
void addtre(int p,int l,int r,int val)
{
	if(l<=l(p)&&r(p)<=r)
	{
		cov(p)+=val;
		pushup(p);
		return;
	}
	int mid=(l(p)+r(p))/2;
	if(mid>=l)addtre(p<<1,l,r,val);
	if(mid<r)addtre(p<<1|1,l,r,val);
	pushup(p);
}

void solve()
{
	sfi(n);
	memset(tre,0,sizeof(tre));
	for(int i=1;i<=n;++i)
	{
		double fx,fy,sx,sy;
		sff(fx);sff(fy);sff(sx);sff(sy);
		
		line[i].rupy=max(fy,sy);line[i].rdowny=min(fy,sy);
		line[i].x=min(fx,sx);line[i].flag=1;
		
		line[i+n].rupy=max(fy,sy);line[i+n].rdowny=min(fy,sy);
		line[i+n].x=max(fx,sx);line[i+n].flag=-1;
		
		ry[i]=fy;ry[i+n]=sy;
	}
	
	sort(line+1,line+1+2*n);
	sort(ry+1,ry+1+2*n);
	build(1,1,2*n);
	for(int i=1;i<=2*n;++i)
	{
		line[i].upy=lower_bound(ry+1,ry+1+2*n,line[i].rupy)-ry;
		line[i].downy=lower_bound(ry+1,ry+1+2*n,line[i].rdowny)-ry;
	}
	double ans=0.00;
	for(int i=1;i<=2*n;++i)
	{
		if(i!=1)ans+=dat2(1)*(line[i].x-line[i-1].x);
		addtre(1,line[i].downy,line[i].upy-1,line[i].flag);
	}
	printf("%.2lf\n",ans+1e-6);
}
int main()
{
	int _;sfi(_);
	while(_--)solve();
	return 0;
}
数学
矩阵快速幂
struct matrix
{
	LL jie[5][5];
	int n,m;
	matrix operator * (const matrix &a)const
	{
		matrix b;b.n=n;b.m=a.m;
		for(int i=1;i<=n;++i)
		for(int j=1;j<=a.m;++j)
		{
			b.jie[i][j]=0;
			for(int k=1;k<=m;++k)
			{
				b.jie[i][j]=((jie[i][k]*a.jie[k][j])+b.jie[i][j]);
			}
		}
		return b;
	}
}t1,t2,t3;
matrix mpow(matrix a,LL b)
{
	matrix kk;kk.n=kk.m=a.n;
	for(int i=1;i<=kk.n;++i)kk.jie[i][i]=1;
	while(b)
	{
		if(b&1)kk=kk*a;
		a=a*a;
		b>>=1;
	}
	return kk;
}
01分数规划
LL a[maxn],b[maxn],n,k;//01分数规划
bool check(double L)
{
	vector<double>V;
	for(int i=1;i<=n;++i)V.push_back((double)((double)a[i]-L*b[i]));
	sort(V.begin(),V.end(),greater<double>());
	double asd=0;
	for(int i=0;i<k;++i)asd+=V[i];
	if(asd>1e-8)return 1;
	return 0;
}
void solve()
{
	n=read();k=read();
	for(int i=1;i<=n;++i){b[i]=read();a[i]=read();}
	double l,r;l=0;r=1e5;
	while(abs(r-l)>1e-4){double mid=(l+r)/2.0;if(check(mid))l=mid;else r=mid;}
	printf("%.2lf\n",l);
}
扩展中国剩余(excrt)、exgcd、同余方程、逆元
LL k[1000222],m[2000222],n;
LL exgcd(LL a,LL b,LL &x,LL &y)
{
	if(!b){x=1;y=0;return a;}
	LL lin=exgcd(b,a%b,y,x);
	y-=(a/b)*x;
	return lin;
}
LL mmul(LL a,LL b,LL mod)
{
	LL kk=0;
	while(b)
	{
		if(b&1)kk=(kk+a)%mod;
		a=(a+a)%mod;
		b>>=1;
	}
	return kk;
}
LL excrt()
{
	LL x,y,a,b,c,M,ans,g;
	M=m[1];ans=k[1];
	for(int i=2;i<=n;++i)
	{
		a=M;b=m[i];
		c=((k[i]-ans)%b+b)%b;
		g=exgcd(a,b,x,y);
		if(c%g!=0)return -1;
		x=mmul(x,c/g,b/g);
		ans=(ans+M*x);M*=(b/g);
		ans=(ans%M+M)%M;
	}
	return (ans%M+M)%M;
}
LL ni(LL a)
{
   LL nia,y;
   exgcd(a,mod,nia,y);
   return (nia%mod+mod)%mod;
}
LL inv[MAXN];
void init()//线性预处理O(P)
{
	inv[1]=1;
	for(int i=2;i<mod;i++)
	{
		inv[i]=(mod-(mod/i))*inv[mod%i]%mod;
	}
}
bool tyfc(LL a,LL b,LL c)//求ax+by==c
{
	LL x,y;
	LL gcdn=exgcd(a,b,x,y);
	if(c%gcdn)return 0;
	x*=(c/gcdn);y*=(c/gcdn);//初始解
	LL mmod=b/gcdn;//通解x的增量
	x=(x%mmod+mmod)%mmod;
	y=(c-a*x)/b;
	if(y<0)return 0;//失败
	while(1)//所有正解
	{
		y=y-a/gcdn;
		if(y<0)break;
		x=x+b/gcdn;
	}
	return 1;
}
莫比乌斯、整除分块与欧拉函数
LL su[maxn],co=0,u[maxn],gao[maxn],sum[maxn],phi[maxn];

bool vissu[maxn];
void init()
{
	u[1]=1;phi[1]=1;
	for(int i=2;i<maxn;++i)
	{
		if(!vissu[i])
		{
			su[++co]=i;phi[i]=i-1;
			u[i]=-1;gao[i]=1;
		}
		for(int j=1;j<=co&&i*su[j]<maxn;++j)
		{
			vissu[i*su[j]]=1;
			if(i%su[j])
			{
				u[i*su[j]]=-u[i];
				gao[i*su[j]]=u[i]-gao[i];//根据题目将某些式子优化成可以预处理的情况
				phi[i*su[j]]=phi[i]*(su[j]-1);
			}
			else
			{
				u[i*su[j]]=0;
				gao[i*su[j]]=u[i];
				phi[i*su[j]]=phi[i]*(su[j]);
				break;
			}
		}
	}
	for(int i=1;i<maxn;++i)
	{
		sum[i]=sum[i-1]+gao[i];
	}
}
void solve()
{
	int n,m;
	n=read();m=read();
	LL asd=0;
	for(int i=1,la;i<=min(n,m);i=la+1)//整除分块
	{
		la=min(n/(n/i),m/(m/i));
		asd+=(sum[la]-sum[i-1])*(n/i)*(m/i);
	}
	printf("%lld\n",asd);
}
常见卷积

μ ( i ) \mu(i) μ(i)莫比乌斯函数
φ ( i ) \varphi(i) φ(i)欧拉函数
i d ( n ) id(n) id(n)单位函数 i d ( n ) = = n id(n)==n id(n)==n
$\sigma(n) $因子和
i d k ( n ) idk(n) idk(n)幂函数 i d k ( n ) = = n k idk(n)==n^k idk(n)==nk
ε ( i ) \varepsilon (i) ε(i)元函数 ε ( i ) = = [ i = = 1 ] \varepsilon (i)==[i==1] ε(i)==[i==1]
性质:
1.若 n = ∑ i = 1 n p i a i n=\sum_{i=1}^n{p_i^{a_i}} n=i=1npiai,那么 f ( n ) = ∏ i = 1 n f ( p i a i ) f(n)=\prod_{i=1}^{n}{f(p_i^{a_i})} f(n)=i=1nf(piai)
2.若f为积性函数且有 f ( p i ) = = f i ( p ) f(p^{i})==f^{i}(p) f(pi)==fi(p)那么f为完全积性函数。
3. μ ∗ i d = = φ \mu*id==\varphi μid==φ
4. φ ∗ I = = i d \varphi*I==id φI==id
5. I ∗ μ = = ε I*\mu==\varepsilon Iμ==ε
6. σ = I ∗ i d = I ∗ ( I ∗ φ ) = ( I ∗ I ) ∗ φ = τ ∗ φ σ=I∗id=I∗(I∗φ)=(I∗I)∗φ=\tau∗φ σ=Iid=I(Iφ)=(II)φ=τφ

杜教筛

为了求积性函数f前缀和s即: S ( n ) = ∑ i = 1 n f [ i ] S(n)=\sum_{i=1}^n{f[i]} S(n)=i=1nf[i]
构造两个积性函数h和g。使得h=f∗g
g ( 1 ) S ( n ) = ∑ i = 1 n h ( i ) − ∑ d = 2 n g [ d ] S ( ⌊ n d ⌋ ) g(1)S(n)=\sum_{i=1}^n{h(i)}−\sum_{d=2}^n{g[d]S(⌊\frac{n}{d}⌋)} g(1)S(n)=i=1nh(i)d=2ng[d]S(dn)
当你的h(i)的前缀和很好求,能在较短的时间内求出,那么我们对后面的式子进行整除分块。
f ( i ) = = μ ( i ) f(i)==\mu(i) f(i)==μ(i)则可让 g ( i ) = = 1 ( g = = I ) g(i)==1(g==I) g(i)==1(g==I)那么 S ( n ) = 1 − ∑ d = 2 n S ( ⌊ n d ⌋ ) S(n)=1−\sum_{d=2}^n{S(⌊\frac{n}{d}⌋)} S(n)=1d=2nS(dn)
f ( i ) = = φ ( i ) f(i)==\varphi(i) f(i)==φ(i)
我们在脑海中找到一个与欧拉函数有关的卷积式子: φ ∗ I = i d \varphi∗I=id φI=id
则可让 g ( i ) = = 1 ( g = = I ) g(i)==1(g==I) g(i)==1(g==I)那么 S ( n ) = ∑ i = 1 n i − ∑ d = 2 n S ( ⌊ n d ⌋ ) S(n)=\sum_{i=1}^n{i}−\sum_{d=2}^n{S(⌊\frac{n}{d}⌋)} S(n)=i=1nid=2nS(dn)
f ( i ) = = i ∗ φ ( i ) f(i)==i*\varphi(i) f(i)==iφ(i)
考虑卷积式 ∑ d ∣ n ( d ⋅ φ ( d ) ) ⋅ g ( n d ) \sum_{d|n}(d⋅\varphi(d))⋅g(\frac{n}{d}) dn(dφ(d))g(dn)
尝试将g配成id。这样就可以把d给弄没
则变成 n ∑ d ∣ n ( φ ( d ) ) ) < = > n φ ∗ I < = > n 2 n\sum_{d|n}(\varphi(d)))<=>n\varphi*I<=>n^2 ndn(φ(d)))<=>nφI<=>n2
那么 S ( n ) = ∑ i = 1 n i 2 − ∑ d = 2 n d S ( ⌊ n d ⌋ ) S(n)=\sum_{i=1}^n{i^2}−\sum_{d=2}^n{dS(⌊\frac{n}{d}⌋)} S(n)=i=1ni2d=2ndS(dn)

以求 f ( i ) = = φ ( i ) f(i)==\varphi(i) f(i)==φ(i) f ( i ) = = μ ( i ) f(i)==\mu(i) f(i)==μ(i)为例

#include<bits/stdc++.h>
//#pragma GCC optimize(2)
using namespace std;
typedef long long LL;
const int MAXN=10000555;
const LL INF=1e9+7;
const LL MOD=1e9+7;
const double PI=acos(-1);
#define pii pair<int,int>
#define mp(x,y) make_pair(x,y)
LL su[MAXN],co=0,u[MAXN],phi[MAXN];
LL sumu[MAXN],sump[MAXN];
bool vis[MAXN];
void init()//小范围线性
{
	u[1]=1;phi[1]=1;
	for(int i=2;i<MAXN;++i)
	{
		if(!vis[i])
		{
			su[++co]=i;
			u[i]=-1;phi[i]=i-1;
		}
		for(int j=1;j<=co&&i*su[j]<MAXN;++j)
		{
			vis[i*su[j]]=1;
			if(i%su[j])
			{
				u[i*su[j]]=-u[i];
				phi[i*su[j]]=phi[i]*phi[su[j]];
			}
			else
			{
				u[i*su[j]]=0;
				phi[i*su[j]]=phi[i]*(su[j]);
				break;
			}
		}
	}
	for(int i=1;i<MAXN;++i)
	{
		sumu[i]=sumu[i-1]+u[i];
		sump[i]=sump[i-1]+phi[i];
	}
}
/*
LL solve(LL x)
{
    if(x<=MAXN-1)return f[x];
    if(F[x])return F[x];
	LL ans=H(x);
    for(long long l=2,r;l<=x;l=r+1)
        r=x/(x/l),ans-=(G(r)-G(l-1))*solve(x/l);
    return F[x]=ans;
}
*/
unordered_map<LL,LL>Sump,Sumu;
LL solveu(LL x)//大范围分块,记忆化
{
    if(x<=MAXN-1)return sumu[x];
    if(Sumu[x])return Sumu[x];
	LL ans=1;
    for(LL l=2,r;l<=x;l=r+1)
        r=x/(x/l),ans-=(r-l+1)*solveu(x/l);
    return Sumu[x]=ans;
}
LL solvep(LL x)
{
    if(x<=MAXN-1)return sump[x];
    if(Sump[x])return Sump[x];
	LL ans=(x+1)*x/2;
    for(LL l=2,r;l<=x;l=r+1)
        r=x/(x/l),ans-=(r-l+1)*solvep(x/l);
    return Sump[x]=ans;
}
void solve()
{
	LL n;scanf("%lld",&n);
	printf("%lld %lld\n",solvep(n),solveu(n));
}
int main()
{
	init();
	int _;scanf("%d",&_);
	while(_--)solve();
	return 0;
}
Min_25筛

显然杜教筛最大的限制就是要找到函数 g,h 使得 f * g = h,这样的函数并不好找。
Min_25筛用来求积性函数前缀和。范围一般是 1 0 11 10^{11} 1011
要求所求积性函数在f(i)已知时,f(p×i){p∈Primes, i×p∈[1,n]} 能快速计算。
因为多项式可以拆成若干个单项式,所以我们只需要考虑求出 f ( p ) = p k f(p)=p^k f(p)=pk的前缀和,然后每一项加起来就行了。
分成两部分,先考虑质数部分
∑ i = 1 n [ i 是 质 数 ] f [ i ] \sum_{i=1}^n{[i是质数]f[i]} i=1n[i]f[i]
g(n,j) 表示在前 j 个质数用埃氏筛筛完后的答案
P j 2 > n P_j^2>n Pj2>n,那么$ P_j$ 便对答案没有贡献了, g ( n , j ) = g ( n , j − 1 ) g(n,j)=g(n,j−1) g(n,j)=g(n,j1)
不然我们在 g(n,j−1) 的基础上要减掉一些东西。
P j 2 < = n P_j^2<=n Pj2<=n,那么 g ( n , j ) = g ( n , j − 1 ) − f ( P j ) × [ g ( n P j , j − 1 ) − ∑ i = 1 j − 1 f ( P i ) ] g(n,j)=g(n,j−1)-f(P_j)×[g(\frac{n}{P_j},j−1)−\sum_{i=1}^{j-1}{f(Pi)}] g(n,j)=g(n,j1)f(Pj)×[g(Pjn,j1)i=1j1f(Pi)]
再考虑合数部分
知道了质数的贡献,而目标函数又是积性函数,所以我们只需要用所有的质数拼出所有的合数就可以计算答案了。令S(n,j)表示所有最小质因子大于Pj的数i的f(i)的和
S ( n , j ) = g ( n , j ) − ∑ i = 1 j − 1 f ( P i ) + ∑ k = j P k 2 < = n ∑ e = 1 P k e + 1 ≤ n [ f ( P k e ) × S ( n P k e , k + 1 ) + f ( P k e + 1 ) ] S(n,j)=g(n,j)−\sum_{i=1}^{j−1}{f(P_i)}+\sum_{k=j}^{P_k^2<=n}{\sum_{e=1}^{P^{e+1}_k≤n}{[f(P^e_k)×S(\frac{n}{P^e_k},k+1)+f(P^{e+1}_k)]}} S(n,j)=g(n,j)i=1j1f(Pi)+k=jPk2<=ne=1Pke+1n[f(Pke)×S(Pken,k+1)+f(Pke+1)]
答案就是S(n,0)+f(1)。

以求 ∑ f [ x ] \sum{f[x]} f[x]其中 f ( p k ) = p k ∗ ( p k − 1 ) f(p^k)=p^k*(p^k-1) f(pk)=pk(pk1)为例

将式子拆分成两部分,一次项和二次项 f ( p k ) = p k ∗ p k − p k f(p^k)=p^k*p^k-p^k f(pk)=pkpkpk

#include<bits/stdc++.h>
//#pragma GCC optimize(2)
//f(p^k)=p^k*(p^k-1)
using namespace std;
typedef long long LL;
const int MAXN=10000555;
const LL MOD=1e9+7;
const LL inv3=333333336;
LL su[MAXN],co=0,num;
LL n,tot,w[MAXN],g1[MAXN],g2[MAXN],ind1[MAXN],ind2[MAXN];
//将题目中的式子拆成两项 1次项和2次项 
//ind1,ind2表示这个数在数组中的位置 
LL sum1[MAXN],sum2[MAXN];
bool vissu[MAXN]; 
LL S(LL x,LL y)
{
    if(su[y]>=x)return 0;
    LL k,ret;
    k=(x<=num?ind1[x]:ind2[n/x]);
	ret=(g2[k]-g1[k]-(sum2[y]-sum1[y])+MOD)%MOD;
    for(LL i=y+1;i<=co&&su[i]*su[i]<=x;i++)
    {
    	for(LL e=1,now=su[i],t;now*su[i]<=x;e++,now=now*su[i])
    	{
    		LL lin=now%MOD;
    		LL nex=now*su[i]%MOD;
    		ret=(ret+(lin)*(lin-1)%MOD*(S(x/now,i))%MOD+nex*(nex-1)%MOD)%MOD;
		}
	}
    return ret%MOD;
} 
void init()
{
	for(LL i=2;i<=num;i++)
	{
        if(!vissu[i])
		{
			su[++co]=i;
			sum1[co]=(sum1[co-1]+i)%MOD;//预处理递推质数时F的前缀和 
			sum2[co]=(sum2[co-1]+i*i)%MOD;
		}
        for(LL j=1;j<=co&&i*su[j]<=num;j++)
		{
            vissu[i*su[j]]=1;
            if(i%su[j]==0)break;
        }
    }
}
int main()
{
    scanf("%lld",&n);
    num=sqrt(n);
    init();
    for(LL l=1,r;l<=n;l=r+1)//预处理g函数j=0的时候的递推边界 
	{
        r=n/(n/l);
        w[++tot]=n/l;
        g1[tot]=w[tot]%MOD;
        g2[tot]=g1[tot]*(g1[tot]+1)/2%MOD;//此时g2是1~w[tot]的平方和 
		g2[tot]=g2[tot]*(2*g1[tot]+1)%MOD*inv3%MOD-1;
        
        g1[tot]=g1[tot]*(g1[tot]+1)/2%MOD-1;//此时是等差数列和 
        if(n/l<=num)ind1[n/l]=tot;
        else ind2[r]=tot;
    }
    for(LL i=1;i<=co;i++)//g函数的1~k的递推 
    for(LL j=1,k;j<=tot&&su[i]*su[i]<=w[j];j++)
	{
		if(w[j]/su[i]<=num)k=ind1[w[j]/su[i]];
		else k=ind2[n/(w[j]/su[i])];
        g1[j]-=su[i]*(g1[k]-sum1[i-1]+MOD)%MOD;
        g2[j]-=su[i]*su[i]%MOD*(g2[k]-sum2[i-1]+MOD)%MOD;
        g1[j]=(g1[j]%MOD+MOD)%MOD;
        g2[j]=(g2[j]%MOD+MOD)%MOD;
    }
    printf("%lld\n",(S(n,0)+1+MOD)%MOD);
    return 0;
}

例二

#include<bits/stdc++.h>
//#pragma GCC optimize(2)
//f(p^k)==p xor k and f(1)==1
using namespace std;
typedef long long LL;
const int MAXN=10000555;
const LL MOD=1e9+7;
const LL inv3=333333336;
LL su[MAXN],co=0,num;
LL n,tot,w[MAXN],g1[MAXN],ind1[MAXN],ind2[MAXN],g2[MAXN];
//ind1,ind2表示这个数在数组中的位置
LL sum1[MAXN],sum2[MAXN];
bool vissu[MAXN];
LL S(LL x,LL y)
{
    if(su[y]>=x||x<=1)return 0;
    LL k,ret;
    k=(x<=num?ind1[x]:ind2[n/x]);
	ret=(g1[k]-g2[k]-(sum1[y]-sum2[y])+MOD)%MOD;
	if(y==0)ret+=2;
    for(LL i=y+1;i<=co&&su[i]*su[i]<=x;i++)
    {
    	for(LL e=1,now=su[i],t;now*su[i]<=x;e++,now=now*su[i])
    	{
    		LL lin=now;
    		LL val=S(x/now,i);
    		ret=(ret+(su[i]^e)*val%MOD+(su[i]^(e+1))%MOD)%MOD;
		}
	}
    return ret%MOD;
}
void init()
{
	//sum1[1]=1;
	for(LL i=2;i<=num;i++)
	{
        if(!vissu[i])
		{
			su[++co]=i;
			sum1[co]=(sum1[co-1]+i)%MOD;//预处理递推质数时F的前缀和
			sum2[co]=(sum2[co-1]+1)%MOD;
		}
        for(LL j=1;j<=co&&i*su[j]<=num;j++)
		{
            vissu[i*su[j]]=1;
            if(i%su[j]==0)break;
        }
    }
}
int main()
{
    scanf("%lld",&n);
    num=sqrt(n);
    init();
    for(LL l=1,r;l<=n;l=r+1)//预处理g函数j=0的时候的递推边界
	{
        r=n/(n/l);
        w[++tot]=n/l;
        g1[tot]=w[tot]%MOD;
        g2[tot]=g1[tot]-1;
        g1[tot]=(g1[tot]*(g1[tot]+1)/2)%MOD-1	;
        if(n/l<=num)ind1[n/l]=tot;
        else ind2[r]=tot;
    }
    for(LL i=1;i<=co;i++)//g函数的1~k的递推
    for(LL j=1,k;j<=tot&&su[i]*su[i]<=w[j];j++)
	{
		if(w[j]/su[i]<=num)k=ind1[w[j]/su[i]];
		else k=ind2[n/(w[j]/su[i])];
        g1[j]-=(su[i])*(g1[k]-sum1[i-1]+MOD)%MOD;
        g2[j]-=(1)*(g2[k]-sum2[i-1]+MOD)%MOD;
        g2[j]=(g2[j]%MOD+MOD)%MOD;
        g1[j]=(g1[j]%MOD+MOD)%MOD;
    }
    printf("%lld\n",(S(n,0)+1+MOD)%MOD);
    return 0;
}

FFT
#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
const int MAXN=4*1e6+2e5;
const double PI=acos(-1.0);
struct complex
{
    double x,y;
    complex (double xx=0,double yy=0){x=xx,y=yy;}
	complex operator+(const complex &b)const{return (complex){x+b.x,y+b.y};}
	complex operator-(const complex &b)const{return (complex){x-b.x,y-b.y};}
	complex operator*(const complex &b)const{return (complex){x*b.x-y*b.y,x*b.y+y*b.x};}
}a[MAXN],b[MAXN];
int rev[MAXN],n,m;
int limit,bit;
void FFT(complex *A,int type)
{
    for(int i=0;i<limit;i++) 
    if(i<rev[i])swap(A[i],A[rev[i]]);
    for(int mid=1;mid<limit;mid<<=1)
    {
        complex Wn;
		Wn=(complex){cos(PI/(double)mid),type*sin(PI/(double)mid)};//单位根
        for(int R=(mid<<1),j=0;j<limit;j+=R)
        {
            complex w;
			w=(complex){1,0};//单位1
            for(int k=0;k<mid;k++,w=w*Wn)
            {
                complex x=A[j+k],y=w*A[j+mid+k];
                A[j+k]=x+y;//蝴蝶变换
                A[j+mid+k]=x-y;
            }
        }
    }
}
void gaofft()//a,b存初始系数,最后的系数放a里面
{
	limit=1;bit=0;
	while(limit<=n+m){limit<<=1;bit++;}
	for(int i=0;i<limit;++i)
	{
		rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
	}
    FFT(a,1);//a做一次逆变换,
    FFT(b,1);//从系数变成点权
    for(int i=0;i<=limit;i++)
    a[i]=a[i]*b[i];
    FFT(a,-1);//从点权变成系数
    for(int i=0;i<=n+m;++i)
    {
    	a[i].x=(int)(a[i].x/limit+0.5);
	}
}
int main()
{
    scanf("%d",&n);scanf("%d",&m);
    for(int i=0;i<=n;i++)
    {
    	scanf("%lf",&a[i].x);
	}
    for(int i=0;i<=m;i++)
	{
		scanf("%lf",&b[i].x);
	}
	gaofft();
    for(int i=0;i<=n+m;i++)printf("%.0f ",a[i].x);
    return 0;
}
博弈
博弈打表
int SG[maxn][maxn];
void get()//以威佐夫博弈为例 打印必败
{
    int n=20;
    for(int i=0;i<=n;i++)
    for(int j=0;j<=n;j++) 
    {
        int win=0,lose=0; 
        for(int k=1;k<=min(i,j);k++)
			if(SG[i-k][j-k])win++; 
			else lose++;
        for(int k=0;k<i;k++)
			if(SG[k][j])win++;
			else lose++;
        for(int k=0;k<j;k++)
			if(SG[i][k])win++;
			else lose++;
			//所有后继都是必胜态,当前是必败态
        if(lose==0)printf("%d %d\n",i,j);
        else SG[i][j]=1; 
    }
}
图论
强联通缩点
#include <bits/stdc++.h>
//luogu P2341
using namespace std;

/*
	scc表示某标号的强连通分量中的点,co表示某个点属于哪个强连通分量
	gao函数是重建图,按照题意寻找有没有链
*/

const int N=10005;

int n,m,x[N*5],y[N*5];
vector<int> e[N],scc[N];
int co[N],color=0;
stack<int> s;
bool vis[N];
int dfn[N],low[N],timer=0;

void tarjan(int u)
{
	dfn[u]=low[u]=++timer;
	s.push(u);
	vis[u]=1;
	for(auto v:e[u])
	{
		if(!dfn[v])
		{
			tarjan(v);
			low[u]=min(low[u],low[v]);
		}
		else if(vis[v]) low[u]=min(low[u],dfn[v]);
	}
	if(low[u]==dfn[u])
	{
		++color;
		int t;
		do
		{
			t=s.top();
			s.pop();
			co[t]=color;
			vis[t]=0;
			scc[color].push_back(t);
		}
		while(u!=t);
	}
}

int f[N];
int _find(int x)
{
	if(x!=f[x]) f[x]=_find(f[x]);
	return f[x];
}
void _merge(int x,int y)
{
	x=_find(x),y=_find(y);
	if(x!=y) f[x]=y;
}

int d[N];
void gao()
{
	for(int i=1;i<=color;++i)
		f[i]=i;
	for(int i=1;i<=m;++i)
	{
		if(co[x[i]]!=co[y[i]])
			_merge(co[x[i]],co[y[i]]),
			d[co[x[i]]]++;
	}
	int F=_find(1);
	for(int i=1;i<=color;++i)
		if(_find(i)!=F) {puts("0");return;}
	int ans=0,tmp=0;
	for(int i=1;i<=color;++i)
	{
		if(d[i]==0)
			ans+=scc[i].size(),tmp++;
	}
	if(tmp>1) ans=0;
	printf("%d",ans);
}

int main()
{
	scanf("%d%d",&n,&m);
	for(int i=1;i<=m;++i)
	{
		scanf("%d%d",&x[i],&y[i]);
		e[x[i]].push_back(y[i]);
	}
	for(int i=1;i<=n;++i)
		if(!dfn[i]) tarjan(i);
	gao();
	return 0;
}
严格次小生成树
#include <bits/stdc++.h>
#define int long long
using namespace std;
const int N=5e5+5;
const int M=1e6+5;
int inf=1e18;
//luogu P1480

/*
	求解严格小于的次小生成树
*/

struct edge
{
	int to,val;
};
vector<edge> e[N];

struct node
{
	int x,y,v;
}eg[M];

void add(int a,int b,int c)
{
	e[a].push_back({b,c});
	e[b].push_back({a,c});
}

int fa[N][35],dep[N],mx[N][35],mn[N][35];
void dfs(int x,int f)
{
	fa[x][0]=f;
	for(auto i:e[x])
	{
		if(i.to==f) continue;
		dep[i.to]=dep[x]+1;
		mx[i.to][0]=i.val;
		mn[i.to][0]=-inf;
		dfs(i.to,x);
	}
}

int n,m;
void cal()
{
	for(int i=1;i<=30;++i)
	{
		for(int j=1;j<=n;++j)
		{
			fa[j][i]=fa[fa[j][i-1]][i-1];
			mx[j][i]=max(mx[j][i-1],mx[fa[j][i-1]][i-1]);
			mn[j][i]=max(mn[j][i-1],mn[fa[j][i-1]][i-1]);
			if(mx[j][i-1]>mx[fa[j][i-1]][i-1])
				mn[j][i]=max(mn[j][i],mx[fa[j][i-1]][i-1]);
			else if(mx[j][i-1]<mx[fa[j][i-1]][i-1])
				mn[j][i]=max(mn[j][i],mx[j][i-1]);
		}
	}
}

int lca(int x,int y)
{
	if(dep[x]<dep[y])
		swap(x,y);
	for(int i=30;i>=0;--i)
		if(dep[fa[x][i]]>=dep[y])
			x=fa[x][i];
	if(x==y) return x;
	for(int i=30;i>=0;--i)
		if(fa[x][i]^fa[y][i])
			x=fa[x][i],
			y=fa[y][i];
	return fa[x][0];
}

int qmax(int u,int v,int val)
{
	int ans=-inf;
	for(int i=30;i>=0;--i)
	{
		if(dep[fa[u][i]]>=dep[v])
		{
			if(val!=mx[u][i])
				ans=max(ans,mx[u][i]);
			else
				ans=max(ans,mn[u][i]);
			u=fa[u][i];
		}
	}
	return ans;
}

bool cmp(node a,node b)
{
	return a.v<b.v;
}

int F[N];
int _find(int x)
{
	if(x!=F[x]) F[x]=_find(F[x]);
	return F[x];
}

bool flg[M];
signed main()
{
	scanf("%lld%lld",&n,&m);
	for(int i=1;i<=m;++i)
	{
		scanf("%lld%lld%lld",&eg[i].x,&eg[i].y,&eg[i].v);
	}
	sort(eg+1,eg+1+m,cmp);
	for(int i=1;i<=n;++i)
	{
		F[i]=i;
	}
	int cnt=0;
	for(int i=1;i<=m;++i)
	{
		int fx=_find(eg[i].x);
		int fy=_find(eg[i].y);
		if(fx!=fy)
		{
			cnt+=eg[i].v;
			F[fx]=fy;
			add(eg[i].x,eg[i].y,eg[i].v);
			flg[i]=1;
		}
	}
	mn[1][0]=-inf;
	dep[1]=1;
	dfs(1,-1);
	cal();
	int ans=inf;
	for(int i=1;i<=m;++i)
	{
		if(flg[i]) continue;
		int x=eg[i].x;
		int y=eg[i].y;
		int v=eg[i].v;
		int l=lca(x,y);
		int mxu=qmax(x,l,v);
		int mxv=qmax(y,l,v);
		ans=min(ans,cnt-max(mxu,mxv)+v);
	}
	printf("%lld",ans);
	return 0;
}
最大流

用链式前向星存的,下面的费用流用vector

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
inline LL read()
{
    LL kk=0,f=1;
    char cc=getchar();
    while(cc<'0'||cc>'9'){if(cc=='-')f=-1;cc=getchar();}
    while(cc>='0'&&cc<='9'){kk=(kk<<1)+(kk<<3)+cc-'0';cc=getchar();}
    return kk*f;
}
void outLL(LL x)
{
	if(x<0){x=~x+1;putchar('-');}
	if(x>9)outLL(x/10);
	putchar(char(x%10+'0'));
}const int maxn=2e5+111;
const LL INF=1e18+7;
LL mod=1e9+7;
const double pi=acos(-1);
#define pii pair<int,int>
#define mp(x,y) make_pair(x,y)
LL head[maxn];
LL n,m;
LL dis[maxn];
queue<LL>Q;
int S,T;
struct Edge
{
	LL to,nex;
	LL flow;
}edge[maxn];
LL cnt=0;
void inittu()
{
	memset(head,-1,sizeof(head));
	memset(edge,0,sizeof(edge));cnt=0;
}

void addedge(LL fr,LL to,LL val)
{
	edge[cnt].to=to;edge[cnt].flow=val;
	edge[cnt].nex=head[fr];head[fr]=cnt;cnt++;
}
bool vis[maxn];
void bfs()
{
	memset(dis,0,sizeof(dis));
	while(!Q.empty())Q.pop();
	dis[S]=0;vis[S]=1;
	Q.push(S);
	while(!Q.empty())
	{
		LL now=Q.front();Q.pop();
		for(int i=head[now];i!=-1;i=edge[i].nex)
		{
			int to=edge[i].to;
			if(!vis[to]&&edge[i].flow)
			{
				Q.push(to);dis[to]=dis[now]+1;
				vis[to]=1;
			}
		}
	}
}
LL dfs(LL now,LL freeflow)
{
	if(now==T)return freeflow;
	LL flow=0,tmp=0;
	for(int i=head[now];i!=-1;i=edge[i].nex)
	{
		LL to=edge[i].to;
		if(dis[to]==(dis[now]+1)&&edge[i].flow)
		{
			tmp=dfs(to,min(edge[i].flow,freeflow));
			if(tmp<=0)continue;
			freeflow-=tmp;
			flow+=tmp;
			edge[i].flow-=tmp;
			edge[i^1].flow+=tmp;
			if(freeflow<=0)break;
		}
	}
	if(flow==0)dis[now]=0;//剪枝
	return flow;
}
LL maxflow()
{
	LL ans=0;
	while(1)
	{
		memset(vis,0,sizeof(vis));
		bfs();
		if(!vis[T])return ans;
		ans+=dfs(S,1e18);
	}
	return ans;
}
int main()
{
	inittu();
	n=read();m=read();S=read();T=read();
	for(int i=1;i<=m;++i)
	{
		LL a,b;LL c;a=read();b=read();c=read();
		addedge(a,b,c);
		addedge(b,a,0ll);//反向弧
	}
	outLL(maxflow());
}
费用流

最小费用流,求最大只需要将费用的边权取负数即可,答案为-mincost

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int maxn=1e7+111;
const LL INF=1e18+7;
LL mod=1e9+7;
const double PI=acos(-1);
#define pii pair<int,int>
#define mp(x,y) make_pair(x,y)
//MCMF
const int N=5e3+15,M=2e6+15;
struct Edge
{
	LL to,flow,cost,rev;//rev是反向边的编号,flow为流上限,cost为花费
};
vector<Edge>edge[N];
inline void addflow(LL x,LL y,LL z,LL w) //加边和反向边,记录对应的反向边的
{
	edge[x].push_back((Edge){y,z,w,edge[y].size()});
	edge[y].push_back((Edge){x,0,-w,edge[x].size()-1});
}
int n,m,s,t;
LL dis[N],maflow,mincost;
bool vis[N];
void init()
{
	for(int i=0;i<=n;++i)edge[i].clear();
	maflow=0;mincost=0;//最小费用最大流
}
bool spfa()
{
    for(int i=0;i<=n;i++)vis[i]=0,dis[i]=INF;
    dis[s]=0;
    queue<int>Q;
    Q.push(s);
    while(!Q.empty())
    {
        int now=Q.front();Q.pop();
        vis[now]=0;
        for(auto e:edge[now])
        if(e.flow&&dis[now]+e.cost<dis[e.to])
        {
            dis[e.to]=dis[now]+e.cost;
            if(!vis[e.to])
            {
            	vis[e.to]=1;Q.push(e.to);
			}
        }
    }
    return dis[t]!=INF;
}
int cur[N];
LL dfs(int now,LL flow)//反向通过増广路更新花费和容量
{
    if(now==t)return flow;
    LL rest=flow,k;
    vis[now]=1;
    for(int i=cur[now];i<edge[now].size();i++)//cur加速
    {
        Edge &e=edge[now][i];
        if(e.flow&&dis[now]+e.cost==dis[e.to]&&!vis[e.to])
        {
            cur[now]=i;
            k=dfs(e.to,min(rest,e.flow));
            e.flow-=k;//容量
            edge[e.to][e.rev].flow+=k;
            mincost+=k*e.cost;//花费
            rest-=k;
        }
    }
    vis[now]=0;
    return flow-rest;
}
void MCMF()
{
    LL delta;
    while(spfa())
    {
        for(int i=0;i<=n;i++)cur[i]=vis[i]=0;
        while(delta=dfs(s,INF))maflow+=delta;
    }
}
//MCMF
int main()
{
	LL u,v,w,f;
	scanf("%d%d%d%d",&n,&m,&s,&t);
	init();
	while(m--)
	{
		scanf("%lld%lld%lld%lld",&u,&v,&w,&f);
		addflow(u,v,w,f);
	}
	MCMF();
	printf("%lld %lld",maflow,mincost);
	return 0;
}

二分图最大匹配(匈牙利)
bool dfs(int x)
{
	for(int i:edge[x])
	{
		if(vis[i])continue;
		vis[i]=1;
		if(!match[i]||dfs(match[i]))
		{
			match[i]=x;return 1;
		}
	}
	return 0;
}
int xiongyali()
{
	int ans=0;
	memset(match,0,sizeof(match));
	for(int i=1;i<=n;++i)
	{
		memset(vis,0,sizeof(vis));
		if(dfs(i))ans++;
	}
	return ans;
}
LCA(倍增)
LL fun[maxn][22],de[maxn];//LL dis[maxn];
void dfsLCA(LL now,LL fa)
{
	fun[now][0]=fa;de[now]=de[fa]+1;
	for(LL i=1;i<=19;++i)
	fun[now][i]=fun[fun[now][i-1]][i-1];
	for(LL i=0;i<edge[now].size();++i)
	{
		LL to=edge[now][i];
		if(to==fa)continue;
		//dis[to]=dis[now]+val;
		dfsLCA(to,now);
	}
}
LL lca(LL x,LL y)
{
	if(de[x]>de[y])swap(x,y);
	for(LL i=19;i>=0;--i)
	if(de[fun[y][i]]>=de[x])y=fun[y][i];
	if(x==y)return x;
	for(LL i=19;i>=0;--i)
	if(fun[x][i]!=fun[y][i])
	{
		x=fun[x][i];y=fun[y][i];
	}
	return fun[x][0];
}
负环SPFA
bool SPFA()
{
	queue<int>Q;
	Q.push(0);vis[0]=1;dis[0]=0;
	while(!Q.empty())
	{
		int now=Q.front();Q.pop();
		vis[now]=0;
		if(++cnt[now]>n||dis[now]<tot)//入队次数和最小负距离作为判断依据
		return 1;
		for(auto to:edge[now])
		{
			if(dis[to.second]>dis[now]+to.first)
			{
				dis[to.second]=dis[now]+to.first;
				ans=min(ans,dis[to.second]);
				if(!vis[to.second])Q.push(to.second);
				vis[to.second]=1;
			}
		}
	}
	return 0;
}
bool SPFA(int sta)//dfs版本,部分情况性能更优
{
	if(dis[sta]<tot)return 1;
	vis[sta]=1;
	for(auto to:edge[sta])
	{
		if(dis[to.second]>dis[sta]+to.first)
		{
			dis[to.second]=dis[sta]+to.first;
			ans=min(ans,dis[to.second]);
			if(vis[to.second])return 1;
			if(SPFA(to.second))return 1;
		}
	}
	vis[sta]=0;
	return 0;
}
字符串
KMP

s1为文本串,s2为查询。

next 数组各值的含义:代表当前字符之前的字符串中,有多大长度的前缀和后缀相同。例如如果next [j] = k,代表j 之前的字符串中有最大长度为k 的前缀和后缀相同。len2-(mnext[len2-1]+1);为最小循环节长度。

char s1[1000222],s2[1000222];int len1,len2;int co=0;
int mnext[1000222];
void cnext()
{
	mnext[0]=-1;int j=0;
	for(int i=1;i<len2;++i)
	{
		j=mnext[i-1];while(j>-1&&s2[j+1]!=s2[i])j=mnext[j];//前结点的后继失配则前溯
		if(s2[j+1]==s2[i])mnext[i]=j+1;
		else mnext[i]=-1;//失配
	}
}
void KMP()
{
	int i=0,j=-1;
	cnext();
	for(;i<len1;++i)
	{
		while(j>-1&&s2[j+1]!=s1[i])j=mnext[j];
		if(s2[j+1]==s1[i])j++;//后继匹配成功
		if(j+1==len2)
		{
			co++;j=mnext[j];
		}
	}
}
int main()
{
	cin>>s1>>s2;
	len1=strlen(s1);len2=strlen(s2);
	KMP();
	cout<<co;
}
Trie

求前缀,等各种操作

void insert(char *str)
{
	int len=strlen(str),p=1;
	for(int i=0;i<len;++i)
	{
		int ch=str[i]-'a';
		if(!trie[p][ch])trie[p][ch]=++tot;
		p=trie[p][ch];
	}
	cnt[p]++;//该处结尾的串
}
int search(char* str)//查询插入的串里面前缀数量
{
	int len=strlen(str),p=1,asd=0;
	for(int i=0;i<len;++i)
	{
		p=trie[p][str[i]-'a'];
		asd+=cnt[p];//加入已经是结尾的串
		if(p==0)break;
	}
	return asd;
}
char ss[maxn],s1[maxn];
int main()
{
	int n,m;n=read();m=read();
	for(int i=1;i<=n;++i)
	{
		scanf("%s",ss);
		insert(ss);
	}
	for(int i=1;i<=m;++i)
	{
		scanf("%s",s1);
		printf("%d\n",search(s1));
	}
}
AC自动机

给定n个模式串和1个文本串,求有多少个模式串在文本串里出现过。查询文本串时,当遇到失配情况时走向fail指向的trie树上的结点继续查询。

struct TRIE
{
	int son[33],fail,cnt;//fail 指针的含义:((最长的(当前字符串的后缀))在Trie上可以查找到)的末尾编号。
	void clear()
	{
		for(int i=0;i<26;++i)son[i]=0;
		fail=0;cnt=0;
	}
}trie[maxn];
int pos=1;
void insert_trie(char *str)//将查询加入trie
{
	int now=1,len=strlen(str);
	for(int i=0;i<len;++i)
	{
		int ch=str[i]-'a';
		if(!trie[now].son[ch])trie[now].son[ch]=++pos;
		now=trie[now].son[ch];
	}
	trie[now].cnt++;//记录这个串出现次数
    //trie[now].cnt=id;
}
void getfail()//取得失配情况
{
	for(int i=0;i<26;++i)trie[0].son[i]=1;//完全失配,回到根
	queue<int>Q;
	Q.push(1);trie[1].fail=0;//0为根的虚拟父亲
	while(!Q.empty())
	{
		int now=Q.front();Q.pop();
		for(int i=0;i<26;++i)
		{
			int to=trie[now].son[i];
			int fafail=trie[now].fail;
			if(!to)//不存在下一个这个节点,让下一个这个结点为父节点失配情况的下一个该节点,统一合法性
			{
				trie[now].son[i]=trie[fafail].son[i];
				continue;
			}
			trie[to].fail=trie[fafail].son[i];//当前失配的情况为父亲失配的相应的下一个节点
			Q.push(to);
		}
	}
}
int query_trie(char *str)
{
	int len=strlen(str),fr,to,asd=0;
	fr=1;
	for(int i=0;i<len;++i)//在trie上遍历文本串
	{
		int ch=str[i]-'a';
		to=trie[fr].son[ch];
        /*
        while(to>1)
		{
			if(trie[to].cnt)vis[trie[to].cnt]++;//统计每种串的出现次数
			to=trie[to].fail;
		}
        */
		while(to>1&&trie[to].cnt!=-1)//除非完全失配或者之前走过了,不然继续
		{
			asd+=trie[to].cnt;
			trie[to].cnt=-1;//本串的所有情况都加过了
			to=trie[to].fail;
		}
		fr=trie[fr].son[ch];
	}
	return asd;
}
char ss[maxn];
int main()
{
	int T=read();
	while(T--)
	{
		scanf("%s",ss);
		insert_trie(ss);
	}
	getfail();
	scanf("%s",ss);
	printf("%d\n",query_trie(ss));
}
其他
莫队
#define pii pair<int,int>
#define mp(x,y) make_pair(x,y)
#define sf(a) scanf("%lld",&a)
int n,m;
int a[maxn];
int belong[maxn];
struct zj
{
	LL l,r,k,id;
	bool operator<(const zj b)const
	{
		if(belong[l]==belong[b.l])return r<b.r;
		return belong[l]<belong[b.l];
	}
}que[maxn];
LL cnt[maxn];
set<LL>S;
vector<LL>V;
bool cmp(zj a,zj b)
{
	return a.id<b.id;
}
void del(int x)
{
	cnt[a[x]]--;
}
void add(int x)
{
	cnt[a[x]]++;
}
LL ans[maxn];
int main()
{
	n=read();m=read();
	for(int i=1;i<=n;++i)a[i]=read(),S.insert(a[i]);
	for(auto it=S.begin();it!=S.end();++it)V.push_back(*it);
	int sn=sqrt(n),len=ceil((double)n/sn);
	for(int i=1;i<=sn;++i)
	for(int j=(i-1)*len+1;j<=min(i*len,n);++j)belong[j]=i;
	
	for(int i=1;i<=m;++i)
	{
		que[i].l=read();que[i].r=read();que[i].k=read();que[i].id=i;
	}
	sort(que+1,que+1+m);
	int l,r;
	l=r=que[1].l;cnt[a[l]]++;S.insert(a[l]);
	for(int i=1;i<=m;++i)
	{
		while(l<que[i].l){del(l++);}
		while(l>que[i].l){add(--l);}
		while(r<que[i].r){add(++r);}
		while(r>que[i].r){del(r--);}
		int co=0;
		for(auto kk:V)
		{
			if(cnt[kk]>0&&__gcd(cnt[kk],que[i].k)==1)
			co++;
		}
		ans[que[i].id]=co;
	}
	for(int i=1;i<=m;++i)
	{
		printf("%lld\n",ans[i]);
	}
}
获取一个排列的序号,从0开始
int GetIndex(int* a)
{
    int res=0;
	for(int i=1;i<8;i++)
	{
		int cnt=0;
		for(int j=i+1;j<=8;j++)
		if(a[i]>a[j])
		cnt++;
		res+=cnt*jie[8-i];
	}
	return res;
}
几何板子
#include<bits/stdc++.h>
using namespace std;

#define db double
const db EPS=1e-9;
inline int sign(db a){return a<-EPS?-1:a>EPS;}
inline int cmp(db a,db b){return sign(a-b);}
struct P
{
    db x,y;
    P(){}
    P(db x,db y):x(x),y(y){}
    P operator+(P p){return {x+p.x,y+p.y};}
    P operator-(P p){return {x-p.x,y-p.y};}
    P operator*(db d){return {x*d,y*d};}
    P operator/(db d){return {x/d,y/d};}
    bool operator<(P p) const
    {
        int c=cmp(x,p.x);
        if(c) return c==-1;
        return cmp(y,p.y)==-1;
    }
    bool operator==(P o) const
    {
        return cmp(x,o.x)==0&&cmp(y,o.y)==0;
    }
    db distTo(P p){return (*this-p).abs();}//点距离
    db alpha(){return atan2(y,x);}
    void read(){scanf("%lf%lf",&x,&y);}
    void write(){printf("(%.10f,%.10f)\n",x,y);}
    db abs(){return sqrt(abs2());}//模长
    db abs2(){return x*x+y*y;}
    P rot90(){return P(-y,x);}//转90度
    P unit(){return *this/abs();}//单位向量
    int quad() const {return sign(y)==1||(sign(y)==0&&sign(x)>=0);}
    db dot(P p){return x*p.x+y*p.y;}//点积
    db det(P p){return x*p.y-y*p.x;}//叉积
    P rot(db an){return {x*cos(an)-y*sin(an),x*sin(an)+y*cos(an)};}//旋转角度
};

int compareAngle(P a,P b)
{
    if(a.quad()!=b.quad()) return a.quad()<b.quad();
    return sign(a.det(b))>0;
}

//For segment
#define cross(p1,p2,p3) ((p2.x-p1.x)*(p3.y-p1.y)-(p3.x-p1.x)*(p2.y-p1.y))
#define crossOp(p1,p2,p3) sign(cross(p1,p2,p3))

bool chkLL(P p1,P p2,P q1,P q2) //0:parallel
{
    db a1=cross(q1,q2,p1),a2=-cross(q1,q2,p2);
    return sign(a1+a2)!=0;
}

P isLL(P p1,P p2,P q1,P q2) //crossover point if chkLL()
{
    db a1=cross(q1,q2,p1),a2=-cross(q1,q2,p2);
    return (p1*a2+p2*a1)/(a1+a2);
}

bool intersect(db l1,db r1,db l2,db r2)//快速排斥实验
{
    if(l1>r1) swap(l1,r1);if(l2>r2) swap(l2,r2);
    return !(cmp(r1,l2)==-1||cmp(r2,l1)==-1);
}

bool isSS(P p1,P p2,P q1,P q2)//快速跨立实验
{
    return intersect(p1.x,p2.x,q1.x,q2.x)&&intersect(p1.y,p2.y,q1.y,q2.y)&&
    crossOp(p1,p2,q1)*crossOp(p1,p2,q2)<=0&&crossOp(q1,q2,p1)*crossOp(q1,q2,p2)<=0;
}

bool isSS_strict(P p1,P p2,P q1,P q2)//判断线段是否严格相交(不包括点在线段上的情形),只用跨立实验即可
{
    return crossOp(p1,p2,q1)*crossOp(p1,p2,q2)<0
    &&crossOp(q1,q2,p1)*crossOp(q1,q2,p2)<0;
}

bool isMiddle(db a,db m,db b)
{
    return sign(a-m)==0||sign(b-m)==0||(a<m!=b<m);
}

bool isMiddle(P a,P m,P b)//判断m点是否在a和b形成的矩形区域内
{
    return isMiddle(a.x,m.x,b.x)&&isMiddle(a.y,m.y,b.y);
}

bool onSeg(P p1,P p2,P q)//q is in the seg of p1-p2
{
    return crossOp(p1,p2,q)==0&&isMiddle(p1,q,p2);
}

bool onSeg_strict(P p1,P p2,P q)
{
    return crossOp(p1,p2,q)==0&&sign((q-p1).dot(p1-p2))*sign((q-p2).dot(p1-p2))<0;
}

P proj(P p1,P p2,P q)//点在线段上的射影
{
    P dir=p2-p1;
    return p1+dir*(dir.dot(q-p1)/dir.abs2());
}

P reflect(P p1,P p2,P q)//点关于线段的对称点
{
    return proj(p1,p2,q)*2-q;
}

db nearest(P p1,P p2,P q)//点到线段的距离
{
    P h=proj(p1,p2,q);
    if(isMiddle(p1,h,p2))
        return q.distTo(h);
    return min(p1.distTo(q),p2.distTo(q));
}

db disSS(P p1,P p2,P q1,P q2) //dist of 2 segments
{
    if(isSS(p1,p2,q1,q2)) return 0;
    return min(min(nearest(p1,p2,q1),nearest(p1,p2,q2)),min(nearest(q1,q2,p1),nearest(q1,q2,p2)));
}

db rad(P p1,P p2)//OP1与OP2的夹角
{
    return atan2l(p1.det(p2),p1.dot(p2));
}

db area(vector<P> ps)//多边形面积
{
    db ret=0;
    for(int i=0;i<ps.size();i++)
        ret+=ps[i].det(ps[(i+1)%ps.size()]);
    return ret/2;
}

int contain(vector<P> ps,P p) //2:inside,1:on_seg,0:outside
{
    int n=ps.size(),ret=0;
    for(int i=0;i<n;i++)
    {
        P u=ps[i],v=ps[(i+1)%n];
        if(onSeg(u,v,p)) return 1;
        if(cmp(u.y,v.y)<=0) swap(u,v);
        if(cmp(p.y,u.y)>0||cmp(p.y,v.y)<=0) continue;
        ret^=crossOp(p,u,v)>0;
    }
    return ret*2;
}

int convexContain(vector<P> &ps,P p) //1:inside|on_seg,0:outside
{
    int n=ps.size(),l=1,r=n-1,mid;
    while(l<r)
    {
        mid=(l+r+1)>>1;
        if(sign(cross(ps[0],ps[mid],p))<0) r=mid-1; else l=mid;
    }
    if(l==1&&sign(cross(ps[0],ps[1],p))<0) return 0;
    if(l==n-1&&onSeg(ps[0],ps[n-1],p)) return 1;
    if(l!=n-1&&sign(cross(ps[l],ps[l+1],p))>=0) return 1;
    return 0;
}

vector<P> convexHull(vector<P> ps)
{
    int n=ps.size();if(n<=1) return ps;
    sort(ps.begin(),ps.end());
    vector<P> qs(n*2);int k=0;
    for(int i=0;i<n;qs[k++]=ps[i++])
        while(k>1&&crossOp(qs[k-2],qs[k-1],ps[i])<=0) --k;
    for(int i=n-2,t=k;i>=0;qs[k++]=ps[i--])
        while(k>t&&crossOp(qs[k-2],qs[k-1],ps[i])<=0) --k;
    qs.resize(k-1);
    return qs;
}

db convexDiameter(vector<P> ps)
{
    int n=ps.size();if(n<=1) return 0;
    int is=0,js=0;
    for(int k=1;k<n;k++) is=ps[k]<ps[is]?k:is,js=ps[js]<ps[k]?js:k;
    int i=is,j=js;
    db ret=ps[i].distTo(ps[j]);
    do{
        if((ps[(i+1)%n]-ps[i]).det(ps[(j+1)%n]-ps[j])>=0) (++j)%=n;
        else (++i)%=n;
        ret=max(ret,ps[i].distTo(ps[j]));
    }while(i!=is||j!=js);
    return ret;
}

struct L // p[0]->p[1]
{
    P p[2];
    L(P k1,P k2){p[0]=k1,p[1]=k2;}
    P& operator [] (int k){return p[k];}
    int include(P k){return sign((p[1]-p[0]).det(k-p[0]))>0;}
    P dir(){return p[1]-p[0];}
    L push(db dis) // push dis (left hand)
    {
        P delta=(p[1]-p[0]).rot90().unit()*dis;
        return {p[0]+delta,p[1]+delta};
    }
};

bool parallel(L l0,L l1){return sign(l0.dir().det(l1.dir()))==0;}

bool sameDir(L l0,L l1){return parallel(l0,l1)&&sign(l0.dir().dot(l1.dir()))==1;}

bool operator < (L l0,L l1)
{
    if(sameDir(l0,l1)) return l1.include(l0[0]);
    return compareAngle(l0.dir(),l1.dir());
}

P isLL(L l0,L l1){return isLL(l0[0],l0[1],l1[0],l1[1]);}

bool check(L u,L v,L w){return w.include(isLL(u,v));}

vector<P> halfPlaneIS(vector<L> &l)
{
    sort(l.begin(),l.end());
    deque<L> q;
    for(int i=0;i<(int)l.size();i++)
    {
        if(i&&sameDir(l[i],l[i-1])) continue;
        while(q.size()>1&&!check(q[q.size()-2],q[q.size()-1],l[i])) q.pop_back();
        while(q.size()>1&&!check(q[1],q[0],l[i])) q.pop_front();
        q.push_back(l[i]);
    }
    while(q.size()>2&&!check(q[q.size()-2],q[q.size()-1],q[0])) q.pop_back();
    while(q.size()>2&&!check(q[1],q[0],q[q.size()-1])) q.pop_front();
    vector<P> ret;
    for(int i=0;i<(int)q.size();i++) ret.push_back(isLL(q[i],q[(i+1)%q.size()]));
    return ret;
}

vector<P> v;
int n,l;

int main()
{
    scanf("%d%d",&n,&l);
    for(int i=0,x,y;i<n;i++)
    {
        scanf("%d%d",&x,&y);
        v.emplace_back(x+l,y);
        v.emplace_back(x,y+l);
        v.emplace_back(x,y-l);
        v.emplace_back(x-l,y);
    }
    v=convexHull(v);
    v.push_back(v[0]);
    double ans=0;
    for(int i=0;i<v.size()-1;i++)
        ans+=v[i].distTo(v[i+1]);
    printf("%.10f",ans);
    return 0;
}

偷来的几何
#define mp make_pair
#define fi first
#define se second
#define pb push_back
typedef double db;
const db eps=1e-6;
const db pi=acos(-1);
int sign(db k){
    if (k>eps) return 1; else if (k<-eps) return -1; return 0;
}
int cmp(db k1,db k2){return sign(k1-k2);}
int inmid(db k1,db k2,db k3){return sign(k1-k3)*sign(k2-k3)<=0;}// k3 在 [k1,k2] 内 
struct point{
    db x,y;
    point operator + (const point &k1) const{return (point){k1.x+x,k1.y+y};}
    point operator - (const point &k1) const{return (point){x-k1.x,y-k1.y};}
    point operator * (db k1) const{return (point){x*k1,y*k1};}
    point operator / (db k1) const{return (point){x/k1,y/k1};}
    int operator == (const point &k1) const{return cmp(x,k1.x)==0&&cmp(y,k1.y)==0;}
    // 逆时针旋转 
    point turn(db k1){return (point){x*cos(k1)-y*sin(k1),x*sin(k1)+y*cos(k1)};}
    point turn90(){return (point){-y,x};}
    bool operator < (const point k1) const{
        int a=cmp(x,k1.x);
        if (a==-1) return 1; else if (a==1) return 0; else return cmp(y,k1.y)==-1;
    }
    db abs(){return sqrt(x*x+y*y);}
    db abs2(){return x*x+y*y;}
    db dis(point k1){return ((*this)-k1).abs();}
    point unit(){db w=abs(); return (point){x/w,y/w};}
    void scan(){double k1,k2; scanf("%lf%lf",&k1,&k2); x=k1; y=k2;}
    void print(){printf("%.11lf %.11lf\n",x,y);}
    db getw(){return atan2(y,x);} 
    point getdel(){if (sign(x)==-1||(sign(x)==0&&sign(y)==-1)) return (*this)*(-1); else return (*this);}
	int getP() const{return sign(y)==1||(sign(y)==0&&sign(x)==-1);}
};
int inmid(point k1,point k2,point k3){return inmid(k1.x,k2.x,k3.x)&&inmid(k1.y,k2.y,k3.y);}
db cross(point k1,point k2){return k1.x*k2.y-k1.y*k2.x;}
db dot(point k1,point k2){return k1.x*k2.x+k1.y*k2.y;}
db rad(point k1,point k2){return atan2(cross(k1,k2),dot(k1,k2));}
// -pi -> pi
int compareangle (point k1,point k2){
    return k1.getP()<k2.getP()||(k1.getP()==k2.getP()&&sign(cross(k1,k2))>0);
}
point proj(point k1,point k2,point q){ // q 到直线 k1,k2 的投影 
    point k=k2-k1; return k1+k*(dot(q-k1,k)/k.abs2());
}
point reflect(point k1,point k2,point q){return proj(k1,k2,q)*2-q;}
int clockwise(point k1,point k2,point k3){// k1 k2 k3 逆时针 1 顺时针 -1 否则 0  
    return sign(cross(k2-k1,k3-k1));
}
int checkLL(point k1,point k2,point k3,point k4){// 求直线 (L) 线段 (S)k1,k2 和 k3,k4 的交点 
    return cmp(cross(k3-k1,k4-k1),cross(k3-k2,k4-k2))!=0;
}
point getLL(point k1,point k2,point k3,point k4){
    db w1=cross(k1-k3,k4-k3),w2=cross(k4-k3,k2-k3); return (k1*w2+k2*w1)/(w1+w2);
}
int intersect(db l1,db r1,db l2,db r2){
    if (l1>r1) swap(l1,r1); if (l2>r2) swap(l2,r2); return cmp(r1,l2)!=-1&&cmp(r2,l1)!=-1;
}
int checkSS(point k1,point k2,point k3,point k4){
    return intersect(k1.x,k2.x,k3.x,k4.x)&&intersect(k1.y,k2.y,k3.y,k4.y)&&
    sign(cross(k3-k1,k4-k1))*sign(cross(k3-k2,k4-k2))<=0&&
    sign(cross(k1-k3,k2-k3))*sign(cross(k1-k4,k2-k4))<=0;
}
db disSP(point k1,point k2,point q){
    point k3=proj(k1,k2,q);
    if (inmid(k1,k2,k3)) return q.dis(k3); else return min(q.dis(k1),q.dis(k2));
}
db disSS(point k1,point k2,point k3,point k4){
    if (checkSS(k1,k2,k3,k4)) return 0;
    else return min(min(disSP(k1,k2,k3),disSP(k1,k2,k4)),min(disSP(k3,k4,k1),disSP(k3,k4,k2)));
}
int onS(point k1,point k2,point q){return inmid(k1,k2,q)&&sign(cross(k1-q,k2-k1))==0;}
struct circle{
    point o; db r;
    void scan(){o.scan(); scanf("%lf",&r);}
    int inside(point k){return cmp(r,o.dis(k));}
};
struct line{
    // p[0]->p[1]
    point p[2];
    line(point k1,point k2){p[0]=k1; p[1]=k2;}
    point& operator [] (int k){return p[k];}
    int include(point k){return sign(cross(p[1]-p[0],k-p[0]))>0;}
    point dir(){return p[1]-p[0];}
    line push(){ // 向外 ( 左手边 ) 平移 eps 
        const db eps = 1e-6;
        point delta=(p[1]-p[0]).turn90().unit()*eps;
        return {p[0]-delta,p[1]-delta};
    }
};
point getLL(line k1,line k2){return getLL(k1[0],k1[1],k2[0],k2[1]);}
int parallel(line k1,line k2){return sign(cross(k1.dir(),k2.dir()))==0;}
int sameDir(line k1,line k2){return parallel(k1,k2)&&sign(dot(k1.dir(),k2.dir()))==1;}
int operator < (line k1,line k2){
    if (sameDir(k1,k2)) return k2.include(k1[0]); 
    return compareangle(k1.dir(),k2.dir());
}
int checkpos(line k1,line k2,line k3){return k3.include(getLL(k1,k2));}
vector<line> getHL(vector<line> &L){ // 求半平面交 , 半平面是逆时针方向 , 输出按照逆时针
    sort(L.begin(),L.end()); deque<line> q;
    for (int i=0;i<(int)L.size();i++){
        if (i&&sameDir(L[i],L[i-1])) continue;
        while (q.size()>1&&!checkpos(q[q.size()-2],q[q.size()-1],L[i])) q.pop_back();
        while (q.size()>1&&!checkpos(q[1],q[0],L[i])) q.pop_front();
        q.push_back(L[i]);
    }
    while (q.size()>2&&!checkpos(q[q.size()-2],q[q.size()-1],q[0])) q.pop_back();
    while (q.size()>2&&!checkpos(q[1],q[0],q[q.size()-1])) q.pop_front();
    vector<line>ans; for (int i=0;i<q.size();i++) ans.push_back(q[i]);
    return ans;
}
db closepoint(vector<point>&A,int l,int r){ // 最近点对 , 先要按照 x 坐标排序 
    if (r-l<=5){
        db ans=1e20;
        for (int i=l;i<=r;i++) for (int j=i+1;j<=r;j++) ans=min(ans,A[i].dis(A[j]));
        return ans;
    }
    int mid=l+r>>1; db ans=min(closepoint(A,l,mid),closepoint(A,mid+1,r));
    vector<point>B; for (int i=l;i<=r;i++) if (abs(A[i].x-A[mid].x)<=ans) B.push_back(A[i]);
    sort(B.begin(),B.end(),[](point k1,point k2){return k1.y<k2.y;});
    for (int i=0;i<B.size();i++) for (int j=i+1;j<B.size()&&B[j].y-B[i].y<ans;j++) ans=min(ans,B[i].dis(B[j]));
    return ans;
}
int checkposCC(circle k1,circle k2){// 返回两个圆的公切线数量
    if (cmp(k1.r,k2.r)==-1) swap(k1,k2);
    db dis=k1.o.dis(k2.o);  int w1=cmp(dis,k1.r+k2.r),w2=cmp(dis,k1.r-k2.r);
    if (w1>0) return 4; else if (w1==0) return 3; else if (w2>0) return 2; 
    else if (w2==0) return 1; else return 0;
}
vector<point> getCL(circle k1,point k2,point k3){ // 沿着 k2->k3 方向给出 , 相切给出两个 
    point k=proj(k2,k3,k1.o); db d=k1.r*k1.r-(k-k1.o).abs2();
    if (sign(d)==-1) return {};
    point del=(k3-k2).unit()*sqrt(max((db)0.0,d)); return {k-del,k+del};
}
vector<point> getCC(circle k1,circle k2){// 沿圆 k1 逆时针给出 , 相切给出两个 
    int pd=checkposCC(k1,k2); if (pd==0||pd==4) return {};
    db a=(k2.o-k1.o).abs2(),cosA=(k1.r*k1.r+a-k2.r*k2.r)/(2*k1.r*sqrt(max(a,(db)0.0)));
    db b=k1.r*cosA,c=sqrt(max((db)0.0,k1.r*k1.r-b*b));
    point k=(k2.o-k1.o).unit(),m=k1.o+k*b,del=k.turn90()*c;
    return {m-del,m+del};
} 
vector<point> TangentCP(circle k1,point k2){// 沿圆 k1 逆时针给出 
    db a=(k2-k1.o).abs(),b=k1.r*k1.r/a,c=sqrt(max((db)0.0,k1.r*k1.r-b*b));
    point k=(k2-k1.o).unit(),m=k1.o+k*b,del=k.turn90()*c;
    return {m-del,m+del};
} 
vector<line> TangentoutCC(circle k1,circle k2){
    int pd=checkposCC(k1,k2); if (pd==0) return {}; 
    if (pd==1){point k=getCC(k1,k2)[0]; return {(line){k,k}};}
    if (cmp(k1.r,k2.r)==0){
        point del=(k2.o-k1.o).unit().turn90().getdel();
        return {(line){k1.o-del*k1.r,k2.o-del*k2.r},(line){k1.o+del*k1.r,k2.o+del*k2.r}};
    } else {
        point p=(k2.o*k1.r-k1.o*k2.r)/(k1.r-k2.r);
        vector<point>A=TangentCP(k1,p),B=TangentCP(k2,p);
        vector<line>ans; for (int i=0;i<A.size();i++) ans.push_back((line){A[i],B[i]}); 
        return ans;
    }
}
vector<line> TangentinCC(circle k1,circle k2){
    int pd=checkposCC(k1,k2); if (pd<=2) return {};
    if (pd==3){point k=getCC(k1,k2)[0]; return {(line){k,k}};} 
    point p=(k2.o*k1.r+k1.o*k2.r)/(k1.r+k2.r);
    vector<point>A=TangentCP(k1,p),B=TangentCP(k2,p);
    vector<line>ans; for (int i=0;i<A.size();i++) ans.push_back((line){A[i],B[i]}); 
    return ans;
}
vector<line> TangentCC(circle k1,circle k2){
    int flag=0; if (k1.r<k2.r) swap(k1,k2),flag=1;
    vector<line>A=TangentoutCC(k1,k2),B=TangentinCC(k1,k2);
    for (line k:B) A.push_back(k); 
    if (flag) for (line &k:A) swap(k[0],k[1]);
    return A;
}
db getarea(circle k1,point k2,point k3){
    // 圆 k1 与三角形 k2 k3 k1.o 的有向面积交
    point k=k1.o; k1.o=k1.o-k; k2=k2-k; k3=k3-k;
    int pd1=k1.inside(k2),pd2=k1.inside(k3); 
    vector<point>A=getCL(k1,k2,k3);
    if (pd1>=0){
        if (pd2>=0) return cross(k2,k3)/2;
        return k1.r*k1.r*rad(A[1],k3)/2+cross(k2,A[1])/2;
    } else if (pd2>=0){ 
        return k1.r*k1.r*rad(k2,A[0])/2+cross(A[0],k3)/2;
    }else {
        int pd=cmp(k1.r,disSP(k2,k3,k1.o));
        if (pd<=0) return k1.r*k1.r*rad(k2,k3)/2;
        return cross(A[0],A[1])/2+k1.r*k1.r*(rad(k2,A[0])+rad(A[1],k3))/2;
    }
}
circle getcircle(point k1,point k2,point k3){
    db a1=k2.x-k1.x,b1=k2.y-k1.y,c1=(a1*a1+b1*b1)/2;
    db a2=k3.x-k1.x,b2=k3.y-k1.y,c2=(a2*a2+b2*b2)/2;
    db d=a1*b2-a2*b1;
    point o=(point){k1.x+(c1*b2-c2*b1)/d,k1.y+(a1*c2-a2*c1)/d};
    return (circle){o,k1.dis(o)};
}
circle getScircle(vector<point> A){
    random_shuffle(A.begin(),A.end());
    circle ans=(circle){A[0],0};
    for (int i=1;i<A.size();i++)
        if (ans.inside(A[i])==-1){
            ans=(circle){A[i],0};
            for (int j=0;j<i;j++)
                if (ans.inside(A[j])==-1){
                    ans.o=(A[i]+A[j])/2; ans.r=ans.o.dis(A[i]);
                    for (int k=0;k<j;k++)
                        if (ans.inside(A[k])==-1)
                            ans=getcircle(A[i],A[j],A[k]);
                }
        }
    return ans;
}
db area(vector<point> A){ // 多边形用 vector<point> 表示 , 逆时针 
    db ans=0;
    for (int i=0;i<A.size();i++) ans+=cross(A[i],A[(i+1)%A.size()]);
    return ans/2;
}
int checkconvex(vector<point>A){
    int n=A.size(); A.push_back(A[0]); A.push_back(A[1]);
    for (int i=0;i<n;i++) if (sign(cross(A[i+1]-A[i],A[i+2]-A[i]))==-1) return 0;
    return 1;
}
int contain(vector<point>A,point q){ // 2 内部 1 边界 0 外部
    int pd=0; A.push_back(A[0]);
    for (int i=1;i<A.size();i++){
        point u=A[i-1],v=A[i];
        if (onS(u,v,q)) return 1; if (cmp(u.y,v.y)>0) swap(u,v);
        if (cmp(u.y,q.y)>=0||cmp(v.y,q.y)<0) continue;
        if (sign(cross(u-v,q-v))<0) pd^=1;
    }
    return pd<<1;
}
vector<point> ConvexHull(vector<point>A,int flag=1){ // flag=0 不严格 flag=1 严格 
    int n=A.size(); vector<point>ans(n*2); 
    sort(A.begin(),A.end()); int now=-1;
    for (int i=0;i<A.size();i++){
        while (now>0&&sign(cross(ans[now]-ans[now-1],A[i]-ans[now-1]))<flag) now--;
        ans[++now]=A[i];
    } int pre=now;
    for (int i=n-2;i>=0;i--){
        while (now>pre&&sign(cross(ans[now]-ans[now-1],A[i]-ans[now-1]))<flag) now--;
        ans[++now]=A[i];
    } ans.resize(now); return ans;
}
db convexDiameter(vector<point>A){
    int now=0,n=A.size(); db ans=0;
    for (int i=0;i<A.size();i++){
        now=max(now,i);
        while (1){
            db k1=A[i].dis(A[now%n]),k2=A[i].dis(A[(now+1)%n]);
            ans=max(ans,max(k1,k2)); if (k2>k1) now++; else break;
        }
    }
    return ans;
}
vector<point> convexcut(vector<point>A,point k1,point k2){
    // 保留 k1,k2,p 逆时针的所有点
    int n=A.size(); A.push_back(A[0]); vector<point>ans;
    for (int i=0;i<n;i++){
        int w1=clockwise(k1,k2,A[i]),w2=clockwise(k1,k2,A[i+1]);
        if (w1>=0) ans.push_back(A[i]);
        if (w1*w2<0) ans.push_back(getLL(k1,k2,A[i],A[i+1]));
    }
    return ans;
}
int checkPoS(vector<point>A,point k1,point k2){
    // 多边形 A 和直线 ( 线段 )k1->k2 严格相交 , 注释部分为线段
    struct ins{
        point m,u,v;
        int operator < (const ins& k) const {return m<k.m;}
    }; vector<ins>B;
    //if (contain(A,k1)==2||contain(A,k2)==2) return 1;
    vector<point>poly=A; A.push_back(A[0]); 
    for (int i=1;i<A.size();i++) if (checkLL(A[i-1],A[i],k1,k2)){
        point m=getLL(A[i-1],A[i],k1,k2); 
        if (inmid(A[i-1],A[i],m)/*&&inmid(k1,k2,m)*/) B.push_back((ins){m,A[i-1],A[i]});
    }
    if (B.size()==0) return 0; sort(B.begin(),B.end()); 
    int now=1; while (now<B.size()&&B[now].m==B[0].m) now++; 
    if (now==B.size()) return 0;
    int flag=contain(poly,(B[0].m+B[now].m)/2);
    if (flag==2) return 1;
    point d=B[now].m-B[0].m;
    for (int i=now;i<B.size();i++){
        if (!(B[i].m==B[i-1].m)&&flag==2) return 1;
        int tag=sign(cross(B[i].v-B[i].u,B[i].m+d-B[i].u));
        if (B[i].m==B[i].u||B[i].m==B[i].v) flag+=tag; else flag+=tag*2;
    }
    //return 0;
    return flag==2;
}
int checkinp(point r,point l,point m){
	if (compareangle(l,r)){return compareangle(l,m)&&compareangle(m,r);}
	return compareangle(l,m)||compareangle(m,r);
}
int checkPosFast(vector<point>A,point k1,point k2){ // 快速检查线段是否和多边形严格相交
	if (contain(A,k1)==2||contain(A,k2)==2) return 1; if (k1==k2) return 0;
	A.push_back(A[0]); A.push_back(A[1]);
	for (int i=1;i+1<A.size();i++)
		if (checkLL(A[i-1],A[i],k1,k2)){
			point now=getLL(A[i-1],A[i],k1,k2);
			if (inmid(A[i-1],A[i],now)==0||inmid(k1,k2,now)==0) continue;
			if (now==A[i]){
				if (A[i]==k2) continue;
				point pre=A[i-1],ne=A[i+1];
				if (checkinp(pre-now,ne-now,k2-now)) return 1;
			} else if (now==k1){
				if (k1==A[i-1]||k1==A[i]) continue;
				if (checkinp(A[i-1]-k1,A[i]-k1,k2-k1)) return 1;
			} else if (now==k2||now==A[i-1]) continue;
			else return 1;
		}
	return 0;
}
// 拆分凸包成上下凸壳 凸包尽量都随机旋转一个角度来避免出现相同横坐标 
// 尽量特判只有一个点的情况 凸包逆时针
void getUDP(vector<point>A,vector<point>&U,vector<point>&D){
    db l=1e100,r=-1e100;
    for (int i=0;i<A.size();i++) l=min(l,A[i].x),r=max(r,A[i].x);
    int wherel,wherer;
    for (int i=0;i<A.size();i++) if (cmp(A[i].x,l)==0) wherel=i;
    for (int i=A.size();i;i--) if (cmp(A[i-1].x,r)==0) wherer=i-1;
    U.clear(); D.clear(); int now=wherel;
    while (1){D.push_back(A[now]); if (now==wherer) break; now++; if (now>=A.size()) now=0;}
    now=wherel;
    while (1){U.push_back(A[now]); if (now==wherer) break; now--; if (now<0) now=A.size()-1;}
}
// 需要保证凸包点数大于等于 3,2 内部 ,1 边界 ,0 外部
int containCoP(const vector<point>&U,const vector<point>&D,point k){
    db lx=U[0].x,rx=U[U.size()-1].x;
    if (k==U[0]||k==U[U.size()-1]) return 1;
    if (cmp(k.x,lx)==-1||cmp(k.x,rx)==1) return 0;
    int where1=lower_bound(U.begin(),U.end(),(point){k.x,-1e100})-U.begin();
    int where2=lower_bound(D.begin(),D.end(),(point){k.x,-1e100})-D.begin();
    int w1=clockwise(U[where1-1],U[where1],k),w2=clockwise(D[where2-1],D[where2],k);
    if (w1==1||w2==-1) return 0; else if (w1==0||w2==0) return 1; return 2;
}
// d 是方向 , 输出上方切点和下方切点
pair<point,point> getTangentCow(const vector<point> &U,const vector<point> &D,point d){
    if (sign(d.x)<0||(sign(d.x)==0&&sign(d.y)<0)) d=d*(-1);
    point whereU,whereD;
    if (sign(d.x)==0) return mp(U[0],U[U.size()-1]);
    int l=0,r=U.size()-1,ans=0;
    while (l<r){int mid=l+r>>1; if (sign(cross(U[mid+1]-U[mid],d))<=0) l=mid+1,ans=mid+1; else r=mid;}
    whereU=U[ans]; l=0,r=D.size()-1,ans=0;
    while (l<r){int mid=l+r>>1; if (sign(cross(D[mid+1]-D[mid],d))>=0) l=mid+1,ans=mid+1; else r=mid;}
    whereD=D[ans]; return mp(whereU,whereD);
}
// 先检查 contain, 逆时针给出
pair<point,point> getTangentCoP(const vector<point>&U,const vector<point>&D,point k){
    db lx=U[0].x,rx=U[U.size()-1].x;
    if (k.x<lx){
        int l=0,r=U.size()-1,ans=U.size()-1;
        while (l<r){int mid=l+r>>1; if (clockwise(k,U[mid],U[mid+1])==1) l=mid+1; else ans=mid,r=mid;}
        point w1=U[ans]; l=0,r=D.size()-1,ans=D.size()-1;
        while (l<r){int mid=l+r>>1; if (clockwise(k,D[mid],D[mid+1])==-1) l=mid+1; else ans=mid,r=mid;}
        point w2=D[ans]; return mp(w1,w2);
    } else if (k.x>rx){
        int l=1,r=U.size(),ans=0;
        while (l<r){int mid=l+r>>1; if (clockwise(k,U[mid],U[mid-1])==-1) r=mid; else ans=mid,l=mid+1;}
        point w1=U[ans]; l=1,r=D.size(),ans=0;
        while (l<r){int mid=l+r>>1; if (clockwise(k,D[mid],D[mid-1])==1) r=mid; else ans=mid,l=mid+1;}
        point w2=D[ans]; return mp(w2,w1);
    } else {
        int where1=lower_bound(U.begin(),U.end(),(point){k.x,-1e100})-U.begin();
        int where2=lower_bound(D.begin(),D.end(),(point){k.x,-1e100})-D.begin();
        if ((k.x==lx&&k.y>U[0].y)||(where1&&clockwise(U[where1-1],U[where1],k)==1)){
            int l=1,r=where1+1,ans=0;
            while (l<r){int mid=l+r>>1; if (clockwise(k,U[mid],U[mid-1])==1) ans=mid,l=mid+1; else r=mid;}
            point w1=U[ans]; l=where1,r=U.size()-1,ans=U.size()-1;
            while (l<r){int mid=l+r>>1; if (clockwise(k,U[mid],U[mid+1])==1) l=mid+1; else ans=mid,r=mid;}
            point w2=U[ans]; return mp(w2,w1);
        } else {
            int l=1,r=where2+1,ans=0;
            while (l<r){int mid=l+r>>1; if (clockwise(k,D[mid],D[mid-1])==-1) ans=mid,l=mid+1; else r=mid;}
            point w1=D[ans]; l=where2,r=D.size()-1,ans=D.size()-1;
            while (l<r){int mid=l+r>>1; if (clockwise(k,D[mid],D[mid+1])==-1) l=mid+1; else ans=mid,r=mid;}
            point w2=D[ans]; return mp(w1,w2);
        }
    }
}
struct P3{
    db x,y,z;
    P3 operator + (P3 k1){return (P3){x+k1.x,y+k1.y,z+k1.z};}
    P3 operator - (P3 k1){return (P3){x-k1.x,y-k1.y,z-k1.z};}
    P3 operator * (db k1){return (P3){x*k1,y*k1,z*k1};}
    P3 operator / (db k1){return (P3){x/k1,y/k1,z/k1};}
    db abs2(){return x*x+y*y+z*z;}
    db abs(){return sqrt(x*x+y*y+z*z);}
    P3 unit(){return (*this)/abs();}
    int operator < (const P3 k1) const{
        if (cmp(x,k1.x)!=0) return x<k1.x;
        if (cmp(y,k1.y)!=0) return y<k1.y;
        return cmp(z,k1.z)==-1;
    }
    int operator == (const P3 k1){
        return cmp(x,k1.x)==0&&cmp(y,k1.y)==0&&cmp(z,k1.z)==0;
    }
    void scan(){
        double k1,k2,k3; scanf("%lf%lf%lf",&k1,&k2,&k3);
        x=k1; y=k2; z=k3;
    }
};
P3 cross(P3 k1,P3 k2){return (P3){k1.y*k2.z-k1.z*k2.y,k1.z*k2.x-k1.x*k2.z,k1.x*k2.y-k1.y*k2.x};}
db dot(P3 k1,P3 k2){return k1.x*k2.x+k1.y*k2.y+k1.z*k2.z;}
//p=(3,4,5),l=(13,19,21),theta=85 ans=(2.83,4.62,1.77)
P3 turn3D(db k1,P3 l,P3 p){
    l=l.unit(); P3 ans; db c=cos(k1),s=sin(k1);
    ans.x=p.x*(l.x*l.x*(1-c)+c)+p.y*(l.x*l.y*(1-c)-l.z*s)+p.z*(l.x*l.z*(1-c)+l.y*s);
    ans.y=p.x*(l.x*l.y*(1-c)+l.z*s)+p.y*(l.y*l.y*(1-c)+c)+p.z*(l.y*l.z*(1-c)-l.x*s);
    ans.z=p.x*(l.x*l.z*(1-c)-l.y*s)+p.y*(l.y*l.z*(1-c)+l.x*s)+p.z*(l.x*l.x*(1-c)+c);
    return ans;
}
typedef vector<P3> VP;
typedef vector<VP> VVP;
db Acos(db x){return acos(max(-(db)1,min(x,(db)1)));}
// 球面距离 , 圆心原点 , 半径 1
db Odist(P3 a,P3 b){db r=Acos(dot(a,b)); return r;}
db r; P3 rnd;
vector<db> solve(db a,db b,db c){
    db r=sqrt(a*a+b*b),th=atan2(b,a);
    if (cmp(c,-r)==-1) return {0};
    else if (cmp(r,c)<=0) return {1};
    else {
        db tr=pi-Acos(c/r); return {th+pi-tr,th+pi+tr};
    }
}
vector<db> jiao(P3 a,P3 b){
    // dot(rd+x*cos(t)+y*sin(t),b) >= cos(r)
    if (cmp(Odist(a,b),2*r)>0) return {0};
    P3 rd=a*cos(r),z=a.unit(),y=cross(z,rnd).unit(),x=cross(y,z).unit();
    vector<db> ret = solve(-(dot(x,b)*sin(r)),-(dot(y,b)*sin(r)),-(cos(r)-dot(rd,b))); 
    return ret;
}
db norm(db x,db l=0,db r=2*pi){ // change x into [l,r)
    while (cmp(x,l)==-1) x+=(r-l); while (cmp(x,r)>=0) x-=(r-l);
    return x;
}
db disLP(P3 k1,P3 k2,P3 q){
    return (cross(k2-k1,q-k1)).abs()/(k2-k1).abs();
}
db disLL(P3 k1,P3 k2,P3 k3,P3 k4){
    P3 dir=cross(k2-k1,k4-k3); if (sign(dir.abs())==0) return disLP(k1,k2,k3);
    return fabs(dot(dir.unit(),k1-k2));
}
VP getFL(P3 p,P3 dir,P3 k1,P3 k2){
    db a=dot(k2-p,dir),b=dot(k1-p,dir),d=a-b;
    if (sign(fabs(d))==0) return {};
    return {(k1*a-k2*b)/d};
}
VP getFF(P3 p1,P3 dir1,P3 p2,P3 dir2){// 返回一条线
    P3 e=cross(dir1,dir2),v=cross(dir1,e);
    db d=dot(dir2,v); if (sign(abs(d))==0) return {};
    P3 q=p1+v*dot(dir2,p2-p1)/d; return {q,q+e};
}
// 3D Covex Hull Template
db getV(P3 k1,P3 k2,P3 k3,P3 k4){ // get the Volume
    return dot(cross(k2-k1,k3-k1),k4-k1);
}
db rand_db(){return 1.0*rand()/RAND_MAX;}
VP convexHull2D(VP A,P3 dir){
    P3 x={(db)rand(),(db)rand(),(db)rand()}; x=x.unit();
    x=cross(x,dir).unit(); P3 y=cross(x,dir).unit();
    P3 vec=dir.unit()*dot(A[0],dir);
    vector<point>B;
    for (int i=0;i<A.size();i++) B.push_back((point){dot(A[i],x),dot(A[i],y)});
    B=ConvexHull(B); A.clear();
    for (int i=0;i<B.size();i++) A.push_back(x*B[i].x+y*B[i].y+vec);
    return A;
}
namespace CH3{
    VVP ret; set<pair<int,int> >e;
    int n; VP p,q;
    void wrap(int a,int b){
        if (e.find({a,b})==e.end()){
            int c=-1;
            for (int i=0;i<n;i++) if (i!=a&&i!=b){
                if (c==-1||sign(getV(q[c],q[a],q[b],q[i]))>0) c=i;
            }
            if (c!=-1){
                ret.push_back({p[a],p[b],p[c]});
                e.insert({a,b}); e.insert({b,c}); e.insert({c,a});
                wrap(c,b); wrap(a,c);
            }
        }
    }
    VVP ConvexHull3D(VP _p){
        p=q=_p; n=p.size();
        ret.clear(); e.clear();
        for (auto &i:q) i=i+(P3){rand_db()*1e-4,rand_db()*1e-4,rand_db()*1e-4};
        for (int i=1;i<n;i++) if (q[i].x<q[0].x) swap(p[0],p[i]),swap(q[0],q[i]);
        for (int i=2;i<n;i++) if ((q[i].x-q[0].x)*(q[1].y-q[0].y)>(q[i].y-q[0].y)*(q[1].x-q[0].x)) swap(q[1],q[i]),swap(p[1],p[i]);
        wrap(0,1);
        return ret;
    }
}
VVP reduceCH(VVP A){
    VVP ret; map<P3,VP> M;
    for (VP nowF:A){
        P3 dir=cross(nowF[1]-nowF[0],nowF[2]-nowF[0]).unit();
        for (P3 k1:nowF) M[dir].pb(k1);
    }
    for (pair<P3,VP> nowF:M) ret.pb(convexHull2D(nowF.se,nowF.fi));
    return ret;
}
//  把一个面变成 ( 点 , 法向量 ) 的形式
pair<P3,P3> getF(VP F){
    return mp(F[0],cross(F[1]-F[0],F[2]-F[0]).unit());
}
// 3D Cut 保留 dot(dir,x-p)>=0 的部分
VVP ConvexCut3D(VVP A,P3 p,P3 dir){
    VVP ret; VP sec;
    for (VP nowF: A){
        int n=nowF.size(); VP ans; int dif=0;
        for (int i=0;i<n;i++){
            int d1=sign(dot(dir,nowF[i]-p));
            int d2=sign(dot(dir,nowF[(i+1)%n]-p));
            if (d1>=0) ans.pb(nowF[i]);
            if (d1*d2<0){
                P3 q=getFL(p,dir,nowF[i],nowF[(i+1)%n])[0];
                ans.push_back(q); sec.push_back(q);
            }
            if (d1==0) sec.push_back(nowF[i]); else dif=1;
            dif|=(sign(dot(dir,cross(nowF[(i+1)%n]-nowF[i],nowF[(i+1)%n]-nowF[i])))==-1);
        }
        if (ans.size()>0&&dif) ret.push_back(ans);
    }
    if (sec.size()>0) ret.push_back(convexHull2D(sec,dir));
    return ret;
}
db vol(VVP A){
    if (A.size()==0) return 0; P3 p=A[0][0]; db ans=0;
    for (VP nowF:A)
        for (int i=2;i<nowF.size();i++)
            ans+=abs(getV(p,nowF[0],nowF[i-1],nowF[i]));
    return ans/6;
}
VVP init(db INF) {
    VVP pss(6,VP(4));
    pss[0][0] = pss[1][0] = pss[2][0] = {-INF, -INF, -INF};
    pss[0][3] = pss[1][1] = pss[5][2] = {-INF, -INF, INF};
    pss[0][1] = pss[2][3] = pss[4][2] = {-INF, INF, -INF};
    pss[0][2] = pss[5][3] = pss[4][1] = {-INF, INF, INF};
    pss[1][3] = pss[2][1] = pss[3][2] = {INF, -INF, -INF};
    pss[1][2] = pss[5][1] = pss[3][3] = {INF, -INF, INF};
    pss[2][2] = pss[4][3] = pss[3][1] = {INF, INF, -INF};
    pss[5][0] = pss[4][0] = pss[3][0] = {INF, INF, INF};
    return pss;
}

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 1
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值