题解：P7504 「HMOI R1」可爱的德丽莎

· · 题解

数论题，先化式子。

\sum_{x = 1}^n\sum_{y = 1}^n\sum_{i = 1}^x[x \bot k_1][i \bot x]\cdot i\cdot \sum_{j = 1}^y[y \bot k_2][j \bot y]\cdot j =\sum_{x = 1}^n[x \bot k_1]\sum_{i = 1}^x[i \bot x]\cdot i\cdot\sum_{y = 1}^n [y \bot k_2]\sum_{j = 1}^y[j \bot y]\cdot j

可以发现左右两个式子是相同的，我们讨论其中一个。

\sum_{x = 1}^n[x \bot k_1]\sum_{i = 1}^x[i \bot x]\cdot i

看到 [i \bot x] 套路的上莫比乌斯反演。

=\sum_{x = 1}^n[x \bot k_1]\sum_{d \mid x}\sum_{i=1}^{ \frac{x}{d} }\mu(d)\cdot i \cdot d

把 i 单独列出去，发现是个求和，可以消去一个 \sum。

=\frac{1}{2} \sum_{x = 1}^n[x \bot k_1]\sum_{d \mid x}\mu(d)\cdot d \cdot \frac{x}{d}\cdot(\frac{x}{d}+1) =\frac{1}{2} \sum_{x = 1}^n[x \bot k_1]\cdot x \cdot \sum_{d \mid x}\mu(d)\cdot(\frac{x}{d}+1)

发现 \sum_{d \mid x} 里面一个是 \mu * id 一个是 \mu * 1，分别是 \varphi 和 \epsilon。

=\frac{1}{2} \sum_{x = 1}^n[x \bot k_1]\cdot x \cdot (\varphi(x)+\epsilon(x))

把 \epsilon 拆出来。

=\frac{1}{2} (1+\sum_{x = 1}^n[x \bot k_1]\cdot x \cdot \varphi(x))

似乎没什么好化的了，现在变成了一个前缀和的形式，于是我们考虑杜教筛。

发现若 g(x)=[x \bot k]\cdot x，则 (f * g)(x)=[x \bot k] \cdot x^2，这个形式看起来很不错。

接下来我们需要求 g(x)=[x \bot k] \cdot x 的前缀和，列出式子：

\sum_{i=1}^{n}[i \bot k] \cdot i =\sum_{i=1}^{n}i\cdot \sum_{d\mid k}\mu (d) =\sum_{d \mid k}\mu(d) \cdot d \cdot \sum_{i=1}^{\lfloor \frac{n}{d} \rfloor}i

发现后面的 \sum 是一个求和，可消去，前面的 \sum 想不到啥方法优化，直接枚举 k 的因数即可。

对于 f*g 的前缀和，与 g 的前缀和同理。

于是我们可以在可接受的时间内计算出 g 和 f*g 的前缀和，接下来杜教筛求出 f 的前缀和即可。

下附代码：

