模板库-数学和计算几何

数学

高斯消元

需要求逆元：

inline ModInt det(int n,ModInt (*a)[N]){
	ModInt ans;ans=1;
	for(int i=1;i<=n;i++){
		if(!a[i][i].x){
			ans=0-ans;
			for(int k=i+1;k<=n;k++)if(a[k][i].x){std::swap(a[i],a[k]);goto CONTINUE;}
			return {0};
		}
	CONTINUE:
		ans*=a[i][i];
		for(int k=i+1;k<=n;k++)if(a[k][i].x){
			ModInt o=a[k][i]/a[i][i];
			for(int j=i;j<=n;j++) a[k][j]-=a[i][j]*o;
		}
	}
	return ans;
}

不需要逆元：

inline ModInt det(){
	ModInt ans;ans=1;
	for(int i=1;i<=n;i++){
		for(int j=i+1;j<=n;j++){
			while(a[j][i]){
				long long k=a[i][i]/a[j][i];
				for(int h=1;h<=n;h++) a[i][h]-=k*a[j][h];
				std::swap(a[i],a[j]);
				ans=-ans;
			}
		}
	}
	for(int i=1;i<=n;i++) ans*=a[i][i];
	return ans;
}

矩阵求逆（高斯消元法）

inline int work(){
	for(int i=1;i<=n;i++){
		if(!a[i][i]){
			for(int k=i+1;k<=n;k++)if(a[k][i]){
				std::swap(a[k],a[i]);goto CONTINUE;
			}
			return 0;
		}
CONTINUE:
		long long inv=power(a[i][i],MOD-2);
		for(int k=i+1;k<=n;k++)if(a[k][i]){
			long long o=a[k][i]*inv%MOD;
			for(int j=i;j<=n<<1;j++) a[k][j]-=a[i][j]*o%MOD,a[k][j]=(a[k][j]+MOD)%MOD;
		}
	}
	for(int i=n;i;i--){
		if(!a[i][i]){
			for(int k=i-1;k;k--)if(a[k][i]){
				std::swap(a[k],a[i]);break;
			}
		}
		long long inv=power(a[i][i],MOD-2);
		for(int k=1;k<i;k++)if(a[k][i]){
			long long o=a[k][i]*inv%MOD;
			for(int j=i;j<=n<<1;j++) a[k][j]-=a[i][j]*o%MOD,a[k][j]=(a[k][j]+MOD)%MOD;
		}
		for(int j=i;j<=n<<1;j++) a[i][j]=a[i][j]*inv%MOD;
	}
	return 1;
}

快速乘 和 minner-rabbin 和 Pollard-Rho

inline unsigned long long Rand(){
	return (unsigned long long)rand()*RAND_MAX*RAND_MAX+(unsigned long long)rand()*RAND_MAX+
	(unsigned long long)rand()*RAND_MAX*RAND_MAX*RAND_MAX+rand();
}
inline long long mul(long long x,long long y,long long mod){
	return (x*y-(long long)((long double)x/mod*y)*mod+mod)%mod;
}
inline long long power(long long a,long long b,long long mod){
	long long ans=1;
	while(b){
		if(b&1) ans=mul(ans,a,mod);
		a=mul(a,a,mod);b>>=1;
	}
	return ans;
}
inline int check(long long x,int p){
	long long d=x-1,r=0;
	while(!(d&1)) d>>=1,r++;
	long long pd=power(p,d,x);
	if(pd==1) return 1;
	while(r--){
		if(pd==x-1) return 1;
		pd=mul(pd,pd,x);
	}
	return 0;
}
const int P[]={2,3,5,7,11,13,17,37};
inline int isPrime(long long x){
	if(x<2) return 0;
	for(int i=0;i<8;i++){
		if(x==P[i]) return 1;
		if(!(x%P[i])) return 0;
		if(!check(x,P[i])) return 0;
	}
	return 1;
}
inline long long getGcd(long long a,long long b){return b?getGcd(b,a%b):a;}
long long divide(long long x){
	if(isPrime(x)) return x;
	long long T;
	int k=log2(x),o=k;
	for(;;){
		T=Rand()%19260917+1;
		long long v0=Rand()%(x-1)+1,v=(mul(v0,v0,x)+T)%x,sum=1;
		while(v0^v){
			sum=mul(sum,lib::abs(v-v0),x);
			if(!sum) break;
			if(!o--){
				long long d=getGcd(sum,x);
				if(d>1&&d<x) return lib::max(divide(d),divide(x/d));
				o=k;
			}
			v=(mul(v,v,x)+T)%x;v=(mul(v,v,x)+T)%x;v0=(mul(v0,v0,x)+T)%x;
		}
	}
}

ntt

常用素数：
P = 1004535809  ====>  pr = 3
P = 998244353  =====>  pr = 3

 

\(n\) 存的是项数而不是次数

