【模板】数论

FFT https://www.luogu.com.cn/problem/P3723

#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;
const double pi=acos(-1.0);
const int N=401000;
//int a[N],b[N],c[N<<1],rev[N<<3],limit;
int a[N],b[N],c[N],rev[N],limit;
struct comp{
	double x,y;
	comp(double xx=0,double yy=0) { x=xx;y=yy; }
	friend comp operator + (const comp &a,const comp &b) { return comp(a.x+b.x , a.y+b.y); }
	friend comp operator - (const comp &a,const comp &b) { return comp(a.x-b.x , a.y-b.y); }
	friend comp operator * (const comp &a,const comp &b) { return comp(a.x*b.x-a.y*b.y , a.y*b.x + a.x*b.y); }
}A[N],B[N];//A[N<<3],B[N<<3];//limit最大6*N 

int minn(int x,int y) { return x<y?x:y; }
void fft(comp *a,double f)
{
	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)
	{
		comp wn=comp(cos(pi/mid),f*sin(pi/mid));
		for(int L=0,R=(mid<<1);L<limit;L+=R){
			comp w=comp(1.0,0.0);
			for(int k=0;k<mid;k++,w=w*wn)
			{
				comp A1=a[L+k],A2=w*a[L+k+mid];
				a[L+k]=A1+A2;
				a[L+k+mid]=A1-A2;
			}
		}
	}
}
int main()
{
	int n,m,a1=0,a2=0,b1=0,b2=0;scanf("%d %d",&n,&m);
	for(int i=1;i<=n;i++) scanf("%d",&a[i]),a1+=a[i]*a[i],a2+=a[i];
	for(int i=1;i<=n;i++) scanf("%d",&b[i]),b1+=b[i]*b[i],b2+=b[i];
	for(int i=1;i<=n;i++) c[i]=a[n-i+1],b[i+n]=b[i];
	
	int nn=3*n,l=0;//c:2n b:n -- limit:3n
	limit=1;
	while(limit<nn) limit<<=1,l++;
	for(int i=1;i<limit;i++) rev[i]=(rev[i>>1]>>1)|(i&1)<<(l-1);
	for(int i=0;i<limit;i++) A[i]=comp(1.0*c[i],0.0);
	for(int i=0;i<limit;i++) B[i]=comp(1.0*b[i],0.0);
	
	fft(A,1);
	fft(B,1);
	for(int i=0;i<limit;i++) A[i]=A[i]*B[i];
	fft(A,-1);

	int ans=2147483647;
	for(int i=1;i<=n;i++)
		for(int c=-m;c<=m;c++)
			ans=minn(ans,a1+b1+n*c*c + 2*c*(a2-b2) -2*(int)(A[i+n].x/limit+0.5));
	printf("%d",ans);
	return 0;
}

原根NTT（非任意模数） https://www.luogu.com.cn/problem/P3723

#include<cstdio>
#include<cstring>
#include<string>
#include<cmath>
#include<algorithm>
#define mod 998244353
using namespace std;

