NKOJ 2841 (SDOI 2014)数表(莫比乌斯反演+树状数组+线性筛)

P2841【SDOI2014 R1D1】数表

问题描述

有一张n*m的数表,其第i行第j列(1<=i<=n,1<=j<=m)的数值为能同时整除i和j的所有自然数之和。给定a,计算数表中不大于a的数之和。

输入格式

输入包含多组数据。
输入的第一行一个整数Q表示测试点内的数据组数,接下来Q行,每行三个整数n,m,a描述一组数据。

输出格式

对每组数据,输出一行一个整数,表示答案模2^31的值。

样例输入

2
4 4 3
10 10 5

样例输出

20
148

提示

对于30%的数据,1<=n,m<=400,1<=Q<=200
对于另外30%的数据,1<=n,m<=10^5,1<=Q<=10
对于100%的数据,1<=n,m<=10^5,1<=Q<=20000,0<=a<=10^9


容易得到对于单次询问
$$
Ans=\sum_{i=1}^{n}\sum_{j=1}^{m}\sigma(gcd(i,j)),[\sigma(gcd(i,j))<=a]
$$
那么,推导一番
$$
Ans=\sum_{d=1}^{N}\sigma(d)\sum_{i=1}^{n}\sum_{j=1}^{m}1[gcd(i,j)=d],N=min(n,m)
$$
然后,反演一下
$$
Ans=\sum_{d=1}^{N}\sigma(d)\sum_{d|K}^{N}\mu(\frac{K}{d})\lfloor{\frac{n}{K}}\rfloor\lfloor{\frac{m}{K}}\rfloor=\sum_{K=1}^{N}\lfloor{\frac{n}{K}}\rfloor\lfloor{\frac{m}{K}}\rfloor\sum_{d|K}\sigma(d)\mu(\frac{K}{d}),\sigma(d)<=a
$$
注意到$\sum_{d|K}\sigma(d)\mu(\frac{K}{d})$是个积性函数,用线性筛预处理,但是只有$\sigma(d)<=a$的项才对答案有贡献,因此离线处理,对询问按照a排序,然后每次将小于等于a的那些项暴力枚举一下倍数,更新到树状数组里面,查询的时候直接用。
复杂度为$O(n\ ln\ n+m\sqrt{n}log_2n)$左右


代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
#include<stdio.h>
#include<iostream>
#include<algorithm>
#include<cstring>
#define ll long long
#define N 100005
using namespace std;
const ll mod=(1ll<<31);
struct node{ll n,m,a,id;}K[N],g[N];
bool cmp(node aa,node bb)
{return aa.a<bb.a;}
ll T,P[N],tot,mu[N],pc[N],pd[N],G[N],Ans[N];
void EU()
{
ll i,j;mu[1]=1;g[1].a=1;
for(i=1;i<N;i++)g[i].n=i;
for(i=2;i<N;i++)
{
if(!g[i].a)P[++tot]=i,mu[i]=-1,g[i].a=i+1,pc[i]=i+1,pd[i]=tot;
for(j=1;j<=tot&&P[j]*i<N;j++)
if(i%P[j])
{
pc[i*P[j]]=P[j]+1;
pd[i*P[j]]=j;
mu[i*P[j]]=-mu[i];
g[i*P[j]].a=g[i].a*pc[i*P[j]]%mod;
}
else
{
pc[i*P[j]]=pd[i]==j?pc[i]*P[j]+1:P[j]+1;
pd[i*P[j]]=j;
g[i*P[j]].a=pd[i]==j?g[i].a/pc[i]*pc[i*P[j]]:g[i].a*pc[i*P[j]];
pc[i*P[j]]%=mod;
g[i*P[j]].a%=mod;
}
}
sort(g+1,g+N,cmp);
}
void MD(ll x,ll d)
{for(ll i=x;i<N;i+=(i&-i))G[i]+=d,G[i]%=mod;}
ll GS(ll x)
{
ll i,sum=0;
for(i=x;i;i-=(i&-i))sum+=G[i],sum%=mod;
return sum;
}
int main()
{
ll i,j,p,k=1,ans;EU();
scanf("%lld",&T);
for(i=1;i<=T;i++)scanf("%lld%lld%lld",&K[i].n,&K[i].m,&K[i].a),K[i].id=i;
sort(K+1,K+T+1,cmp);
for(i=1;i<=T;i++)
{
while(g[k].a<=K[i].a)
{
for(j=g[k].n;j<N;j+=g[k].n)MD(j,g[k].a*mu[j/g[k].n]%mod);
k++;
}
ans=0;
for(j=1;j<=K[i].n&&j<=K[i].m;j=p+1)
{
p=min(K[i].n/(K[i].n/j),K[i].m/(K[i].m/j));
ans+=(K[i].n/j)*(K[i].m/j)%mod*(GS(p)-GS(j-1))%mod;ans%=mod;
}
Ans[K[i].id]=(ans+mod)%mod;
}
for(i=1;i<=T;i++)printf("%lld\n",Ans[i]);
}