LGP6788题解
太慢了！太慢了！我的替身 【The World】 是最强的替身！
\(O(n^{\frac 2 3})\) 的解法！不清楚用 sbt 能不能更快一些，可能会吧。灵感来源于BZOJ4176，同时也可看到我也是 BZOJ4176 的最优解。理所当然地，我也是 P6788 的最优解
首先看着这个柿子：
\[\prod_{i=1}^n\prod_{d|i}\frac {d^{\sigma_0(d)}} {\prod_{k|d} (k+1)^2}
\]
考虑将 \(d^{\sigma_0(d)}\) 展开：
\[\prod_{i=1}^n\prod_{d|i}\frac {\prod_{k|d}d} {\prod_{k|d}(k+1)^2}
\]
设 \(kx=d\)：
\[\prod_{i=1}^n\prod_{d|i}\prod_{k|d}\frac {kx} {(k+1)^2}
\]
容易发现 \(\prod_{k|d}k \times (\frac d k)=\prod_{k|d} k \times \prod_{k|d}k=(\prod_{k|d}k)^2\)
所以：
\[(\prod_{i=1}^n\prod_{d|i}\prod_{k|d}\frac k {k+1})^2
\]
\[(\prod_{k=1}^n(\frac k {k+1})^{\prod_{i=1}^{\lfloor \frac n k \rfloor}\lfloor \frac n {ki} \rfloor})^2
\]
大多数同学是使用整除分块暴力计算 \(\sum_{i=1}^n \lfloor \frac n i \rfloor\) 而达到 \(O(n^{\frac 3 4})\) 的复杂度，但是这玩意儿其实有性质。
\[\sum_{i=1}^n\sigma_0(i) = \sum_{i=1}^n \sum_{d|i}1 = \sum_{d=1}^n \lfloor \frac n d \rfloor
\]
所以这玩意儿相当于求 \(\sigma_0\) 的块筛。求块筛的常见做法是使用杜教筛或挖掘性质，这里考虑杜教筛。
因为 \(\sigma_0 = 1 * 1\)，所以考虑配对一个 \(\mu\) 上去，使其变为 \(1\)。
只需要同时筛 \(mu\) 和 \(\sigma_0\) 即可。
没有必要对 \(n^{\frac 2 3}\) 以下的 \(\sigma_0\) 使用整除分块计算前缀和，因为在筛 \(\mu\) 的同时把 \(\sigma_0\) 也筛了，这样反而会增加常数。同样也没有必要使用 sbt。
upd:我麻烦了，不需要使用杜教筛，小部分线性筛大部分整除分块即可。其他的好像都没这个快。
#include<cstdio>
#include<cmath>
typedef unsigned ui;
typedef unsigned long long ull;
const int M=2e6;
ui n,mod,lim,top,pri[149000],idx[M+5],d[M+5];
double inv[200005];
inline ui Pow(ui a,ui b=mod-2){
ui ans(1);
for(;b;b>>=1,a=1ull*a*a%mod)if(b&1)ans=1ull*ans*a%mod;
return ans;
}
inline void sieve(const ui&M){
register ui i,j,x;d[1]=1;
for(i=2;i<=M;++i){
if(!d[i])pri[++top]=i,d[i]=2,idx[i]=1;
for(j=1;j<=top&&(x=i*pri[j])<=M;++j){
if(!(i%pri[j])){
idx[x]=idx[i]+1;d[x]=((ui)(d[i]*inv[idx[x]]))*(idx[x]+1);break;
}
idx[x]=1;d[x]=d[i]*2;
}
}
for(i=1;i<=M;++i)d[i]+=d[i-1];
}
inline ui GetSd(const ui&n){
if(n<=lim)return d[n];
ull ans(0);ui i;
for(i=1;i*i<=n;++i)ans+=n*inv[i];
return ((ans<<1)-(i-1)*(i-1))%(mod-1);
}
signed main(){
register ui i,x,L,R,ans=1;
scanf("%u%u",&n,&mod);
for(i=1;(i-1)*(i-1)<=n;++i)inv[i]=1./i+1e-15;sieve(lim=ceil(pow(n,2./3)));
for(L=1;L*L<=n;++L)ans=1ull*ans*Pow(1ull*L*Pow(L+1)%mod,GetSd(n*inv[L]))%mod;
for(x=n*inv[L];L<=n;L=R+1){
if(x*L>n)--x;R=n*inv[x];
ans=1ull*ans*Pow(1ull*L*Pow(R+1)%mod,GetSd(x))%mod;
}
printf("%u",1ull*ans*ans%mod);
}
