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。

数学问题 莫比乌斯反演 脑洞题

$ans(n)=sum_{i=1}^{n} sum_{j=1}^{n} f(i,j)$
$=sum_{i=1}^{n} sum_{j=1}^{n} sum_{k=1}^{n^2} [frac{k}{gcd(i,j)}|j]$
(枚举gcd)
$=sum_{d=1}^{n} d sum_{i=1}^{lfloor frac{n}{d} floor}sum_{j=1}^{lfloor frac{n}{d} floor}sum_{k=1}^{lfloor frac{n^2}{d} floor} [k|j][gcd(k,i)=1] $
(消掉j)
$=sum_{d=1}^{n} d sum_{i=1}^{lfloor frac{n}{d} floor}sum_{k=1}^{lfloor frac{n^2}{d} floor} d lfloor frac{n}{dk} floor [gcd(k,i)=1] $
(约掉一个d）
$=sum_{d=1}^{n} sum_{i=1}^{lfloor frac{n}{d} floor}sum_{k=1}^{lfloor frac{n^2}{d} floor} lfloor frac{n}{k} floor [gcd(k,i)=1] $

$=sum_{d=1}^{n} sum_{i=1}^{lfloor frac{n}{d} floor}sum_{k=1}^{lfloor frac{n^2}{d} floor} lfloor frac{n}{k} floor sum_{t|(ik)} mu(t)$

$= sum_{t=1}^{1} mu(t) sum_{d=1}^{n} sum_{i=1}^{lfloor frac{n}{dt} floor}sum_{k=1}^{lfloor frac{n^2}{dt} floor} lfloor frac{n}{kt} floor$
(注意到i对后面没有影响，把求和变为乘)
$= sum_{t=1}^{1} mu(t) sum_{d=1}^{n} lfloor frac{n}{dt} floorsum_{k=1}^{lfloor frac{n^2}{dt} floor} lfloor frac{n}{kt} floor$
(注意到d在大于$ lfloor frac{n}{t} floor$ 时，后面全是0，所以缩小求和上界，后面的k同理)
$= sum_{t=1}^{1} mu(t) sum_{d=1}^{lfloor frac{n}{t} floor} lfloor frac{n}{dt} floorsum_{k=1}^{lfloor frac{n}{t} floor} lfloor frac{n}{kt} floor$

$= sum_{t=1}^{1} mu(t) (sum_{d=1}^{lfloor frac{n}{t} floor} lfloor frac{n}{dt} floor)^2 $

$mu$的前缀和可以用bzoj3944的方法筛出来，后面的平方项可以分块计算。
也就是说分块套分块就可以愉快地出解了。

隔壁Sfailsth告诉我可以直接用Bzoj3994的结论来做
↓也就是从这个式子开始：
$sum_{i=1}^{N}sum_{j=1}^{M}d(ij)=sum_{i=1}^{N}sum_{j=1}^{M}lfloorfrac{N}{i} floorlfloorfrac{M}{j} floorlbrack gcd(i,j)=1 brack $
确实很方便

（做完这道题之后，陪伴了我们一年的权限号过期了）

#include<iostream>
#include<algorithm>
#include<cstring>
#include<cstdio>
#include<cmath>
#include<map>
#define LL long long
using namespace std;
const int mod=1e9+7;
const int mxn=5000010;
int read(){
    int x=0,f=1;char ch=getchar();
    while(ch<'0' || ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0' && ch<='9'){x=x*10+ch-'0';ch=getchar();}
    return x*f;
}
int pri[mxn],mu[mxn],cnt=0;
bool vis[mxn];
void init(){
    mu[1]=1;
    for(int i=2;i<mxn;i++){
        if(!vis[i]){
            pri[++cnt]=i;
            mu[i]=-1;
        }
        for(int j=1;j<=cnt && (LL)pri[j]*i<mxn;j++){
            int tmp=pri[j]*i;
            vis[tmp]=1;
            if(i%pri[j]==0){vis[tmp]=1;mu[tmp]=0;break;}
            mu[tmp]=-mu[i];
        }
    }
    for(int i=2;i<mxn;i++){mu[i]=mu[i-1]+mu[i];}
    return;
}
map<int,LL>mpmu;
int F(int x){
    if(x<mxn)return mu[x];
    if(mpmu.count(x))return mpmu[x];
    int res=1;
    for(int i=2,pos;i<=x;i=pos+1){
        int y=x/i;
        pos=x/y;
        res=((LL)res-(LL)(pos-i+1)*F(y)%mod+mod)%mod;
    }
    mpmu[x]=res;
    return res;
}
int calc(int x){
    int res=0;
    for(int i=1,pos;i<=x;i=pos+1){
        int y=x/i;
        pos=x/y;
        res=(res+(pos-i+1)*(LL)y%mod)%mod;
    }
    return (LL)res*res%mod;
}
int main(){
    init();
    int n=read();
    LL ans=0;
    for(int i=1,pos;i<=n;i=pos+1){
        int y=n/i;pos=n/y;
        ans=((LL)ans+(((LL)F(pos)-F(i-1)%mod))*calc(y)%mod)%mod;
    }
    ans=(ans%mod+mod)%mod;
    printf("%lld
",ans);
    return 0;
}