bzoj 4176: Lucas的数论 莫比乌斯反演+杜教筛

题意

设$f(i)$表示i的约数个数。给出n，求$\sum_{i=1}^n\sum_{j=1}^nf(ij)$。答案模1e9+7.
n<=10^9

分析

按套路来，有

$$f(nm)=\sum_{i|n}\sum_{j|m}[(i,j)=1]$$
把这条式子代进原式，反演一波，就可以得到 $$ans=\sum_{d=1}^n\mu(d)(\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\lfloor\frac{n}{id}\rfloor)^2$$
考虑如何求 $$\sum_{i=1}^n\lfloor\frac{n}{i}\rfloor$$
把根号下取整看成n里面i的倍数数量，那么就有 $$\sum_{i=1}^n\lfloor\frac{n}{i}\rfloor=\sum_{i=1}^nf(i)$$
当n小于阈值时预处理，大于阈值时用分块来做。

代码

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#include<map>
using namespace std;

typedef long long LL;

const int N=2000005;
const int M=2000000;
const int MOD=1000000007;

int f[N],mu[N],low[N],prime[N],tot;
bool not_prime[N];
map<int,int> w;

void get_prime(int n)
{
    f[1]=mu[1]=1;
    for (int i=2;i<=n;i++)
    {
        if (!not_prime[i]) prime[++tot]=i,mu[i]=-1,f[i]=2,low[i]=i;
        for (int j=1;j<=tot&&i*prime[j]<=n;j++)
        {
            not_prime[i*prime[j]]=1;
            if (i%prime[j]==0)
            {
                low[i*prime[j]]=low[i]*prime[j];
                f[i*prime[j]]=f[i]+f[i/low[i]];
                break;
            }
            low[i*prime[j]]=prime[j];
            f[i*prime[j]]=f[i]*2;
            mu[i*prime[j]]=-mu[i];
        }
    }
    for (int i=1;i<=n;i++) mu[i]+=mu[i]<0?MOD:0,mu[i]+=mu[i-1],mu[i]-=mu[i]>=MOD?MOD:0;
    for (int i=1;i<=n;i++) f[i]+=f[i-1],f[i]-=f[i]>=MOD?MOD:0;
}

int get_mu(int n)
{
    if (n<=M) return mu[n];
    if (w[n]) return w[n];
    int ans=1;
    for (int i=2,last;i<=n;i=last+1)
    {
        last=n/(n/i);
        ans+=MOD-(LL)(last-i+1)*get_mu(n/i)%MOD;
        ans-=ans>=MOD?MOD:0;
    }
    return ans;
}

int get_f(int n)
{
    if (n<=M) return f[n];
    int ans=0;
    for (int i=1,last;i<=n;i=last+1)
    {
        last=n/(n/i);
        ans+=(LL)(last-i+1)*(n/i)%MOD;
        ans-=ans>=MOD?MOD:0;
    }
    return ans;
}

int solve(int n)
{
    int ans=0;
    for (int i=1,last;i<=n;i=last+1)
    {
        last=n/(n/i);int w=get_f(n/i);
        ans+=(LL)(get_mu(last)+MOD-get_mu(i-1))*w%MOD*w%MOD;
        ans-=ans>=MOD?MOD:0;
    }
    return ans;
}

int main()
{
    get_prime(M);
    int n;
    scanf("%d",&n);
    printf("%d",solve(n));
    return 0;
}