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最後面，因為有兩部分完全相同，所以這個式子就變成了：

$\sum\limits_{d=1}^n \mu(d)(\sum\limits_{i=1}^{\lfloor {n\over i}\rfloor} {\lfloor {n\over i*d}\rfloor})^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\n",ans);
}
 

【bzoj4176】Lucas的數論 莫比烏斯反演+杜教篩