题解：P8570 [JRKSJ R6] 牵连的世界

· · 题解

好纯粹的推柿子题啊，主要借鉴了 @littlez_meow 的题解，写得太好了，很清晰，但作为蒟蒻，我来说一下我的理解。

进入正题

首先，这两个函数的乘积一定是不好维护的，所以对这两个函数进行变形。

\varphi(ij)

证明：\varphi(ij)=\dfrac{\varphi(i)\varphi(j)}{\varphi(\gcd(i,j))}\times\gcd(i,j)

考虑 i=p_{i,1}^{k_{i,1}}p_{i,2}^{k_{i,2}}\cdots p_{i,n}^{k_{i,n}} 和 j=p_{j,1}^{k_{j,1}}p_{j,2}^{k_{j,2}}\cdots p_{j,m}^{k_{j,m}}，p 都为质数，由定义则 \varphi(i)=i\times \frac{p_{i,1}-1}{p_{i,1}}\times \frac{p_{i,2}-1}{p_{i,2}}\times\cdots\times \frac{p_{i,n}-1}{p_{i,n}}，\varphi(j) 同理。

因为每一个 p 若同为 i 和 j 的因数，则 \varphi(i)\times\varphi(j) 会多乘上一个 \frac{p-1}{p}，所以 \varphi(i)\varphi(j) 除以 \varphi(\gcd(i,j))，即除掉了所有重复的 \frac{p-1}{p}，但是会多除掉一个 \gcd(i,j)，所以乘回来就好了。

\sigma_0(ij)

证明：\sigma_0(ij)=\sum_{x\mid i}\sum_{y\mid j}[\gcd(x,y)=1]

考虑把 i 分解为 p_{i,1}^{k_{i,1}}p_{i,2}^{k_{i,2}}\cdots p_{i,n}^{k_{i,n}}，把 j 也分解了，计算每一个因数的贡献。

若一个因数 p 不属于 i，则 j 中的 k_j 个 p 可随意取，产生了 k_j 个因数，而不属于 j 同理。若都有 p 这个因数，则能产生 k_i+k_j 个因数，数量上等于使 j 取 0 个 p，再从 i 里面取 p 的 k_i 种可能加 i 取 0 个 p，再取 j 中的 k_j 种可能，用数学方式表达就是 \sum_{x\mid p^{k1}}\sum_{y\mid p^{k2}}[\gcd(x,y)=1]，这样就能完成上面的答案。再考虑从单一的 p 转移到 i 或 j，发现可以直接用乘法原理，直接将答案乘算，得 \sum_{x\mid i}\sum_{y\mid j}[\gcd(x,y)=1] 等于右边。

证明：\sigma_0(ij)=\sum_{d\mid i,d\mid j}\mu(d)\sigma_0(\dfrac{i}{d})\sigma_0(\dfrac{j}{d})

因为已知 \sigma_0(ij)=\sum_{x\mid i}\sum_{y\mid j}[\gcd(x,y)=1]，所以经典技巧把右侧 \gcd 做莫比乌斯反演，得 \sigma_0(ij)=\sum_{x\mid i}\sum_{y\mid j}\sum_{d\mid x,d\mid y}\mu(d)，把式子改为枚举 d，得 \sum_{d\mid i,d\mid j}\mu(d)\sum_{dx\mid i}\sum_{dy\mid j}1，再用 x 替换 dx，则原式为 \sum_{d\mid i,d\mid j}\mu(d)(\sum_{x\mid\frac{i}{d}}1)(\sum_{y\mid\frac{j}{d}}1)，发现后面两个求和分别为 \sigma_0(\dfrac{i}{d}) 和 \sigma_0(\dfrac{j}{d})，得证。

好了，推完这两个难搞的东西以后，就可以套用进题目给定的式子里面了，这里保证 n<m。

\sum_{i=1}^{n}\sum_{j=1}^{m}\varphi(ij)\sigma_0(ij)

变为：

\sum_{i=1}^{n}\sum_{j=1}^{m}\dfrac{\varphi(i)\varphi(j)}{\varphi(\gcd(i,j))}\times\gcd(i,j)\sum_{r\mid i,r\mid j}\mu(r)\sigma_0(\dfrac{i}{r})\sigma_0(\dfrac{j}{r})

前面改为枚举 \gcd，得：

