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