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