\sum_{d=1}^{n}\frac{d}{\varphi(d)}\sum_{i=1}^{n}\sum_{j=1}^{m}\varphi(i)\varphi(j)[\gcd(i,j)=d]\sum_{r\mid d}\mu(r)\sigma_0(\dfrac{i}{r})\sigma_0(\dfrac{j}{r})

经典套路要把 \gcd 转为莫反：

\sum_{d=1}^{n}\frac{d}{\varphi(d)}\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{m}{d}}\varphi(id)\varphi(jd)[\gcd(i,j)=1]\sum_{r\mid d}\mu(r)\sigma_0(\dfrac{id}{r})\sigma_0(\dfrac{jd}{r}) \sum_{d=1}^{n}\frac{d}{\varphi(d)}\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{m}{d}}\varphi(id)\varphi(jd)\sum_{k\mid i,k\mid j}\mu(k)\sum_{r\mid d}\mu(r)\sigma_0(\dfrac{id}{r})\sigma_0(\dfrac{jd}{r})

这个 k 很烦，继续枚举 k，

\sum_{d=1}^{n}\frac{d}{\varphi(d)}\sum_{k=1}^{\frac{n}{d}}\mu(k)\sum_{i=1}^{\frac{n}{kd}}\sum_{j=1}^{\frac{m}{kd}}\varphi(ikd)\varphi(jkd)\sum_{r\mid d}\mu(r)\sigma_0(\dfrac{ikd}{r})\sigma_0(\dfrac{jkd}{r})

此时，像合并同类项一般把两个 \varphi 放在同一位置，

\sum_{d=1}^{n}\frac{d}{\varphi(d)}\sum_{k=1}^{\frac{n}{d}}\mu(k)\sum_{r\mid d}\mu(r)[\sum_{i=1}^{\frac{n}{kd}}\sigma_0(\dfrac{ikd}{r})\varphi(ikd)][\sum_{j=1}^{\frac{m}{kd}}\sigma_0(\dfrac{jkd}{r})\varphi(jkd)]

发现后两项属于同一形式，那记 f_{a,b,c}=\sum_{i=1}^{\frac{a}{b}}\sigma_0(\dfrac{ib}{c})\varphi(ib)，则式子变为：

\sum_{d=1}^{n}\frac{d}{\varphi(d)}\sum_{k=1}^{\frac{n}{d}}\mu(k)\sum_{r\mid d}\mu(r)f_{n,kd,r}f_{m,kd,r}

好的，变得简洁多了，发现还有 r 没办法枚举，继续往前面提：

\sum_{r=1}^{n}\mu(r)r\sum_{d=1}^{\frac{n}{r}}\frac{d}{\varphi(dr)}\sum_{k=1}^{\frac{n}{dr}}\mu(k)f_{n,kdr,r}f_{m,kdr,r}

现在，f 还是没办法快速地计算，因为需要变量在两个维度，不如修改一下 f 的定义，使变量仅存在于一维之中，

f_{a,b,c}=\sum_{i=1}^{\frac{a}{bc}}\varphi(ibc)\sigma_0(ib)

此时，需要传入的数即为 f_{n,kd,r}，仅一维有变量，可以快速地计算。到此这道题就基本上完成了，只需要预处理出 \varphi，\mu，\sigma_0 就可以了，并且每个传入的数都小于 n，直接线性筛就好了。

终于好了，这是笔者第一次写这么大段的公式的题解，若有错误请指出，也请大家点点赞吧，上代码。

