隐藏
「SDOI2014」数表 - 莫比乌斯反演+树状数组 | Bill Yang's Blog

路终会有尽头,但视野总能看到更远的地方。

0%

「SDOI2014」数表 - 莫比乌斯反演+树状数组

题目大意

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


题目分析

首先我们忽略$a$的限制。
令$f(x)$为$x$的约数之和,$sum(x)$为$x=\gcd(i,j),1\le i\le n,1\le j\le m$的$(i,j)$对数。

故:

令$g(x)=\sum\limits_{d\mid x}^{\min(n,m)}f(d)\mu(\frac xd)$
则:

对于$g(x)=\sum\limits_{d\mid x}^{\min(n,m)}f(d)\mu(\frac xd)$,因为$x$的约数只有$\sqrt n$个,故可以暴力求。

现在考虑题目$a$的限制。
也就是说当$d\gt a$时,$f(d)=0$。
那么也就是$g(x)$的值与限制$a$有关。

我们可以将询问离线下来,按照$a$升序排序。
每次加入$d\le a$的$f(d)$对$g$的影响,然后分块计算$ans$。
因此这是一个单点修改,询问前缀和的问题,可以采用树状数组实现。

如何快速求出$f(x)$?
可以通过预处理筛出$f$,可以枚举倍数暴力筛,也可以使用线性筛,约数和线性筛见这里


代码

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
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
#include<algorithm>
#include<iostream>
#include<iomanip>
#include<cstring>
#include<cstdlib>
#include<climits>
#include<vector>
#include<cstdio>
#include<cmath>
#include<queue>
using namespace std;

#define mp make_pair
#define pii pair<int,int>

inline const int Get_Int() {
int num=0,bj=1;
char x=getchar();
while(x<'0'||x>'9') {
if(x=='-')bj=-1;
x=getchar();
}
while(x>='0'&&x<='9') {
num=num*10+x-'0';
x=getchar();
}
return num*bj;
}

const int maxn=100005,maxq=20005;

int t,cnt=0,Max,u[maxn],F[maxn],Min_Fac_last[maxn],Min_Fac_sum[maxn],Prime[maxn],ans[maxn];
pii f[maxn];
bool vst[maxn];

void Sieve(int n) {
u[1]=F[1]=1;
for(int i=2; i<=n; i++) {
if(!vst[i]) {
Prime[++cnt]=i;
u[i]=-1;
Min_Fac_last[i]=i;
F[i]=Min_Fac_sum[i]=i+1;
}
for(int j=1; j<=cnt&&i*Prime[j]<=n; j++) {
vst[i*Prime[j]]=1;
if(i%Prime[j]==0) {
u[i*Prime[j]]=0;
Min_Fac_last[i*Prime[j]]=Min_Fac_last[i]*Prime[j];
Min_Fac_sum[i*Prime[j]]=Min_Fac_sum[i]+Min_Fac_last[i*Prime[j]];
F[i*Prime[j]]=F[i]/Min_Fac_sum[i]*Min_Fac_sum[i*Prime[j]];
break;
}
u[i*Prime[j]]=-u[i];
F[i*Prime[j]]=F[i]*(Prime[j]+1);
Min_Fac_last[i*Prime[j]]=Prime[j];
Min_Fac_sum[i*Prime[j]]=Prime[j]+1;
}
}
for(int i=1; i<=n; i++)f[i]=mp(F[i],i);
}

struct BIT {
int c[maxn];
#define lowbit(x) x&(-x)
void add(int x,int v) {
for(int i=x; i<=Max; i+=lowbit(i))c[i]+=v;
}
int sum(int x) {
int ans=0;
for(int i=x; i>=1; i-=lowbit(i))ans+=c[i];
return ans;
}
} bit;

struct Query {
int n,m,a,id;
bool operator < (const Query& b) const {
return a<b.a;
}
} q[maxq];

void Solve(int num) {
int n=q[num].n,m=q[num].m;
for(int i=1,next; i<=n; i=next+1) {
next=min(n/(n/i),m/(m/i));
ans[q[num].id]+=(n/i)*(m/i)*(bit.sum(next)-bit.sum(i-1));
}
}

int main() {
t=Get_Int();
for(int i=1; i<=t; i++) {
q[i].n=Get_Int();
q[i].m=Get_Int();
q[i].a=Get_Int();
q[i].id=i;
if(q[i].n>q[i].m)swap(q[i].n,q[i].m);
Max=max(Max,q[i].n);
}
Sieve(Max);
sort(q+1,q+t+1);
sort(f+1,f+Max+1);
int now=0;
for(int i=1; i<=t; i++) {
while(now<Max&&f[now+1].first<=q[i].a)
for(int j=f[++now].second; j<=Max; j+=f[now].second)bit.add(j,f[now].first*u[j/f[now].second]);
Solve(i);
}
for(int i=1; i<=t; i++)printf("%d\n",ans[i]&INT_MAX);
return 0;
}
姥爷们赏瓶冰阔落吧~