bzoj 4176 Lucas的数论

Posted 2021-02-15 jklover

bzoj 4176 Lucas的数论

• 和约数个数和那题差不多.只不过那个题是多组询问,这题只询问一次,并且 (n) 开到了 (10^9).
  [ egin{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 floor} sum_{y=1}^{lfloor frac M d floor}lfloor frac {N}{dx} floor lfloor frac {N}{dy} floor\\&= sum_{d=1}^N mu(d)cdot sum_{x=1}^{lfloor frac N d floor}lfloor frac {N}{dx} floorcdot sum_{y=1}^{lfloor frac N d floor}lfloor frac {N}{dy} floor. end{align*} ]

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

• 则答案为

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

• (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;
}