【省选模拟】20/04/07 羊（莫比乌斯反演）（杜教筛 / min25）

最新推荐文章于 2022-06-23 22:46:43 发布

FSYo

最新推荐文章于 2022-06-23 22:46:43 发布

阅读量203

点赞数

分类专栏：筛法杜教筛莫比乌斯反演

本文链接：https://blog.csdn.net/sslz_fsy/article/details/105368765

版权

莫比乌斯反演同时被 3 个专栏收录

28 篇文章 0 订阅

订阅专栏

筛法

9 篇文章 0 订阅

订阅专栏

杜教筛

9 篇文章 1 订阅

订阅专栏

今天被打爆了。。。题解单独写

题意：求
$\sum_{k=1}^n\sum_{i=1}^k\sum_{j=1}^kgcd(i,j,k)$
简单反演
$Ans=\sum_{k}\sum_{d|k}d\sum_{i}\sum_{j}[ gcd(i,j,k)=d]\\=\sum_k\sum_{d|k}d\sum_i^{k/d}\sum_{j}^{k/d}[gcd(i,j,k/d)=1]\\=\sum_k\sum_{d|k}d\sum_{l|\frac{k}{d}}\mu(l)(\frac{k}{dl})^2\\=\sum_k\sum_{T|k}(\frac{k}{T})^2\sum_{l|T}d\mu(\frac{T}{l})\\ \sum_k\sum_{T|k}(\frac{k}{T})^2\varphi(T)\\=\sum_{T}T^2\sum_{k}^{n/T}\varphi(k)$
直接杜教筛，我这辈子干过最蠢的事就是考试的时候在构造 $\varphi * (Id\cdot Id)$ 的杜教筛
为了雪耻，用 $m i n 25$ 筛了一下，考虑质数和质数次幂的取值即可，最慢的不到 $0.1 s$

#include<bits/stdc++.h>
#define cs const
#define pb push_back
#define poly vector<int> 
using namespace std;
int read(){
	int cnt=0, f=1; char ch=0;
	while(!isdigit(ch)){ ch=getchar(); if(ch=='-') f=-1; }
	while(isdigit(ch)) cnt=cnt*10+(ch-'0'), ch=getchar();
	return cnt*f;
}
typedef long long ll;
int Mod;
int add(int a, int b){ return a + b >= Mod ? a + b - Mod : a + b; }
int dec(int a, int b){ return a - b < 0 ? a - b + Mod : a - b; }
int mul(int a, int b){ return 1ll * a * b % Mod; }
void Add(int &a, int b){ a = add(a,b); }
void Dec(int &a, int b){ a = dec(a,b); }
void Mul(int &a, int b){ a = mul(a,b); }
int sgn(int a){ return a & 1 ? Mod - 1 : 1; }
int ksm(int a, int b){ int as=1; for(;b;b>>=1,a=mul(a,a)) if(b&1) as=mul(as,a); return as; }
cs int N = 1e5 + 50;
int n, iv6; bool isp[N]; int prm[N], pc;
vector<int> coef[N];
int Sp[3][N], f[3][N]; ll S[N];
int clc(int r){ return mul(mul(r,r+1),mul(r+r+1,iv6)); }
poly operator * (poly a, poly b){
	int deg=a.size()+b.size()-1; poly c(deg,0);
	for(int i=0; i<(int)a.size(); i++)
	for(int j=0; j<(int)b.size(); j++)
	Add(c[i+j],mul(a[i],b[j])); return c;
}
void linear_sieve(int n){
	for(int i=2; i<=n; i++){
		if(!isp[i]){
			prm[++pc] = i;
			Sp[0][pc]=Sp[0][pc-1]+1;
			Sp[1][pc]=Sp[1][pc-1]+i;
			Sp[2][pc]=add(Sp[2][pc-1],mul(i,i));
		}
		for(int j=1; j<=pc; j++) if(prm[j]*i>n) break;
		else{ isp[prm[j]*i]=true; if(i%prm[j]==0) break; }
	} 
}
int sz, nd, A[N], B[N], vl[N]; bool vs[N];
int idx(int x){ return x<=sz?A[x]:B[n/x]; }
int work(int x, int n){
	if(n<=1||prm[x]>n) return 0;
	int as = dec(f[2][idx(n)],Sp[2][x-1]);
	Add(as, dec(f[1][idx(n)],Sp[1][x-1]));
	Dec(as, dec(f[0][idx(n)],Sp[0][x-1]));
	for(int i=x; i<=pc; i++) {
		if(prm[i]*prm[i]>n) break;
		for(int e=1,t=prm[i];(ll)t*prm[i]<=n;t*=prm[i],++e){
			Add(as,mul(work(i+1,n/t),coef[i][e]));
			Add(as,coef[i][e+1]);
		}
	} return as;
}
int main(){
	n=read(); Mod=read(); 
	linear_sieve(sz=sqrt(n)); iv6=ksm(6,Mod-2); 
	for(int i=1; i<=pc; i++){
		int deg=0; ll t=1; while(t<=n) t*=prm[i], ++deg;
		poly a(deg+1,0), b(deg+1,0);
		a[0]=1; for(int j=1,coe=1; j<=deg; j++) 
		a[j]=mul(prm[i]-1,coe), Mul(coe,prm[i]);
		for(int j=0,coe=1; j<=deg; j++) b[j]=mul(coe,coe), Mul(coe,prm[i]);
		coef[i]=a*b; coef[i].resize(deg+1); 
	}
	for(int l=1,r; l<=n; l=r+1){
		int v=n/l; r=n/v;
		if(v<=sz) A[v]=++nd; else B[n/v]=++nd; vl[nd]=v;
		f[0][nd]=v-1; f[1][nd]=(((ll)v*(v+1)>>1)-1)%Mod; f[2][nd]=clc(v)-1;
	}
	for(int i=1; i<=pc; i++){
		for(int j=1; j<=nd; j++){
			if(prm[i]*prm[i]>vl[j]) break;
			int k=idx(vl[j]/prm[i]);
			Dec(f[0][j],dec(f[0][k],Sp[0][i-1]));
			Dec(f[1][j],mul(prm[i],dec(f[1][k],Sp[1][i-1])));
			Dec(f[2][j],mul(mul(prm[i],prm[i]),dec(f[2][k],Sp[2][i-1])));
		}
	} 
	cout << add(1,work(1,n)); return 0;
}