#include<iostream>
#include<cstdio>
#include<unordered_map>
#define int long long
using namespace std;
const int mod=998244353;
int n,k,inv6,mu[1000005],p[1000005],phi[1000005],tot,f[1000005],g[1000005],fg[1000005],v[1000005],mul[1000005],cur;
long long ans=1;
bool isp[1000005];
inline int f_pow(long long a,int b){
    long long res=1;
    while(b>0){
        if(b&1) res*=a,res%=mod;
        a*=a,a%=mod,b>>=1;
    }
    return res;
}
inline int Sum1(int x){
    return 1ll*x*(x+1)/2%mod;
}
inline int Sum2(int x){
    return 1ll*x*(x+1)%mod*(2*x+1)%mod*inv6%mod;
}
inline int getmu(int x){
    int res=0;
    for(int i=2;i<=x/i;i++){
        if(x%i==0){
            int cnt=0;
            while(x%i==0) cnt++,x/=i;
            if(cnt>=2) return 0;
            res++;
        }
    }
    if(x>1) res++;
    return (res&1?-1:1);
}
unordered_map<int,int>mpf,mpg,mpfg;
inline void getprime(){
    tot=0;mu[1]=1;phi[1]=1;
    for(int i=2;i<=1000000;i++) isp[i]=1;
    for(int i=2;i<=1000000;i++){
        if(isp[i]) p[++tot]=i,mu[i]=-1,phi[i]=i-1;
        for(int j=1;j<=tot&&p[j]*i<=1000000;j++){
            isp[i*p[j]]=0;
            mu[i*p[j]]=mu[i]*mu[p[j]];
            phi[i*p[j]]=phi[i]*phi[p[j]];
            if(i%p[j]==0){
                mu[i*p[j]]=0;
                phi[i*p[j]]=phi[i]*p[j];
                break;
            }
        }
    }
}
inline void getfg(){
    cur=0;f[1]=1,g[1]=1,fg[1]=1;
    for(int i=2;i<=1000000;i++) f[i]=0,g[i]=0,fg[i]=0;
    for(int i=2;i<=1000000;i++){
        if(isp[i]&&k%i!=0) f[i]=1ll*i*phi[i]%mod,g[i]=i,fg[i]=1ll*i*i%mod;
        for(int j=1;j<=tot&&p[j]*i<=1000000;j++){
            if(f[i]&&k%p[j]!=0) f[i*p[j]]=1ll*i*p[j]*phi[i*p[j]]%mod;
            g[i*p[j]]=1ll*g[i]*g[p[j]]%mod;
            fg[i*p[j]]=1ll*fg[i]*fg[p[j]]%mod;
            if(i%p[j]==0){
                break;
            }
        }
    }
    for(int i=1;i<=1000000;i++){
        f[i]+=f[i-1];f[i]%=mod;
        g[i]+=g[i-1];g[i]%=mod;
        fg[i]+=fg[i-1];fg[i]%=mod;
    }
    for(int i=1;i<=k/i;i++){
        if(k%i==0){//提前记录 k 的因数
            v[++cur]=i;mul[cur]=mu[i];
            if(1ll*i*i!=k) v[++cur]=k/i,mul[cur]=getmu(k/i);
        }
    }
}
inline int getgsum(int n){
    if(n<=1000000) return g[n];
    if(mpg.find(n)!=mpg.end()) return mpg[n];
    long long sum=0;
    for(int i=1;i<=cur;i++){
        sum+=1ll*mul[i]*v[i]*Sum1(n/v[i])%mod;
        // sum=(sum%mod+mod)%mod;
    }
    sum=(sum%mod+mod)%mod;
    mpg[n]=sum;
    return sum;
}
inline int getfgsum(int n){
    if(n<=1000000) return fg[n];
    if(mpfg.find(n)!=mpfg.end()) return mpfg[n];
    long long sum=0;
    for(int i=1;i<=cur;i++){
        sum+=1ll*mul[i]*v[i]*v[i]%mod*Sum2(n/v[i])%mod;
        // sum=(sum%mod+mod)%mod;
    }
    sum=(sum%mod+mod)%mod;
    mpfg[n]=sum;
    return sum;
}
inline int getfsum(int n){
    if(n<=1000000) return f[n];
    if(mpf.find(n)!=mpf.end()) return mpf[n];
    long long sum=getfgsum(n);
    int l=2,r=0;
    while(l<=n){
        r=n/(n/l);
        sum-=1ll*(getgsum(r)-getgsum(l-1))*getfsum(n/l)%mod;
        // sum=(sum%mod+mod)%mod;
        l=r+1;
    }
    sum=(sum%mod+mod)%mod;
    mpf[n]=sum;
    return sum;
}
inline int crychic(){
    mpf.clear();mpg.clear();mpfg.clear();getfg();
    return 1ll*(getfsum(n)+1)*f_pow(2,mod-2)%mod;
}
signed main(){
    cin>>n>>k;inv6=f_pow(6,mod-2);getprime();
    ans*=crychic();
    cin>>k;ans*=crychic();ans%=mod;
    cout<<ans<<endl;
    return 0;
}