inline int init(int n){
	int max=1;while(max<n) max<<=1;
	for(int i=0;i<max;i++) rev[i]=rev[i>>1]>>1,rev[i]|=(i&1)?(max>>1):0;
	return max;
}
inline void ntt(int n,long long *a,int type){
	for(int i=0;i<n;i++)if(rev[i]<i) lib::swap(a[i],a[rev[i]]);
	for(int h=1;h<n;h<<=1){
		long long gn=power(G,(mod-1)/(h<<1)),g,o;
		if(!type) gn=power(gn,mod-2);
		for(int i=0,j;i<n;i+=h<<1){
			for(g=1,j=i;j<i+h;j++,g=g*gn%mod){
				o=g*a[j+h]%mod;
				a[j+h]=(a[j]-o+mod)%mod;a[j]=(a[j]+o)%mod;
			}
		}
	}
	if(!type){
		long long inv=power(n,mod-2);
		for(int i=0;i<n;i++) a[i]=a[i]*inv%mod;
	}
}

fwt

inline void Or(int n,long long *f){for(int h=1;h<n;h<<=1)for(int i=0;i<n;i+=h<<1)for(int j=0;j<h;j++) f[i+j+h]+=f[i+j];}
inline void IOr(int n,long long *f){for(int h=1;h<n;h<<=1)for(int i=0;i<n;i+=h<<1)for(int j=0;j<h;j++) f[i+j+h]-=f[i+j];}
inline void And(int n,long long *f){for(int h=1;h<n;h<<=1)for(int i=0;i<n;i+=h<<1)for(int j=0;j<h;j++) f[i+j]+=f[i+j+h];}
inline void IAnd(int n,long long *f){for(int h=1;h<n;h<<=1)for(int i=0;i<n;i+=h<<1)for(int j=0;j<h;j++) f[i+j]-=f[i+j+h];}
inline void Xor(int n,long long *f){
	for(int h=1;h<n;h<<=1)for(int i=0;i<n;i+=h<<1)for(int j=0;j<h;j++)
		f[i+j]+=f[i+j+h],f[i+j+h]=f[i+j]-f[i+j+h]-f[i+j+h];
}
inline void IXor(int n,long long *f){
	for(int h=1;h<n;h<<=1)for(int i=0;i<n;i+=h<<1)for(int j=0;j<h;j++)
		f[i+j]+=f[i+j+h],f[i+j+h]=f[i+j]-f[i+j+h]-f[i+j+h],f[i+j]*=inv2,f[i+j+h]*=inv2;
}

lu div

inline void luDiv(int n,ModInt (*l)[N],ModInt (*u)[N],ModInt (*a)[N]){
	for(int i=1;i<=n;i++){
		for(int j=i;j<=n;j++){
			u[i][j]=a[i][j];
			for(int k=1;k<i;k++) u[i][j]-=l[i][k]*u[k][j];
			l[j][i]=a[j][i];
			for(int k=1;k<i;k++) l[j][i]-=l[j][k]*u[k][i];
			l[j][i]/=u[i][i];
		}
		l[i][i]=1;
	}
}

自适应辛普森

inline double simpson(double l,double r){
	return (r-l)*(calc(l)+4*calc((l+r)/2)+calc(r))/6;
}
double asr(double l,double r,double eps,double ans){
	double mid=(l+r)/2;
	double fl=simpson(l,mid),fr=simpson(mid,r);
	if(std::fabs(fl+fr-ans)<=15*eps) return fl+fr+(fl+fr-ans)/15;
	return asr(l,mid,eps/2,fl)+asr(mid,r,eps/2,fr);
}

计算几何

lib

struct Vector{
	double x,y;
	inline Vector(){x=y=0;}
	inline Vector(double _x,double _y){x=_x;y=_y;}
	inline double len(){return __builtin_sqrt(x*x+y*y);}
	inline double len2(){return x*x+y*y;}
	inline void operator += (const Vector &a){x+=a.x;y+=a.y;}
	inline void operator -= (const Vector &a){x-=a.x;y-=a.y;}
	inline void operator *= (const double &a){x*=a;y*=a;}
	inline void operator /= (const double &a){x/=a;y/=a;}
	inline Vector operator + (const Vector &a)const{return {x+a.x,y+a.y};}
	inline Vector operator - (const Vector &a)const{return {x-a.x,y-a.y};}
	inline Vector operator * (const double &a)const{return {x*a,y*a};}
	inline Vector operator / (const double &a)const{return {x/a,y/a};}
	inline double operator * (const Vector &a)const{return x*a.x+y*a.y;}
	inline double operator ^ (const Vector &a)const{return x*a.y-y*a.x;}
};
struct Line{
	Vector way,o;
	inline Line(){}
	inline Line(double a,double b,double c){//ax+by+c=0
		if(std::abs(a)>=eps) o=Vector(-c/a,0);
		else o=Vector(0,-c/b);
		way=Vector(b,-a);
	}
	inline Line(const Vector &a,const Vector &b){o=a;way=b-a;}
};
inline Vector getIntersection(const Line &a,const Line &b){
	double x=(b.way^(a.o-b.o))/(a.way^b.way);
	return a.o+a.way*x;
}

