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

题目链接
题意：求∑ni=1∑nj=1f(ij)$\sum _{i=1}^n\sum _{j=1}^{n}f(ij)$模1000000007$1000000007$的值。f(n)$f(n)$为n的约数个数。
题解：我们有一个结论：f(nm)=∑i|n∑j|m(gcd(i,j)==1)$f(nm)=\sum _{i|n}\sum _{j|m}(gcd(i,j)==1)$
这是为什么呢？
考虑到nm$nm$的任何一个约数都可以表示成满足i|n$i|n$，j|m$j|m$的i∗mj$i\ast \frac{m}{j}$的形式，上面的式子就是枚举i，j的值的过程。那为什么要求gcd(i,j)==1$gcd(i,j)==1$呢？如果有d满足d|i，d|j$d|i，d|j$那么i∗mj=ip∗mjp$i\ast \frac{m}{j}=\frac{i}{p}\ast \frac{m}{\frac{j}{p}}$，重复计算了。所以我们有了上面的式子。
ans$ans$
=∑ni=1∑nj=1f(ij)$=\sum _{i=1}^n\sum _{j=1}^{n}f(ij)$
=∑ni=1∑nj=1∑a|i∑b|j(gcd(a,b)==1)$=\sum _{i=1}^n\sum _{j=1}^n\sum _{a|i}\sum _{b|j}(gcd(a,b)==1)$
=∑ni=1∑nj=1∑a|i∑b|j∑k|gcd(a,b)μ(k)$=\sum _{i=1}^n\sum _{j=1}^n\sum _{a|i}\sum _{b|j}\sum _{k|gcd(a,b)}\mu (k)$
=∑ni=1∑nj=1∑a|i∑b|j∑k|a,k|bμ(k)$=\sum _{i=1}^n\sum _{j=1}^n\sum _{a|i}\sum _{b|j}\sum _{k|a,k|b}\mu (k)$
=∑nk=1μ(k)∑⌊nk⌋a=1∑⌊nk⌋b=1∑⌊nak⌋i=1∑⌊nbk⌋j=1$=\sum _{k=1}^n\mu (k)\sum _{a=1}^{\lfloor \frac{n}{k}\rfloor }\sum _{b=1}^{\lfloor \frac{n}{k}\rfloor }\sum _{i=1}^{\lfloor \frac{n}{ak}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{bk}\rfloor }$
=∑nk=1μ(k)(∑⌊nk⌋i=1⌊nik⌋)2$=\sum _{k=1}^n\mu (k)(\sum _{i=1}^{\lfloor \frac{n}{k}\rfloor }\lfloor \frac{n}{ik}\rfloor {)}^2$
于是我们用杜教筛算出μ$\mu$的前缀和，再分块计算即可。
杜教筛
∑ni=1∑j|iμ(j)=1$\sum _{i=1}^n\sum _{j|i}\mu (j)=1$，因为只有i=1时有贡献1。
=>∑nj=1∑⌊nj⌋i=1μ(j)=1$\sum _{j=1}^n\sum _{i=1}^{\lfloor \frac{n}{j}\rfloor }\mu (j)=1$
=>∑ni=1∑⌊ni⌋j=1μ(j)=1$\sum _{i=1}^n\sum _{j=1}^{\lfloor \frac{n}{i}\rfloor }\mu (j)=1$
=>∑nj=1μ(j)=1−∑ni=2∑⌊ni⌋j=1μ(j)$\sum _{j=1}^n\mu (j)=1-\sum _{i=2}^n\sum _{j=1}^{\lfloor \frac{n}{i}\rfloor }\mu (j)$，即把i=1$i=1$带入。
于是我们可以通过记忆化搜索求出前缀和。如果预先线性筛出较小的前缀和，时间复杂度大概为O(n23)$O(n^{\frac{2}{3}})$。我也不会算
代码：

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
const int N=5000005,M=1000033;
const ll mod=1000000007;
int n,p[N/10];
ll ans,tmp,mu[N];
bool vis[N];
struct HashSet{
    int cnt,head[M],to[M],nxt[M];
    ll v[M];
    bool count(int f){
        int x=f%M;
        for(int i=head[x];i;i=nxt[i]){
            if(to[i]==f){
                return true;
            }
        }
        return false;
    }
    ll get(int f){
        int x=f%M;
        for(int i=head[x];i;i=nxt[i]){
            if(to[i]==f){
                return v[i];
            }
        }
    }
    void set(int f,ll val){
        int x=f%M;
        to[++cnt]=f;
        nxt[cnt]=head[x];
        v[cnt]=val;
        head[x]=cnt;
    }
}s;
ll solve(int n){
    if(n<=5000000){
        return mu[n];
    }
    if(s.count(n)){
        return s.get(n);
    }
    ll res=1;
    for(int i=2,last;i<=n;i=last+1){
        last=n/(n/i);
        res-=solve(n/i)*(last-i+1);
        res%=mod;
    }
    s.set(n,res);
    return res;
}
ll calc(int n){
    ll res=0;
    for(int i=1,last;i<=n;i=last+1){
        last=n/(n/i);
        res+=(n/i)*(last-i+1)%mod;
        res%=mod;
    }
    return res;
}
int main(){
    mu[1]=1;
    for(int i=2;i<=5000000;i++){
        if(!vis[i]){
            p[++p[0]]=i;
            mu[i]=-1;
        }
        for(int j=1;j<=p[0]&&i*p[j]<=5000000;j++){
            vis[i*p[j]]=true;
            if(i%p[j]){
                mu[i*p[j]]=-mu[i];
            }else{
                break;
            }
        }
    }
    for(int i=1;i<=5000000;i++){
        mu[i]+=mu[i-1];
        mu[i]%=mod;
    }
    scanf("%d",&n);
    for(int i=1,last;i<=n;i=last+1){
        last=n/(n/i);
        tmp=calc(n/i);
        ans+=tmp*tmp%mod*(solve(last)-solve(i-1))%mod;
        ans%=mod;
    }
    printf("%lld
",(ans+mod)%mod);
    return 0;
}