题解 P7575【「PMOI-3」公约数】

· · 题解

P7575 「PMOI-3」公约数解题报告：

更好的阅读体验

题意

给定n,m以及一个长度为n-1的序列x（保证x_i互不相同），求：

\sum_{i_1=1}^m\sum_{i_2=1}^m\cdots\sum_{i_n=1}^m[\gcd(i_1,i_2)=x_1][\gcd(i_2,i_3)=x_2]\cdots[\gcd(i_{n-1},i_n)=x_{n-1}]

分析

好久没有反演了，反演一波。

这个式子看上去要硬上莫反，但由于x_i互不相同，因此不好直接反演。

考虑dp，设f_{i,j}表示考虑前i个数，且第i个数强制设为j的方案数，那么最后的答案为\sum_{i=1}^m f_{n,i}，且很容易列出转移方程：

f_{i,j}=\sum_{k=1}^m[\gcd(j,k)=x_i]f_{i-1,k}

这个式子可以上莫反了，化一化就可以得到：

f_{i,j\times x_i}=\sum_{k=1}^{\lfloor\frac{m}{x_i}\rfloor}f_{i-1,k\times x_i}\sum_{r\mid\gcd(j,k)}\mu(r)=\sum_{r\mid j}\mu(r)\sum_{k=1}^{\lfloor\frac{m}{x_i\times r}\rfloor}f_{i-1,k\times x_i\times r}

设g_{i,j}=\sum_{k=1}^{\lfloor\frac{m}{j}\rfloor}f_{i,k\times j}，那么上式可化简：

f_{i,j\times x_i}=\sum_{r\mid j}\mu(r)\times g_{i-1,x_i\times r}

式子仍然棘手，再次设f'_{i,j}=f_{i,j\times x_i},g'_{i,j}=g_{i,j\times x_i}，那么会有两个式子成立（分别从上面定义式与原式的化简式得来）：

g'_{i,j}=\sum_{j\mid k}^{\lfloor\frac{m}{x_i}\rfloor}f'_{k}\\f'_{i,j}=\sum_{r\mid j}\mu(r)\times g'_{i-1,r}

设h_{i,j}=\mu(j)\times g_{i,j}，那么我们惊奇地发现g'是f'的狄利克雷后缀和，f'是h的狄利克雷前缀和。

我们知道最开始的f_1，可以直接推出f'_1，然后可以通过狄利克雷后缀和推出g'_1，然后用h的定义推出h_1，再次用狄利克雷前缀和推出f'，然后推出f_2……

分析一下时间复杂度，由于狄利克雷前/后缀和的复杂度是O(n\log \log n)的，而我们进行狄利克雷前/后缀和的总数组大小为\sum_{i=1}^n \frac{m}{x_i}\leqslant \sum_{i=1}^m\frac{m}{i}=m\ln m（由于x_i互不相同），所以复杂度应该为O(m\log m\log \log m)。（n与m同阶）

代码

代码具体内容不需要定义这么多数组，两个就够用了。

目前最优解，欢迎吊打。

#include<stdio.h>
const int maxn=1000005,mod=998244353;
int n,m,cnt,ans;
int p[maxn],c[maxn],miu[maxn],x[maxn],f[maxn],g[maxn];
inline int add(int x,int y){
    return x+y>=mod? x+y-mod:x+y;
}
void sieve(int n){
    c[1]=miu[1]=1;
    for(int i=2;i<=n;i++){
        if(c[i]==0)
            p[++cnt]=i,miu[i]=mod-1;
        for(int j=1;j<=cnt;j++){
            if(i*p[j]>n)
                break;
            c[i*p[j]]=1;
            if(i%p[j]==0){
                miu[i*p[j]]=0;
                break;
            }
            miu[i*p[j]]=miu[i]==0? 0:(mod-miu[i]);
        }
    }
}
int main(){
    scanf("%d%d",&n,&m),sieve(m);
    x[0]=1;
    for(int i=1;i<n;i++)
        scanf("%d",&x[i]);
    for(int i=1;i<=m;i++)
        f[i]=1;
    for(int i=1;i<n;i++){
        for(int j=1;j*x[i]<=m;j++)
            g[j]=f[j*x[i]];
        for(int j=1;j*x[i-1]<=m;j++)
            f[j*x[i-1]]=0;
        for(int j=1;j<=cnt&&p[j]<=m/x[i];j++)
            for(int k=(m/x[i])/p[j];k>=1;k--)
                g[k]=add(g[k],g[k*p[j]]);
        for(int j=1;j*x[i]<=m;j++)
            g[j]=1ll*miu[j]*g[j]%mod;
        for(int j=1;j<=cnt&&p[j]<=m/x[i];j++)
            for(int k=1;k*p[j]<=m/x[i];k++)
                g[k*p[j]]=add(g[k*p[j]],g[k]);
        for(int j=1;j*x[i]<=m;j++)
            f[j*x[i]]=g[j];
    }
    for(int i=1;i<=m;i++)
        ans=add(ans,f[i]);
    printf("%d\n",ans);
    return 0;
}