[Luogu P3327] [SDOI2015]约数个数和

题面:

传送门:洛咕


Solution

首先,我们需要一个结论:
$\large d(i,j)=\sum_{x|i}\sum_{y|j}[gcd(x,y)=1]$

证明
理性证明请看这篇博客的例五
本蒟蒻提供一个感性证明的方法:如果$x*y$是$i*j$的因数,我们必须有$x|i,y|j$,而后面那个$gcd(x,y)$是用来去重的

有了这个柿子之后,我们之后的推导就比较套路了:
为了方便讨论,之后的柿子均默认$m>n$
为了方便书写,之后的除法默认向下取整
原式:
$\large \sum_{i=1}^n\sum_{j=1}^md(i*j)$

把我们上面的结论代进去
$\large \sum_{i=1}^n\sum_{j=1}^m\sum_{x|i}\sum_{y|j}[gcd(x,y)=1]$

根据套路,这里的$x|i$与$y|j$应该写成枚举的形式:
$\large \sum_{i=1}^n\sum_{j=1}^m\sum_{x=1}^{n}[x|i]\sum_{y=1}^{m}[y|j][gcd(x,y)=1]$

这里显然可以把$x,y$的和式写到最前面去:
$\large \sum_{x=1}^{n}[x|i]\sum_{y=1}^{m}[y|j]\sum_{i=1}^n\sum_{j=1}^m[gcd(x,y)=1]$
然后就可以去掉$x,y$后面的那两个判断式啦
$\large \sum_{x=1}^{n}\sum_{y=1}^{m}\sum_{i=1}^{n/x}\sum_{j=1}^{m/y}[gcd(x,y)=1]$

然后我们再去套路后面那个$gcd$,根据莫比乌斯函数的性质$[x=1]=\sum_{d|x}\mu(d)$,我们就把$gcd$带入得
$\large \sum_{x=1}^{n}\sum_{y=1}^{m}\sum_{i=1}^{n/x}\sum_{j=1}^{m/y}\sum_{d|gcd(x,y)}\mu(d)$
然后去枚举$d$
$\large \sum_{x=1}^{n}\sum_{y=1}^{m}\sum_{i=1}^{n/x}\sum_{j=1}^{m/y}\sum_{d=1}^{n}\mu(d)[d|gcd(x,y)]$
套路地把$\mu$的和式丢到最前面,化简一下就有:
$\large \sum_{d=1}^{n}\mu(d)\sum_{x=1}^{n/d}\sum_{y=1}^{m/d}\sum_{i=1}^{n/(x*d)}\sum_{j=1}^{m/(y*d)}1$
然后有
$\large \sum_{d=1}^{n}\mu(d)\sum_{x=1}^{n/d}\sum_{y=1}^{m/d}\frac{n}{x*d}\frac{m}{y*d}$

移项整理一下:
$\large \sum_{d=1}^{n}\mu(d)\sum_{x=1}^{n/d}\frac{n}{x*d}\sum_{y=1}^{m/d}\frac{m}{y*d}$

好了,到这里我就不会推了
kAl7Ct.jpg

之后的内容感谢@Maxwei_wzj的教学
事实上,这个柿子我们已经可以算了。
$\large \sum_{d=1}^{n}\mu(d)\sum_{x=1}^{n/d}\frac{n}{x*d}\sum_{y=1}^{m/d}\frac{m}{y*d}$
首先,我们有一个结论(我不会证)
$\large x/(y*z)=(x/y)/z$ (这里的除法向下取整)

我们的柿子就可以变为
$\large \sum_{d=1}^{n}\mu(d)\sum_{x=1}^{n/d}\frac{\frac{n}{d}}{x}\sum_{y=1}^{m/d}\frac{\frac{m}{d}}{y}$

然后我们再设一个方程:
$f(x)=\sum_{i=1}^x\frac{x}{i}$
我们柿子就可以写为:
$\large \sum_{d=1}^{n}\mu(d)f(\frac{n}{d})f(\frac{m}{d})$

诶?我们好像又能整除分块了。
没错,我们只需要先用$O(n*\sqrt n)$的整除分块预处理出$f$
然后再每次$O(\sqrt n)$整除分块算出那个柿子就好。

时间复杂度$O(n*\sqrt n)$
完结撒花✿✿ヽ(°▽°)ノ✿


Code

//Luogu  P3327 [SDOI2015]约数个数和
//Jan,22ed,2019
//莫比乌斯反演
#include<iostream>
#include<cstdio>
using namespace std;
long long read()
{
    long long x=0,f=1; char c=getchar();
    while(!isdigit(c)){if(c=='-') f=-1;c=getchar();}
    while(isdigit(c)){x=x\*10+c-'0';c=getchar();}
    return x\*f;
}
const int N=50000+1000;
const int M=50000;
int prime[N],p_cnt,mu[N];
bool noPrime[N];
void GetPrime(int n)
{
    noPrime[1]=true,mu[1]=1;
    for(int i=2;i<=n;i++)
    {
        if(noPrime[i]==false)
            prime[++p_cnt]=i,mu[i]=-1;
        for(int j=1;j<=p_cnt and i\*prime[j]<=n;j++)
        {
            noPrime[i\*prime[j]]=true;
            if(i%prime[j]==0)
            {
                mu[i\*prime[j]]=0;
                break;
            }
            mu[i\*prime[j]]=mu[i]\*mu[prime[j]];
        }
    }
}
long long f[N],pre_mu[N];
int main()
{
    GetPrime(M);
    for(int i=1;i<=M;i++)
        for(int l=1,r;l<=i;l=r+1)
        {
            r=i/(i/l);
            f[i]+=(r-l+1)\*(i/l);
        }
    for(int i=1;i<=M;i++)
        pre_mu[i]=pre_mu[i-1]+mu[i];

    int T=read();
    for(;T>0;T--)
    {
        long long n=read(),m=read();
        if(n>m) swap(n,m);

        long long ans=0;
        for(int l=1,r;l<=n;l=r+1)
        {
            r=min(n/(n/l),m/(m/l));
            ans+=(pre_mu[r]-pre_mu[l-1])\*f[n/l]\*f[m/l];
        }

        printf("%lld\n",ans);
    }
    return 0;
}

点赞

发表评论

电子邮件地址不会被公开。必填项已用 * 标注