const int N=3000010;inline char cn(){
	static char buf[1000010],*p1=buf,*p2=buf;
	return p1==p2&&(p2=(p1=buf)+fread(buf,1,1000000,stdin),p1==p2)?EOF:*p1++;
}
#define cn getchar
template<class TY>void in(TY &x){
	x=0;int f1=1;char ch=cn();
	while(ch<'0'||ch>'9'){if(ch=='-')f1=-1;ch=cn();}
	while(ch>='0'&&ch<='9')x=x*10+(ch-'0'),ch=cn(); x*=f1;
}
char OP[1000030],Zhan[30];int op=0,Top;
void fcheck(int type=0){if(type||op>=1000000)fwrite(OP,1,op,stdout),op=0;}
template<class TY>inline void write(TY x){
	if(x==0){OP[op++]='0';fcheck();return;}
	Top=0;while(x)Zhan[++Top]=x%10+'0',x/=10;
	while(Top)OP[op++]=Zhan[Top--];fcheck();
}
#define K_G OP[op++]=' '
//----------------------------------------------------------------------------------
int n,m;
int rev[N],limit=1,l=0;
int ksm(int b,int p){int ans=1;for(;(p&1?ans=1ll*ans*b%mod:0),p;p>>=1,b=1ll*b*b%mod);return ans;}
int w[N],inv[N];
void prep(){
	inv[1]=1;for(int i=2;i<=limit;i++)inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
	for(int mid=1,wn;mid<limit;mid<<=1){
		w[mid]=1;wn=ksm(3,(mod-1)/(mid<<1));
		for(int k=1;k<mid;k++)w[mid+k]=1ll*w[mid+k-1]*wn%mod;
	}
}
int add(int x) {return x>=mod?x-mod:x;}
int dec(int x) {return x<0?x+mod:x;}
void ntt(int *a,int flag){
	if(flag==-1) reverse(a+1,a+limit);
	for(int i=1;i<limit;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
	for(int mid=1;mid<limit;mid<<=1) for(int L=0;L<limit;L+=(mid<<1)) for(int k=0;k<mid;k++)
	{int y=1ll*a[L+k+mid]*w[mid+k]%mod;a[L+k+mid]=dec(a[L+k]-y);a[L+k]=add(a[L+k]+y);}
}
int a[N],b[N];
int main()
{
	in(n),in(m);
	for(int i=0;i<=n;i++) in(a[i]);
	for(int i=0;i<=m;i++) in(b[i]);
	while(limit<n+m+1) limit<<=1,l++;
	for(int i=0;i<limit;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
	prep();
	ntt(a,1);ntt(b,1);
	for(int i=0;i<limit;i++) a[i]=1ll*a[i]*b[i]%mod;
	ntt(a,-1);
	for(int i=0;i<=n+m;i++) a[i]=1ll*a[i]*inv[limit]%mod,write(a[i]),K_G;
	fcheck(1);return 0;
}

多项式求幂

莫比乌斯反演 线性筛（杜教筛在下面)  https://www.luogu.com.cn/problem/P2522

#include<cstdio>
int mu[50100],sum[50100];
int li[50100];//质数队列 
bool p[50100];//质数桶 
void m_mu()//得到莫比乌斯函数 
{
	int l=0,now;
	mu[0]=0,mu[1]=1;
	for(int i=2;i<=50010;i++)
	{
		if(!p[i]) { mu[i]=-1; li[++l]=i; } //一个质数 
		for(int j=1;j<=l;j++)
		{
			now=i*li[j];
			if(now>50010) break;
			p[now]=true;//将i的所有质数倍标记为非质数 
			mu[now]=0;
			if(i%li[j]==0) break; //重复质数,一个数必能被分解质因数 
			else mu[now]=mu[i]*(-1);
		}}
 	for(int i=1;i<=50010;i++) sum[i]=sum[i-1]+mu[i];
}
int t,a,b,c,d,k;
int mymin(int x,int y) { return x<y?x:y; }
int slove(int a,int b) 
{
	/*这里完全可以 a/=k,b/=k;*/
	int minn=mymin(a,b),ans=0; 
	//计算 "sigma_t (n/t)向下取整" ， 显然许多  "(n/t)向下取整"是相同的 
 	//其最后一个数就是n/(n/t)
 	//所以不必枚举每一个t，用l就可以代替一堆的t 
  	for(int l=1,r;l<=minn;l=r+1){
  		r=mymin(a/(a/l),b/(b/l));	
  		ans+=( sum[r]-sum[l-1] ) * ( a/(k*l) ) * ( b/(k*l) );
  	}
 	return ans;
}
int main()
{
	m_mu();
	scanf("%d",&t);
	while(t--){
		scanf("%d %d %d %d %d",&a,&b,&c,&d,&k);
		printf("%d\n",slove(b,d)-slove(a-1,d)-slove(b,c-1)+slove(a-1,c-1));
	} //容斥原理      (1~b,1~d)   (1~a-1,1~d)   (1~b,1~c-1)   (1~a-1, 1~c-1)
	return 0;
}

杜教筛 https://www.luogu.com.cn/problem/P4213

#include <cstdio>
#include <cstring>
#include <map>
#include <algorithm>
using namespace std;
#define N 5000010
#define uint unsigned int
#define LL long long
uint T,n; LL ans1,ans2;
int prime[N],t=0;
LL mu[N];bool v[N];
void work(){//预处理 
	mu[1]=1;
	for(int i=2;i<=N-10;i++){
		if(!v[i]) prime[++t]=i,mu[i]=-1;
		for(int j=1;j<=t && i*prime[j]<=N-10;j++){
			v[i*prime[j]]=true;
			if(i%prime[j]==0) break;
			mu[i*prime[j]]=-mu[i];
		}
		mu[i]+=mu[i-1];
	}
}
map<uint,LL> S;
LL solvemu(uint x){
	if(x<=N-10) return mu[x];
	if(S.count(x)) return S[x];
	LL ans=1;
	for(uint l=2,r;l<=x;l=r+1)
	r=x/(x/l),ans-=(r-l+1)*solvemu(x/l);
	S[x]=ans; return ans;
}
LL solvephi(uint x){
	LL ans=0;
	for(uint l=1,r;l<=x;l=r+1)
	r=x/(x/l), ans+=1LL*(l+r)*(r-l+1)/2*solvemu(x/l);
	return ans;
}
int main()
{
	scanf("%u",&T);work();
	while(T--){
		scanf("%u",&n);
		ans2=solvemu(n);
		ans1=solvephi(n);
		printf("%lld %lld\n",ans1,ans2);
	}
	return 0;
}

矩阵乘法 https://www.luogu.com.cn/problem/P6569
mat(LL t=0){ 
		memset(a,0,sizeof(a)); 
		a[1][1]=a[2][2]=a[3][3]=a[4][4]=a[5][5]=t; 
}t=1是 任一正方形矩阵 与 这行代码相乘等于原矩阵 矩阵快速幂中会用到
	
#include<cstdio>
#include<cstring>
struct mat{
	int a[110][110],n,m;
	mat(){ n=m=0;memset(a,0,sizeof(a)); }
	mat operator * (const mat &B) const{
		mat C; 
		C.n=n; 
		C.m=B.m; //C(n,m)=A(n,k)*B(k,m)
		for (int i=1;i<=n;i++) {
			for (int j=1;j<=B.m;j++) {
				for (int k=1;k<=m;k++) {
					C.a[i][j]^=(a[i][k]*B.a[k][j]);
				}
			}
		}return C;
	}
};
mat ans,ma[35];
int n,m,q,f[110];
int main() 
{
	int m,q,x,y;scanf("%d %d %d",&n,&m,&q);
	for(int i=1;i<=n;i++) scanf("%d",&f[i]);
	for(int i=1;i<=m;i++){
		scanf("%d %d",&x,&y);
		ma[0].a[x][y]=ma[0].a[y][x]=1;
	}
	
	ma[0].n=ma[0].m=n;//预处理
	for(int i=1;i<=31;i++) ma[i]=ma[i-1]*ma[i-1];
	
	ans.n=1,ans.m=n;
	for(int i=1;i<=q;i++){
		scanf("%d",&x);
		for(int j=1;j<=n;j++) ans.a[1][j]=f[j];
		for(int t=0;t<=31;t++)
			if(x&(1<<t)) ans=ans*ma[t];
		printf("%lld\n",ans.a[1][1]);
	}
	return 0;
}