bzoj 4176 Lucas的数论 题解【莫比乌斯反演】【杜教筛】

点击量：219

怎么又有我不知道的东西啊..自闭了

目录

Description

去年的 Lucas 非常喜欢数论题，但是一年以后的 Lucas 却不那么喜欢了。

在整理以前的试题时，发现了这样一道题目“求 $\sum f(i)$，其中 $1\le i\le n$”，其中 $f(i)$ 表示i的约数个数。他现在长大了，题目也变难了。

求如下表达式的值：
$$
\sum_{i=1}^n\sum_{j=1}^nf(ij)
$$
其中 $f(ij)$ 表示 $ij$ 的约数个数。

他发现答案有点大，只需要输出模 $1000000007$ 的值。

Input

第一行一个整数 $n$。

Output

一行一个整数 $ans$，表示答案模 $1000000007$ 的值。

Sample Input

2

Sample Output

8

HINT

对于 $100\%$ 的数据 $n\le 10^9$。

题解：

这个题是 SDOI2015 约数个数和（题面 题解） 的升级版。

首先化式子，根据上一道题推出的结论
$$
\begin{aligned}
&\sum_{i=1}^n\sum_{j=1}^nd(ij)\\
=&\sum_{i=1}^n\sum_{j=1}^n\sum_{x|i}\sum_{y|j}[\gcd(x,y)=1]
\end{aligned}
$$
（题面中的 $f$ 已转成 $d$）令
$$
f(k)=\sum_{i=1}^n\sum_{j=1}^n\sum_{x|i}\sum_{y|j}[\gcd(x,y)=k]\\
\begin{aligned}
g(k)&=\sum_{k|d}f(d)\\
&=\sum_{i=1}^n\sum_{j=1}^n\sum_{x|i}\sum_{y|j}[k|\gcd(x,y)]\\
&=\sum_{i=1}^\left\lfloor\frac nk\right\rfloor\sum_{j=1}^\left\lfloor\frac nk\right\rfloor\sum_{x|i}\sum_{y|j}[1|\gcd(x,y)]\\
&=\sum_{i=1}^\left\lfloor\frac nk\right\rfloor\sum_{j=1}^\left\lfloor\frac nk\right\rfloor d(i)d(j)
\end{aligned}
\\
\begin{aligned}
f(k)&=\sum_{k|d}\mu\left(\frac dk\right)g(d)\\
&=\sum_{i=1}^\left\lfloor\frac nk\right\rfloor\mu(i)g(ki)
\end{aligned}
$$
定义 $sumd(k)=\sum_{i=1}^kd(i)$，则 $g(k)=sumd^2\left(\left\lfloor \frac nk\right\rfloor\right)$。因为我们要求的答案是 $f(1)$，所以代入 $k=1$，答案是
$$
\sum_{i=1}^n\mu(i)sumd^2\left(\left\lfloor\frac ni\right\rfloor\right)
$$
现在我们可以 $O(n)$ 求它了。

但是这还不够。看上去可以数论分块把枚举的复杂度降低到 $O(\sqrt n)$，$\mu(i)$ 可以杜教筛。但是我们不知道计算 $sumd$ 的复杂度。

//以下部分感谢 _rqy 的指导

根据定义，我们把 $sumd$ 拆开：
$$
\begin{aligned}
sumd(k)&=\sum_{i=1}^k\sum_{j|i}1\\
&=\sum_{j=1}^k\sum_{j|i}1\\
&=\sum_{j=1}^k\left\lfloor\frac kj\right\rfloor
\end{aligned}
$$
因此计算 $sumd(k)$ 的复杂度是 $O(\sqrt k)$ 的，如果利用杜教筛的预处理+记忆化，计算 $\forall 1\le i\le n,sumd(i)$ 的时间复杂度是 $O\left(n^\frac 23\right)$。

因此总时间复杂度为 $O\left(n^\frac 23\right)$。

在这里留下一个杜教筛的查错备忘，因为小数据（能手算的）会直接通过形如 if(x<=3000000) 这种判断跑过去，因此线筛写对了就可以过样例。样例一定要把 if 删了多跑几遍，而且这种题对拍应该挺方便。

Code：

#include<cstdio>
#include<cstring>
#define ll long long
#define p 1000000007
int mu[3000100],gmu[5000];
ll d[3000100],f[3000100],gd[5000];
bool is[3000100],usedmu[5000],usedd[5000];
int pri[3000100],cnt=0,n;
int calcmu(int x)
{
    if(x<=3000000)
        return mu[x];
    int t=n/x;
    if(usedmu[t])
        return gmu[t];
    usedmu[t]=1;
    ll ans=1;
    for(int i=2;i<=x;++i)
    {
        int nxt=x/(x/i);
        ans=(ans-(ll)calcmu(x/i)*(nxt-i+1)%p+p)%p;
        i=nxt;
    }
    return gmu[t]=ans;
}
int calcd(int x)
{
    if(x<=3000000)
        return d[x];
    int t=n/x;
    if(usedd[t])
        return gd[t];
    usedd[t]=1;
    ll ans=0;
    for(int i=1;i<=x;++i)
    {
        int nxt=x/(x/i);
        ans=(ans+(ll)x/i*(nxt-i+1))%p;
        i=nxt;
    }
    return gd[t]=ans;
}
int main()
{
    is[0]=is[1]=1;
    mu[1]=1;
    d[1]=1;
    for(int i=2;i<=3000000;++i)
    {
        if(!is[i])
        {
            pri[++cnt]=i;
            mu[i]=-1;
            f[i]=2;
            d[i]=2;
        }
        for(int j=1;j<=cnt&&i*pri[j]<=3000000;++j)
        {
            is[i*pri[j]]=1;
            if(i%pri[j])
            {
                mu[i*pri[j]]=-mu[i];
                f[i*pri[j]]=2;
                d[i*pri[j]]=2*d[i];
            }
            else
            {
                mu[i*pri[j]]=0;
                f[i*pri[j]]=f[i]+1;
                d[i*pri[j]]=d[i]/f[i]*f[i*pri[j]];
                break;
            }
        }
        d[i]=(d[i-1]+d[i])%p;
        mu[i]=((mu[i-1]+mu[i])%p+p)%p;
    }
    scanf("%d",&n);
    ll ans=0;
    for(int i=1;i<=n;++i)
    {
        int nxt=n/(n/i);
        ll k=calcd(n/i);
        ans=(ans+k*k%p*(calcmu(nxt)-calcmu(i-1))%p+p)%p;
        i=nxt;
    }
    printf("%lld\n",ans);
    return 0;
}