bzoj 4176 Lucas的数论

• 和约数个数和那题差不多.只不过那个题是多组询问,这题只询问一次,并且 \(n\) 开到了 \(10^9\).
  \[ \begin{align*} \sum_{i=1}^N \sum_{j=1}^N f(ij)&= \sum_{i=1}^N \sum_{j=1}^N \sum_{x|i} \sum_{y|j}[gcd(x,y)=1]\\&= \sum_{i=1}^N \sum_{j=1}^N \sum_{x|i} \sum_{y|j} \sum_{d|gcd(x,y)}\mu(d)\\&= \sum_{d=1}^N \mu(d)\sum_{x=1}^{\lfloor \frac N d \rfloor} \sum_{y=1}^{\lfloor \frac M d \rfloor}\lfloor \frac {N}{dx} \rfloor \lfloor \frac {N}{dy} \rfloor\\&= \sum_{d=1}^N \mu(d)\cdot \sum_{x=1}^{\lfloor \frac N d \rfloor}\lfloor \frac {N}{dx} \rfloor\cdot \sum_{y=1}^{\lfloor \frac N d \rfloor}\lfloor \frac {N}{dy} \rfloor. \end{align*} \]

• 记 \(f'(n)=\sum_{i=1}^n \lfloor \frac n i \rfloor\).

• 则答案为

\[ \sum_{d=1}^N \mu(d) \cdot f'(\lfloor\frac N d\rfloor)\cdot f'(\lfloor \frac N d \rfloor). \]

• \(N\) 是 \(10^9\) 级别,所以用杜教筛求 \(\mu\) 的前缀和.然后套两个整除分块,外层算答案,里层算 \(f'\) 即可.

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
inline int read()
{
    int out=0,fh=1;
    char jp=getchar();
    while ((jp>'9'||jp<'0')&&jp!='-')
        jp=getchar();
    if (jp=='-')
        fh=-1,jp=getchar();
    while (jp>='0'&&jp<='9')
        out=out*10+jp-'0',jp=getchar();
    return out*fh;
}
const int P=1e9+7;
const int inv2=(P+1)>>1;
inline int add(int a,int b)
{
    return (a + b) % P;
}
inline int mul(int a,int b)
{
    return 1LL * a * b % P;
}
inline int sub(int a,int b)
{
    return add(a,P-b);
}
const int MAXN=3e6+10;
int n,ans=0;
int f[MAXN],prime[MAXN],cnt=0,mu[MAXN],ism[MAXN],summu[MAXN];
int calc_F(int i)
{
    int res = 0;
    for(int l=1,r; l<=i; l=r+1)
    {
        r = i/(i/l);
        res = add(res,mul(r-l+1,(i/l)));
    }
    return res;
}
void init(int N)
{
    ism[1] = 1;
    mu[1] = 1;
    for(int i=2; i<=N; ++i)
    {
        if(!ism[i])
        {
            prime[++cnt] = i;
            mu[i] = -1;
        }
        for(int j=1; j<=cnt; ++j)
        {
            ll num = i * prime[j];
            if(num > N)
                break;
            ism[num] = 1;
            if(i % prime[j] == 0)
                break;
            else
                mu[num] = -mu[i];
        }
    }
    for(int i=1; i<=N; ++i)
        summu[i] = add(summu[i-1],P+mu[i]);
}
int AP(int x)
{
    return mul(mul(x,x+1),inv2);
}
map<int,int> mp;
const int sqN=31200;
int calc(int x)
{
    if(x<=sqN)
        return summu[x];
    if(mp.find(x)!=mp.end())
        return mp[x];
    int res=1;
    for(int l=2,r;l<=x;l=r+1)
    {
        r=x/(x/l);
        int tmp=mul(r-l+1,calc(x/l));
        res=add(res,P-tmp);
    }
    return mp[x]=res;
}
void solve()
{
    init(sqN);
    for(int l=1,r;l<=n;l=r+1)
    {
        r=n/(n/l);
        int tmp=add(calc(r),P-calc(l-1));
        tmp=mul(tmp,mul(calc_F(n/l),calc_F(n/l)));
        ans=add(tmp,ans);
    }
    cout<<ans<<endl;
}
int main()
{
    freopen("math.in","r",stdin);
    freopen("math.out","w",stdout);
    n=read();
    solve();
    return 0;
}