题解 P5221【Product】

· · 题解

一次偶然听到大家说反演就是过不了这题？

出题人标程是反演，但是是 O(n\log n) 的？

有 O(n) 的做法，但不是反演？

…………………………

那我来写 O(n) 的反演做法吧。

由 \operatorname{lcm}(i,j)=\frac{ij}{\gcd(i,j)} 得：

\prod\limits^n_{i=1}\prod\limits^n_{j=1}\frac{ij}{\gcd^2(i,j)} \frac{\prod\limits^n_{i=1}\prod\limits^n_{j=1}ij}{\prod\limits^n_{i=1}\prod\limits^n_{j=1}\gcd^2(i,j)}

分子就是 (n!)^2。分母是下面这个式子的平方：

\prod\limits^n_{i=1}\prod\limits^n_{j=1}\gcd(i,j)

枚举 \gcd，有多少个数对 (i,j) 满足 \gcd(i,j)=d，那么答案的贡献就是 d 的多少次方。

\prod\limits^n_{d=1}d^{\sum\limits^n_{i=1}\sum\limits^n_{j=1}[\gcd(i,j)=d]}

枚举 id 和 jd：

\prod\limits^n_{d=1}d^{\sum\limits^{\lfloor n/d\rfloor}_{i=1}\sum\limits^{\lfloor n/d\rfloor}_{j=1}[\gcd(i,j)=1]}

由 \sum\limits_{g|n}\mu(g)=[n=1] 得：

\prod\limits^n_{d=1}d^{\sum\limits^{\lfloor n/d\rfloor}_{i=1}\sum\limits^{\lfloor n/d\rfloor}_{j=1}\sum\limits_{g|\gcd(i,j)}\mu(g)} \prod\limits^n_{d=1}d^{\sum\limits^{\lfloor n/d\rfloor}_{g=1}\mu(g)\lfloor\frac{n}{dg}\rfloor^2}

这个指数可以整除分块。（注意，mod 是个质数，根据费马小定理，计算指数时应该对 mod-1 取模）

出题人的做法是一个个枚举 d 然后对每个 d 整除分块。

但是我们可以注意到，指数中 \sum 的上界是 \lfloor\frac{n}{d}\rfloor，明显可以对 d 再来一次整除分块！

然后开心的把程序交了上去……

MLE？？？！！！

这就是出题人恶毒所在了。如果我们用两次整除分块，就需要再多一个阶乘数组，分分钟被卡……

那也没事。用尺取法的思想，可以弄两个指针，每次分块的 [l,r] 变化时，指针就移动到相应的位置，边移边对阶乘做乘法。时间复杂度 O(n)。

还是有点说不清楚……看代码吧。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn=1000100,mod=104857601;
#define FOR(i,a,b) for(int i=(a);i<=(b);i++)
#define ROF(i,a,b) for(int i=(a);i>=(b);i--)
#define MEM(x,v) memset(x,v,sizeof(x))
inline int read(){
    char ch=getchar();int x=0,f=0;
    while(ch<'0' || ch>'9') f|=ch=='-',ch=getchar();
    while(ch>='0' && ch<='9') x=x*10+ch-'0',ch=getchar();
    return f?-x:x;
}
int n,mu[maxn],pr[maxn/10],pl;
bool vis[maxn];
void init(){
    mu[1]=1;
    FOR(i,2,n){
        if(!vis[i]) mu[pr[++pl]=i]=-1;
        FOR(j,1,pl){
            if(i*pr[j]>n) break;
            vis[i*pr[j]]=true;
            if(i%pr[j]) mu[i*pr[j]]=-mu[i];
            else{mu[i*pr[j]]=0;break;}
        }
    }
    FOR(i,2,n) mu[i]+=mu[i-1];
}
int calc(int n){    //内层整除分块，上界是n
    int ans=0;
    for(int l=1,r;l<=n;l=r+1){
        r=n/(n/l);
        ans+=1ll*(n/l)*(n/l)%(mod-1)*(mu[r]-mu[l-1]+mod-1)%(mod-1);
        ans%=mod-1; //注意对mod-1取模！
    }
    return ans;
}
int qpow(int a,int b){
    int ans=1;
    for(;b;b>>=1,a=1ll*a*a%mod) if(b&1) ans=1ll*ans*a%mod;
    return ans;
}
int main(){
    n=read();
    init();
    int ans=1;
    FOR(i,1,n) ans=1ll*ans*i%mod;
    ans=qpow(ans,2*n);
    int f1=1,f2=1,c1=1,c2=1;    //两个阶乘和两个指针
    for(int l=1,r;l<=n;l=r+1){  //外层整除分块
        r=n/(n/l);
        while(c1<=r) f1=1ll*f1*c1%mod,c1++; //把指针移到相应的位置
        while(c2<=l-1) f2=1ll*f2*c2%mod,c2++;   //同时修改阶乘
        ans=1ll*ans*qpow(1ll*qpow(f1,mod-2)*f2%mod,2*calc(n/l))%mod;
        //答案除以l到r所有数的积的2*calc(n/l)次方
    }
    printf("%d\n",ans);
}