二维凸包

inline int graham(int n,Vector *p,Vector *stack){
	for(int i=2;i<=n;i++) p[i].ang=atan2(p[i].y-p[1].y,p[i].x-p[1].x);
	std::sort(p+2,p+1+n);
	stack[1]=p[1];int top=1;
	for(int i=2;i<=n;i++){
		while(top>1&&cross(p[i]-stack[top],stack[top]-stack[top-1])<=eps) top--;
		stack[++top]=p[i];
	}
	stack[++top]=p[1];
	return top;
}

半平面交

int left,right;
Line que[N];
inline int halfPlane(int n,Line *a,Vector *p){
	std::sort(a+1,a+1+n,cmp);
	left=right=0;que[0]=a[1];
	for(int i=2;i<=n;i++){
		while(left<right&&onRight(a[i],p[right])) right--;
		while(left<right&&onRight(a[i],p[left+1])) left++;
		que[++right]=a[i];
		if(abs(cross(que[right].way,que[right-1].way))<=eps){//平行
			if(onRight(que[right],que[right-1].p)&&dot(que[right].way,que[right-1].way)<=-eps) return 0;
			right--;
			if(!onRight(que[right],a[i].p)) que[right]=a[i];
		}
		if(left<right) p[right]=intersect(que[right],que[right-1]);
	}
	while(left<right&&onRight(que[left],p[right])) right--;
	if(right-left<=1) return 0;
	p[left]=intersect(que[left],que[right]);
	return 1;
}

旋转卡壳

inline int nex(int i){return i==n?1:(i+1);}
int main(){
	if(n<=3) return printf("%lld\n",(q[2]-q[1]).len()),0;
	long long ans=0;
	for(int i=1,j=3;i<n;i++){
		while(abs((q[j]-q[i])^(q[j]-q[i+1]))<=abs((q[nex(j)]-q[i])^(q[nex(j)]-q[i+1]))) j=nex(j);
		ans=std::max(ans,std::max((q[j]-q[i]).len(),(q[j]-q[i+1]).len()));
	}
	printf("%lld\n",ans);
	return 0;
}

字符串

kmp

manacher

AC 自动机

struct Node{
	Node *tr[26],*fail;
	int cnt,in;
}dizhi[N],*root=&dizhi[0],*which[N];
int tot;
inline void insert(char *s,int id){
	Node *tree=root;
	for(int num,i=1;s[i];i++){
		num=s[i]-'a';
		if(!tree->tr[num]) tree->tr[num]=&dizhi[++tot];
		tree=tree->tr[num];
	}
	which[id]=tree;
}
int left,right;
Node *que[N];
inline void build(){
	left=0;right=-1;
	root->fail=root;
	for(int i=0;i<26;i++){
		if(root->tr[i]) root->tr[i]->fail=root,que[++right]=root->tr[i];
		else root->tr[i]=root;
	}
	Node *u;
	while(left<=right){
		u=que[left++];
		for(int i=0;i<26;i++){
			if(u->tr[i]){
				u->tr[i]->fail=u->fail->tr[i];
				que[++right]=u->tr[i];u->fail->tr[i]->in++;
			}
			else u->tr[i]=u->fail->tr[i];
		}
	}
}

SAM

struct Node{
	Node *son[26];
	int len,cnt;
	Node *link;
}dizhi[N*2],*root=&dizhi[0],*last=root;
int tot;
inline void add(int c){
	Node *p=last,*now=&dizhi[++tot];
	now->len=p->len+1;now->cnt=1;
	for(;p&&!p->son[c];p=p->link) p->son[c]=now;
	if(!p) now->link=root;
	else{
		Node *q=p->son[c];
		if(q->len==p->len+1) now->link=q;
		else{
			Node *clone=&dizhi[++tot];*clone=*q;
			clone->len=p->len+1;clone->cnt=0;
			q->link=now->link=clone;
			for(;p&&p->son[c]==q;p=p->link) p->son[c]=clone;
		}
	}
	last=now;
}

PAM

struct Node{
	Node *son[26],*fail;
	int len,cnt;
}dizhi[N],*root0=&dizhi[0],*root1=&dizhi[1],*last;
int tot;
inline void init(){
	last=root1;
	root0->fail=root1;root1->fail=root0;
	root1->len=-1;
	tot=1;
}
char s[N];
inline int cnt(Node *u){
	if(u==root0||u==root1) return 0;
	if(u->cnt) return u->cnt;
	return u->cnt=1+cnt(u->fail);
}
inline int insert(int i){
	while(s[i-last->len-1]^s[i]) last=last->fail;
	if(!last->son[s[i]-'a']){
		Node *u=&dizhi[++tot];
		u->len=last->len+2;
		reg Node *j=last->fail;
		while(s[i-j->len-1]^s[i]) j=j->fail;
		u->fail=j->son[s[i]-'a'];
		if(!u->fail) u->fail=root0;
		last->son[s[i]-'a']=u;
	}
	last=last->son[s[i]-'a'];
	return cnt(last);
}