BZOJ-4196 Lucas的数论(莫比乌斯反演+杜教筛)

题目描述

  计算：

\[\sum_{i=1}^{n}\sum_{j=1}^{n}\sigma_0(ij) \]

  其中 \(\sigma_0(x)\) 表示 \(x\) 的约数个数，\(n\leq 10^9\)，答案对 \(10^9+7\) 取模。

分析

  首先有一个结论：

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

  证明：

  在等号左边，对于 \(ij\) 的每个质因子 \(p\)，设 \(i=i'\times p^{k_1},j=j'\times p^{k_2}\)（其中 \(k_1,k_2\) 可以为 \(0\)），则 \(ij=i'\times j'\times p^{k_1+k_2}\)，因此 \(p\) 对 \(d(ij)\) 的贡献为 \(k_1+k_2+1\)；在等号右边，设 \(x=x'\times p^{k_x},y=y'\times p^{k_y}\)，为了满足 \(\gcd(x,y)=1\)，有 \(\gcd(p^{k_x},p^{k_y})=1\)，要么 \(k_x=0,k_y\in [0,k_2]\)，共 \(k_2+1\) 种情况，要么 \(k_y=0,k_x\in [0,k_1]\)，共 \(k_1+1\) 种情况，其中需要去掉 \(k_x=0,k_y=0\) 重复的一种情况，即一共 $k_1+k_2+1 $ 种情况，与等号左边相同。

  因此：

\[\begin{aligned}&\sum_{i=1}^{n}\sum_{j=1}^{m}\sigma_0(ij)\\ =&\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{x\mid i}\sum_{y\mid j}[\gcd(x,y)=1]\\ =&\sum_{x=1}^{n}\sum_{y=1}^{m}\Big\lfloor\frac{n}{x}\Big\rfloor\Big\lfloor\frac{m}{y}\Big\rfloor[\gcd(x,y)=1]\\ =&\sum_{i=1}^{n}\sum_{j=1}^{m}\Big\lfloor\frac{n}{i}\Big\rfloor\Big\lfloor\frac{m}{j}\Big\rfloor[\gcd(i,j)=1]\\ =&\sum_{i=1}^{n}\sum_{j=1}^{m}\Big\lfloor\frac{n}{i}\Big\rfloor\Big\lfloor\frac{m}{j}\Big\rfloor\sum_{d\mid \gcd(i,j)}\mu(d)\\ =&\sum_{d=1}^{\min(n,m)}\mu(d)\sum_{i=1}^{n}\sum_{j=1}^{m}[d\mid \gcd(i,j)]\Big\lfloor\frac{n}{i}\Big\rfloor\Big\lfloor\frac{m}{j}\Big\rfloor\\ =&\sum_{d=1}^{\min(n,m)}\mu(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}\Big\lfloor\frac{n}{di}\Big\rfloor\Big\lfloor\frac{m}{dj}\Big\rfloor\\ =&\sum_{d=1}^{\min(n,m)}\mu(d)\Big(\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\Big\lfloor\frac{n}{di}\Big\rfloor\Big)\Big(\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}\Big\lfloor\frac{m}{dj}\Big\rfloor\Big) \end{aligned} \]

  \(\mu\) 的前缀和用杜教筛解决，后面两部分数论分块处理，然后分块套分块。

代码

#include<bits/stdc++.h>
using namespace std;
const int N=2e6+10;//n^(2/3)
const int mod=1e9+7;
long long prime[N+10],vis[N+10],mu[N+10],sum[N+10],cnt;
void init()
{
    mu[1]=1;
    for(int i=2;i<=N;i++)
    {
        if(!vis[i])
        {
            prime[++cnt]=i;
            mu[i]=-1;
        }
        for(int j=1;j<=cnt&&i*prime[j]<=N;j++)
        {
            vis[i*prime[j]]=1;
            if(i%prime[j]==0)
            {
                mu[i*prime[j]]=0;
                break;
            }
            else
                mu[i*prime[j]]=-mu[i];
        }
    }
    for(int i=1;i<=N;i++)
        sum[i]=sum[i-1]+mu[i];
}
map<long long,long long> mp;
long long S_mu(long long n)
{
    if(n<N)
        return sum[n];
    if(mp[n])
        return mp[n];
    long long ans=1;
    for(long long l=2,r;l<=n;l=r+1)
    {
        r=n/(n/l);
        ans=((ans-(r-l+1)*S_mu(n/l)%mod)%mod+mod)%mod;
    }
    mp[n]=(ans+mod)%mod;
    return mp[n];
}
long long solve(long long n)
{
    long long ans=0;
    for(long long l=1,r;l<=n;l=r+1)
    {
        r=n/(n/l);
        ans=(ans+(r-l+1)*(n/l)%mod)%mod;
    }
    return (ans+mod)%mod;
}
int main()
{
    init();
    long long n;
    cin>>n;
    long long ans=0;
    for(long long l=1,r;l<=n;l=r+1)
    {
        r=n/(n/l);
        ans=(ans+((S_mu(r)-S_mu(l-1))*solve(n/l)%mod*solve(n/l)%mod)%mod)%mod;
    }
    printf("%lld\n",(ans+mod)%mod);
    return 0;
}