#include<bits/stdc++.h>
#define int long long
#define up(i,x,y) for(register int i=x;i<=y;i=-~i)
#define dn(i,y,x) for(register int i=y;i>=x;i--)
#define mst(x,y) memset(x,y,sizeof x)
#define D(x) cout<<#x<<": "<<x<<endl;
#define DE(x) cout<<#x<<": "<<x<<" ";
#define MAXSIZE 1<<21
// #define faster 1
#define pl p<<1
#define pr p<<1|1
using namespace std;
namespace lry{
    char buf[MAXSIZE],*p1=buf,*p2=buf;
    char pbuf[MAXSIZE],*pp=pbuf;
    inline int gc(){
        #if faster
            return p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++;
        #else
            return getchar();
        #endif
    }
    inline void pc(const char &c){
        #if faster
            if(pp-pbuf==MAXSIZE) fwrite(pbuf,1,MAXSIZE,stdout),pp=pbuf;*pp++=c;
        #else
            putchar(c);
        #endif
    }
    inline bool blank(const char x){return !(x^32)||!(x^10)||!(x^13)||!(x^9);}
    inline void read(double &rdx){double t=0;int x=0,s=0,f=1;char c;do c=gc();while(!isdigit(c)&&c!='-'&&c!='.');if (c=='-')f=-1,c=gc();while(isdigit(c)&&c!='.')x=(c^48)+(x<<1)+(x<<3),c=gc();if(c=='.')c=gc();else{rdx=x*f;return;}while(c>='0'&&c<='9')t=t*10+(c^48),s=-~s,c=gc();while(s--)t/=10;rdx=(x+t)*f;}
    inline void read(string &x){x="";char a=gc();for(;blank(a)&&(a^-1);a=gc());for(;!blank(a)&&(a^-1);a=gc())x+=a;}
    inline void read(char &x){char a=gc();for(;blank(a)&&(a^-1);a=gc());if(!blank(a)&&(a^-1))x=gc();}
    inline void read(int &x){int ret=0,f=0;char ch=gc();while(!isdigit(ch)){if(ch=='-')f=1;ch=gc();}while(isdigit(ch)){ret=(ret<<1)+(ret<<3)+(ch^48);ch=gc();}x=(f?-ret:ret);}
    template<typename t,typename ...T> inline void read(t &x,T&...y){read(x),read(y...);}
    inline void write(int x){if(x<0)pc('-'),x=-x;if(x>9)write(x/10);pc(x%10+48);}
    inline void write(const char *x){fputs(x,stdout);}
    inline void write(const string &x){fputs(x.c_str(),stdout);}
    inline void writeln(int x){write(x),pc('\n');}
    inline void writetr(int x){write(x),pc(' ');}
    inline void writeln(string x){write(x),pc('\n');}
    inline void writetr(string x){write(x),pc(' ');}
    template<typename t,typename ...T> inline void write(t x,T...y){write(x),write(y...);}
}
using namespace lry;
constexpr int N=3000000+10,inf=0x3f3f3f3f3f3f3f3f,mod=1000000007;
int prime[N],mu[N],phi[N],sum[N],num[N],cnt;
int tmp1[N],tmp2[N];
bool vis[N];
int n,m;
inline void init(){
    mu[1]=phi[1]=sum[1]=1;
    for(int i=2;i<N;i++){
        if(!vis[i]) prime[++cnt]=i,mu[i]=-1,phi[i]=i-1,sum[i]=2,num[i]=1;
        for(int j=1;j<=cnt&&i*prime[j]<N;j++){
            vis[i*prime[j]]=1;
            phi[i*prime[j]]=phi[i]*phi[prime[j]];
            mu[i*prime[j]]=-mu[i];
            num[i*prime[j]]=1;
            sum[i*prime[j]]=sum[i]*2;
            if(i%prime[j]==0){
                phi[i*prime[j]]=phi[i]*prime[j];
                mu[i*prime[j]]=0;
                num[i*prime[j]]=num[i]+1;
                sum[i*prime[j]]=sum[i]/(num[i]+1)*(num[i*prime[j]]+1);
                break;
            }
        }
    }
}
inline int ksm(int x,int a){
    int mul=1;
    while(a){
        if(a&1) mul=mul*x%mod;
        x=x*x%mod; a>>=1;
    }
    return mul;
}
signed main(){
    init();
    read(n,m);
    if(n>m) swap(n,m);
    int ans=0;
    for(int r=1;r<=n;r++){
        if(!mu[r]) continue;
        for(int d=1;d<=n/r;d++){
            tmp1[d]=tmp2[d]=0;
            for(int k=1;k<=n/d/r;k++) tmp1[d]=(tmp1[d]+phi[r*d*k]*sum[d*k]%mod)%mod;
            for(int k=1;k<=m/d/r;k++) tmp2[d]=(tmp2[d]+phi[r*d*k]*sum[d*k]%mod)%mod;
        }
        for(int d=1;d<=n/r;d++){
            int inv=ksm(phi[d*r],mod-2);
            for(int k=1;k<=n/d/r;k++){
                ans=(ans+mu[r]*r%mod*d%mod*inv%mod*mu[k]%mod*tmp1[k*d]%mod*tmp2[k*d]%mod+mod)%mod;
            }
        }
    }
    writeln(ans);
    #if faster
        fwrite(pbuf,1,pp-pbuf,stdout);
        pp=pbuf;
    #endif
    return 0;
}