题解-洛谷P6788 「EZEC-3」四月樱花

题面

洛谷P6788 「EZEC-3」四月樱花

给定 (n,p)，求：

[ans=left(prod_{x=1}^nprod_{y|x}frac{y^{d(y)}}{prod_{z|y}(z+1)^2} ight)mod p ]

数据范围：(1le nle 2.5cdot 10^9)，(9.9cdot 10^8<p<1.1cdot 10^9)。

---

蒟蒻语

一道题撑起一场月赛，良心又劲爆。

膜拜出题人 @SOSCHINA，@muxii。

---

蒟蒻解

开局一波猛操作：

[y^{d(y)}=prod_{z|y}y=prod_{z|y}zcdotfrac{y}{z}=prod_{z|y}z^2 ]

[s=prod_{x=1}^nprod_{y|x}frac{y^{d(y)}}{prod_{z|y}(z+1)^2}=prod_{x=1}^nprod_{y|x}prod_{z|y}left(frac{z}{z+1} ight)^2=^{color{#dd6622}{(1)}}prod_{z=1}^nleft(frac{z}{z+1} ight)^{2sumd(lfloorfrac{n}{z} floor)} ]

(color{#dd6622}{(1)}) 的原理就是 (sum_{z|y,y|n}=d(frac{n}{z}))，其中 (sumd(n)=sum_{i=1}^n d(i))。

---

然后就是要求：

[prod_{z=1}^nleft(frac{z}{z+1} ight)^{2sumd(lfloorfrac{n}{z} floor)} ]

很明显 (frac{z}{z+1}) 的前缀积是可以 (Theta(log n)) 求的，问题是怎么求 (sumd(lfloorfrac{n}{z} floor))。

其实是可以分块套分块的，时间复杂度 (Theta(n^{frac 34}+sqrt{n}log n))，勉强卡得过去。

但是有两种时间复杂度 (Theta(n^{frac 23}+sqrt{n}log n)) 的方法：

第一种： 由 @alpha1022 巨佬提供，先筛出 (n^{frac 23}) 的 (sumd)，然后剩下分块套分块。

第二种：

蒟蒻的做法，看到数据范围和 (Theta(n^{frac 23})) 想到杜教筛。

很明显 (d) 这个东西不能直接筛，但是有一个炫酷的魔术：杜教套杜教。

首先 (f=d=1*1)，所以可以令 (g=mu)，(f*g=1*1*mu=1)，满足 (f*g) 前缀和可以速速求，问题是要求 (mu) 的前缀和。

于是再来一次：(f=mu)，(g=1)，(f*g=mu*1=epsilon)，就是杜教筛模板，随意筛。

至于具体怎么套可以看代码，考虑到这题只需要求 (n) 及 (n) 的根号分块的前缀和，所以可以预处理形杜教筛。

---

代码

//Data
using mint=unsigned int;
mint n,nn,mod,ans=1;
mint m(mint x){(x>=mod)&&(x%=mod);return x;}
mint p(mint x){(x>=mod-1)&&(x%=(mod-1));return x;}
void mm(mint&x){(x>=mod)&&(x%=mod);}
void pm(mint&x){(x>=mod-1)&&(x%=(mod-1));}
mint mt(mint x,mint y){return 1ll*x*y%mod;}
mint pt(mint x,mint y){return 1ll*x*y%(mod-1);}
mint Pow(mint a,mint x){mint res=1;for(;x;a=mt(a,a),x>>=1)if(x&1) res=mt(res,a);return res;}

//Sieve
const mint N=5841399+1;
bitset<N> np;
vector<mint> prime;
mint mc[N],d[N],mu[N];
void Sieve(){
	np[1]=true,mc[1]=0,d[1]=mu[1]=1;
	R(i,2,nn){
		if(!np[i]) prime.pb(i),mc[i]=1,d[i]=2,mu[i]=mod-2;
		for(mint p:prime){
			if(!(i*p<nn)) break; np[i*p]=true;
			if(i%p==0){mc[i*p]=mc[i]+1,d[i*p]=d[i]/(mc[i]+1)*(mc[i*p]+1),mu[i*p]=0;break;}
			d[i*p]=d[i]*d[p],mu[i*p]=pt(mu[i],mu[p]),mc[i*p]=1;
		}
	}
	R(i,2,nn) pm(d[i]+=d[i-1]),pm(mu[i]+=mu[i-1]);
}

//DuSieve
const mint iN=427+1;
mint dud[iN],dumu[iN];
bitset<iN> vis;
mint D(mint i){return i<nn?d[i]:dud[n/i];}
mint Mu(mint i){return i<nn?mu[i]:dumu[n/i];}
void DuSieve(mint i){
	if(i<nn||vis[n/i]) return; vis[n/i]=true;
	for(mint l=1,r;l<=i;l=r+1) r=i/(i/l),DuSieve(i/l),
	pm(dumu[n/i]+=p(mod-1-pt(p(r-l+1),Mu(i/l)))); pm(dumu[n/i]+=1);
	for(mint l=2,r;l<=i;l=r+1) r=i/(i/l),
	pm(dud[n/i]+=p(mod-1-pt(p(mod-1+Mu(r)-Mu(l-1)),D(i/l)))); pm(dud[n/i]+=p(i));
	//必须先筛 mu，筛 d 时会用到 mu(i)
}

//Main
int main(){
	read(n),read(mod),nn=1+pow(n,0.72),Sieve(),DuSieve(n); //事实证明0.72最快，可以卡进1s
	for(mint l=1,r;l<=n;l=r+1) r=n/(n/l),ans=mt(ans,Pow(mt(m(l),Pow(m(r+1),mod-2)),D(n/l)));
	write(mt(ans,ans)),putchar(10);
	return 0;
}

---

祝大家学习愉快！