【bzoj4176】Lucas的数论 莫比乌斯反演+杜教筛

Description

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

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

求如下表达式的值：

其中 表示ij的约数个数。

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

Input

第一行一个整数n。

Output

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

Sample Input

2

Sample Output

8

HINT

对于100%的数据n <= 10^9。

Sol

这个题的最大难点，在于对(f(ij))的变形，只要这个变形正确了，后面的就是更换求和指标经典套路+数论分块经典套路了。

(f(ij)=sum_{x|i}sum_{y|j}[(x,y)=1])

证明：ij的某个因子一定是i的某个因子*j的某个因子乘起来的，我们不妨设为i和(frac{j}{y})的某个因子，那么设(p=(x,y))，那么你在x中包括了p这个因子，又在(frac{j}{y})中把它消掉，就没意义了，也就会重复统计，所以只有((x,y)=1)的时候才会有合法的贡献。

然后有了((x,y)=1)这个条件，直接上莫比乌斯函数：

(sum^{n}_{i=1}sum_{j=1}^{n}f(i,j)=sum^{n}_{i=1}sum_{j=1}^{n}sum_{x|i}sum_{y|j}sum_{d|x,d|y}mu(d))

后面一步就是更换求和指标啦，把d提到最前面，x,y其次，i,j最后面，因为有两部分完全相同，所以这个式子就变成了：

$sumlimits_{d=1}^n mu(d)(sumlimits_{i=1}^{lfloor {nover i} floor} {lfloor {nover i*d} floor})^2 $

先对于n/i分块，然后对于n/(i/d)分块，前面的莫比乌斯函数根据n/i的分块范围在线使用杜教筛计算。

时间复杂度$O(n^{3/4}+n{2/3}logn) $

Code

#include <bits/stdc++.h>
#define ll long long
using namespace std;
map<ll,ll>mmp;int sum[1000005],pri[1000005],vis[1000005],tot,mu[1000005],n,ls,ans,P=1000000007;
ll djs(int x)
{
    if(x<=1e6) return sum[x];
    if(mmp.count(x)) return mmp[x];
    int ans=1,ls;
    for(int i=2;i<=x;i=ls+1) ls=x/(x/i),ans=(ans-1ll*(ls-i+1)*djs(x/i)%P+P)%P;
    return mmp[x]=ans;
}
ll cal(ll x)
{
    ll ans=0,ls;
    for(int i=1;i<=x;i=ls+1) ls=x/(x/i),ans=(ans+1ll*(ls-i+1)*(x/i)%P)%P;
    return 1ll*ans*ans%P;
}
int main()
{
    mu[1]=sum[1]=1;
    for(int i=2;i<=1000000;i++)
    {
        if(!vis[i]) pri[++tot]=i,mu[i]=-1;
        for(int j=1;j<=tot&&i*pri[j]<=1000000;j++)
        {
            vis[i*pri[j]]=1;
            if(i%pri[j]==0){mu[i*pri[j]]=0;break;}
            mu[i*pri[j]]=-mu[i];
        }
        sum[i]=(sum[i-1]+mu[i]+P)%P;
    }
    scanf("%d",&n);
    for(int i=1;i<=n;i=ls+1) ls=n/(n/i),ans=(ans+1ll*(djs(ls)-djs(i-1)+P)%P*cal(n/i)%P)%P;
    printf("%d
",ans);
}