題面

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

給定 \(n,p\)，求：

\[ans=\left(\prod_{x=1}^n\prod_{y|x}\frac{y^{d(y)}}{\prod_{z|y}(z+1)^2}\right)\bmod p \]

資料範圍：\(1\le n\le 2.5\cdot 10^9\)，\(9.9\cdot 10^8<p<1.1\cdot 10^9\)。

蒟蒻語

一道題撐起一場月賽，良心又勁爆。

致謝出題人 @SOSCHINA，@muxii。

蒟蒻解

開局一波猛操作：

\[y^{d(y)}=\prod_{z|y}y=\prod_{z|y}z\cdot\frac{y}{z}=\prod_{z|y}z^2 \]

\[s=\prod_{x=1}^n\prod_{y|x}\frac{y^{d(y)}}{\prod_{z|y}(z+1)^2}=\prod_{x=1}^n\prod_{y|x}\prod_{z|y}\left(\frac{z}{z+1}\right)^2=^{\color{#dd6622}{(1)}}\prod_{z=1}^n\left(\frac{z}{z+1}\right)^{2sumd(\lfloor\frac{n}{z}\rfloor)} \]

\(\color{#dd6622}{(1)}\) 的原理就是 \(\sum_{z|y,y|n}=sumd(\lfloor\frac{n}{z}\rfloor)\)

然後就是要求：

\[\prod_{z=1}^n\left(\frac{z}{z+1}\right)^{2sumd(\lfloor\frac{n}{z}\rfloor)} \]

很明顯 \(\frac{z}{z+1}\) 的字首積是可以 \(\Theta(\log n)\) 求的，問題是怎麼求 \(sumd(\lfloor\frac{n}{z}\rfloor)\)。

其實是可以分塊套分塊的，時間複雜度 \(\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=2e6+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=2e3+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));
}

//Main
int main(){
	read(n),read(mod),nn=1+pow(n,0.67),Sieve(),DuSieve(n);
	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;
}

祝大家學習愉快！