[SDOI2018]旧试题_牛客博客

题目：

$\left(\sum_{i=1}^{A} \sum_{j=1}^{B} \sum_{k=1}^{C} d(i j k)\right) \bmod \left(10^{9}+7\right)$

其中 $d(ijk)$ 表示 $i \times j \times k$ 的约数个数。

$T \leq 10, A, B, C \leq 10^{5}, \sum \max (A, B, C) \leq 2 \times 10^{5}$

题解：

首先有 $\sum_{x | i} \sum_{y | j} \sum_{z | k}[(x, y)=1][(y, z)=1][(x, z)=1]$

具体证明我不会，见他人博客

$\sum_{i=1}^{A} \sum_{j=1}^{B} \sum_{k=1}^{C} d(i j k)=\sum_{i=1}^{A} \sum_{j=1}^{B} \sum_{k=1}^{C} \sum_{x | i} \sum_{y | j} \sum_{z | k}[(x, y)=1][(y, z)=1][(x, z)=1]$ $=\sum_{x=1}^{A} \sum_{y=1}^{B} \sum_{z=1}^{C}[(x, y)=1][(y, z)=1][(x, z)=1]\left[\frac{A}{x}\right\rfloor\left\lfloor\frac{B}{y}\right\rfloor\left\lfloor\frac{C}{z}\right\rfloor$

设 $a_{i}=\left\lfloor\frac{A}{i}\right\rfloor, b_{i}=\left\lfloor\frac{B}{i}\right\rfloor, c_{i}=\left\lfloor\frac{C}{i}\right\rfloor, f_{i}=[i=1]$

$\sum_{x=1}^{A} \sum_{y=1}^{B} \sum_{z=1}^{C} a_{x} b_{y} c_{z} f_{(x, y)} f_{(y, z)} f_{(x, z)}$

发现了什么，这个式子是不是就是 #2476. 「2018 集训队互测 Day 3」蒜头的奖杯中的

$\sum_{i=1}^{n} \sum_{j=1}^{n} \sum_{k=1}^{n} A_{i} B_{j} C_{k} D_{\operatorname{gcd}(i, j)} E_{\mathrm{gcd}(i, k)} F_{\mathrm{gcd}(j, k)}$

然后，我们来推一下这个式子，感觉这个式子还是蛮重要的。

首先:

$a_x=\sum_{i=1}^{n}a_i[i=x]$ $=\sum_{i | x} a_{i}\left[\frac{x}{i}=1\right]$

常规莫比乌斯反演, $[x=1]= \sum_{p | x} \mu(p)$

$=\sum_{i | x} a_{i} \sum_{p | \frac{x}{i}} \mu(p)$

令 $p=pi$

$=\sum_{p | x} \sum_{i | p} a_{i} \mu\left(\frac{p}{i}\right)$

这一步我不是很看得懂，手动模拟后发现就是 $a$ 和 $\mu$ 的卷积后的 $p$ 位之和

$=\sum_{p | x}(a * \mu)(p)$

然后我们设 $f(a)=a * \mu$ ，原式就变成：

$\sum_{k=1}^{n} C_{k} \sum_{i=1}^{n} \sum_{j=1}^{n} A_{i} B_{j} D_{(i, j)} \sum_{p |(i, k)} f(E)(p) \sum_{q |(j, k)} f(F)(q)$

看那几个 $A_i,B_j,C_k$ 不爽（因为莫比乌斯反演喜欢 $\sum_{x | i}$ 的形式）,设 $g(a)(x)=\sum_{x | i} a_{i}$

枚举 $p,q$ 然后 $k={lcm(p,q) | i}$ ，所以 $\sum_{k=1}^{n} C_{k}=g(C)(\operatorname{lcm}(p, q))$

因为 $i$ 一定是 $p$ 的倍数， $j$ 一定是 $q$ 的倍数，我们再把 $i,j$ 换成 $i/p,j/q$ ，即枚举 $i,j$ ，A,B现在的坐标为 $ip,q$

$\sum_{\operatorname{lcm}(p, q) \leq n} f(E)(p) \times f(F)(q) \times g(C)(\operatorname{lcm}(p, q)) \sum_{i p \leq n} A_{i p} \sum_{j q \leq n} B_{j q} D_{(i p, j q)}$

我们在把 $lcm$ 换成 $gcd$ ， $d=gcd,x=q/d,y=p/d$

$\sum_{d=1}^{n}\sum_{xy \leq\left\lfloor\frac{n}{d}\right\rfloor,(x,y)=1} f(E)(xd) \times f(F)(yd) \times g(C)(xyd) \sum_{i x \leq \left\lfloor\frac{n}{d}\right\rfloor} \sum_{j y \leq \left\lfloor\frac{n}{d}\right\rfloor} A_{ixd} B_{jyd} D_{(ixd, jyd)}$

这是当 $d$ 是已知的话（即枚举 $d$ ），

设 $P_{i}=f(E)(i d), Q_{i}=f(F)(i d), R_{i}=g(C)(i d), S_{i}=A_{i d}, T_{i}=B_{i d}, W_{i}=D_{i d},m=\left\lfloor\frac{n}{d}\right\rfloor$

$\sum_{d=1}^{n} \sum_{x y \leq m,(x, y)=1} P_{x} Q_{y} R_{x y} \sum_{i x \leq m} \sum_{j y \leq m} S_{i x} T_{j y} W_{(i x, j y)}$

如果 $x$ 与 $y$ 已知（同样枚举）那么 $P_{x} Q_{y} R_{x y}$ 就得到了，然后考虑求后面的一部分

$\sum_{i x \leq m} \sum_{j y \leq m} S_{i x} T_{j y} W_{(i x, j y)}$

这 $W_{(i x, j y)}$ 很难受啊，改一下用上面推出 $f$ 的方法，推出 $W_{(i x, j y)}=\sum_{r|(i x, j y)} f(W)(r)$

$=\sum_{i x \leq m} S_{i x} \sum_{j y \leq m} T_{j y} \sum_{r|(i x, j y)} f(W)(r)$

这步应该懂吧， $jy$ 知道了， $r$ 也知道了，然后 $ix$ 就知道了

$=\sum_{j y \leq m} T_{j y} \sum_{r | j y} f(W)(r) \sum_{r | i x} S_{i x}$

设 $u=ix,v=jy$

$\sum_{y | v} T_{v} \sum_{r | v} f(W)(r) \sum_{r | u}[x | u] S_{u}$

经过wzp巨佬的证明，时间复杂度: $O\left(n^{1.5} \log \log n\right)$

代码：

#include
using namespace std;
#define next Next
/*char buf[1<<21],*p1=buf,*p2=buf;
inline int gc(){return p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++;}*/
#define gc getchar
inline int read()
{
    int ret=0,f=0;char c=gc();
    while(!isdigit(c)){if(c=='-')f=1;c=gc();}
    while(isdigit(c)){ret=ret*10+c-48;c=gc();}
    if(f)return -ret;return ret;
}
const int mod=1e9+7;
const int N=1e5+5;
int a[N],b[N],c[N],d[N],e[N],f[N],p[N],q[N],r[N],s[N],t[N],w[N],u[N],v[N];
int n,prim[N],cnt;
bool vis[N];
inline void add(int &a,int b)
{
    a+=b;
    if(a>=mod)a-=mod;
}
inline void jian(int &a,int b)
{
    a-=b;
    if(a<0)a+=mod;
}
inline void F(int*a,int n)//f(a)=a*u
{
    for(int i=1;i<=n;i++)
        for(int j=i+i;j<=n;j+=i)jian(a[j],a[i]);
}
inline void G(int*a,int n)//g(a)(x)=Σai[x|i] 
{
    for(int i=1;i<=cnt&&prim[i]<=n;i++)
        for(int j=n/prim[i];j;--j)
            add(a[j],a[j*prim[i]]);
}
inline void H(int*a,int n)
{
    for(int i=1;i<=cnt&&prim[i]<=n;i++)
        for(int j=1;j*prim[i]<=n;j++)
            add(a[j*prim[i]],a[j]);
}
inline void work(int n)
{
    for(int i=2;i<=n;i++)
    {
        if(!vis[i])prim[++cnt]=i;
        for(int j=1;j<=cnt;j++)
        {
            if(prim[j]*i>n)break;
            vis[prim[j]*i]=1;
            if(i%prim[j]==0)break;
        }
    }
}
int main()
{
    work(100000);
    int T=read();
    while(T--)
    {
        int A=read(),B=read(),C=read();
        n=max(A,max(B,C));
        for(int i=1;i<=n;i++)
        {
            a[i]=A/i;
            b[i]=B/i;
            c[i]=C/i;
            d[i]=e[i]=f[i]=(i==1);
        }
        F(e,n);
        F(f,n);
        G(c,n);
        int ans=0;
        for(int i=1;i<=n;i++)//枚举gcd 
        {
            int m=n/i;//[n/d]
            for(int j=1;j<=m;j++)
            {
                p[j]=e[j*i];//pj=f(E)(jd)
                q[j]=f[j*i];//qj=f(F)(jd)
                r[j]=c[j*i];//ri=g(C)(id)
                s[j]=a[j*i];//si=Aid
                t[j]=b[j*i];//ti=Bid
                w[j]=d[j*i];//wi=Did
            }
            F(w,m);
            for(int x=1;x*x<=m;x++)//枚举x 
            {
                fill(u+1,u+1+m,0);
                 fill(v+1,v+1+m,0);
                for(int j=x;j<=m;j+=x)
                {
                    u[j]=s[j];
                    v[j]=t[j];
                }
                G(u,m);
                G(v,m);
                for(int j=1;j<=m;j++)
                {
                    u[j]=1ll*u[j]*w[j]%mod;
                    v[j]=1ll*v[j]*w[j]%mod;
                }
                H(u,m);
                H(v,m);
                for(int j=1;j<=m;j++)
                {
                    u[j]=1ll*u[j]*t[j]%mod;
                    v[j]=1ll*v[j]*s[j]%mod;
                }
                for(int y=x;x*y<=m;y++)//枚举y 
                {
                    if(__gcd(x,y)==1)
                    {
                        int s1=0,s2=0;
                        for(int j=y;j<=m;j+=y)
                        {
                            add(s1,u[j]);//Tjy[jy<=m]
                            add(s2,v[j]);//Six[ix<=m];
                        }
                        add(ans,1ll*s1*p[x]%mod*q[y]%mod*r[x*y]%mod);
                        if(x!=y)add(ans,1ll*s2*p[y]%mod*q[x]%mod*r[x*y]%mod);
                    }
                }
            }
        }
        printf("%d\n",ans);
    }
    return 0;
}