链接
http://www.lydsy.com/JudgeOnline/problem.php?id=3529
题解
令
g(i)
表示
i
的所有因数之和,这个可以
f(i)
表示
x∈[1,N],y∈[1,M]
中有多少数对
(x,y)
满足
gcd(x,y)=i
。如果不考虑
≤a
的限制,那就成了求表格里所有数的和,那就要枚举最大公约数
x
,看有多少数对的
ans=∑xng(x)f(x)
反演
ans=∑x=1ng(x)∑x|dnμ(dx)⌊nd⌋⌊md⌋
这样是没法优化的,这里有一个技巧,我把它叫做 改变枚举顺序,即有原来的先枚举 x 再枚举其倍数
当 d 确定时,由于我们没有考虑
ans=∑d=1ns[d]⌊nd⌋⌊md⌋
就可以分块统计答案了,时间复杂度 O(TN−−√)
可是…有 ≤a 这个限制。可以把 g(x)>0 的数的 g(x) 暂时变成0,暴力重构前缀和,这样对时间的牺牲很大。
离线,把所有询问按照 a 排个序,所有的数按照
最终时间复杂度
代码
//莫比乌斯反演
#include <cstdio>
#include <algorithm>
#define maxn 100000
#define ll long long
#define lowbit(x) (x&-x)
#define mod ((ll)1<<31)
using namespace std;
ll prime[maxn+10], mark[maxn+10], mu[maxn+10], s[maxn+10], g[maxn+10], numg[maxn+10],
quiry[maxn][4], T, numq[maxn];
void add(ll pos, ll v){for(;pos<=maxn;pos+=lowbit(pos))s[pos]+=v;}
ll sum(ll pos)
{
ll ans=0;
for(;pos;pos-=lowbit(pos))ans+=s[pos];
return ans;
}
bool cmp1(ll a, ll b){return g[a]<g[b];}
bool cmp2(ll a, ll b){return quiry[a][3]<quiry[b][3];}
void init()
{
ll i, j;
mu[1]=1;
for(i=2;i<=maxn;i++)
{
if(!mark[i])prime[++prime[0]]=i,mu[i]=-1;
for(j=1;j<=prime[prime[0]] and i*prime[j]<=maxn;j++)
{
mark[i*prime[j]]=1;
if(i%prime[j]==0){mu[i*prime[j]]=0;break;}
mu[i*prime[j]]=-mu[i];
}
}
for(i=1;i<=maxn;i++)
for(j=i;j<=maxn;j+=i)g[j]+=i;
for(i=1;i<=maxn;i++)numg[i]=i;
sort(numg+1,numg+maxn+1,cmp1);
scanf("%lld",&T);
for(i=1;i<=T;i++)scanf("%lld%lld%lld",quiry[i]+1,quiry[i]+2,quiry[i]+3);
for(i=1;i<=T;i++)numq[i]=i;
sort(numq+1,numq+T+1,cmp2);
}
ll calc(ll n, ll m)
{
ll i, last, ans=0;
if(n>m)swap(n,m);
for(i=1;i<=n;i=last+1)
{
last=min(n/(n/i),m/(m/i));
ans+=(sum(last)-sum(i-1))*(n/i)*(m/i);
}
return ans;
}
void work()
{
ll i, p=1, j;
for(i=1;i<=T;i++)
{
for(;g[numg[p]]<=quiry[numq[i]][3] and p<=maxn;p++)
for(j=numg[p];j<=maxn;j+=numg[p])add(j,g[numg[p]]*mu[j/numg[p]]);
quiry[numq[i]][0]=calc(quiry[numq[i]][1],quiry[numq[i]][2]);
}
}
int main()
{
init();
work();
for(ll i=1;i<=T;i++)printf("%lld\n",quiry[i][0]%mod);
return 0;
}