JZOJ 4379 盈虚方田

这不原题???
手动左转「『SDOI2014』数表」……

首先单纯度可以在欧拉筛的时候顺便筛出来。
然后设 \(n\) 的单纯度为 \(f(n)\)

\[\begin{align*} & \sum\limits_{i = 1}^n \sum\limits_{j = 1}^m [f(\gcd(i,j)) \le p] \\ = & \sum\limits_{d = 1}^{\min(n,m)} [f(d) \le p] \sum\limits_{id \le \min(n,m)} \mu(i) \lfloor \dfrac n {id} \rfloor \lfloor \dfrac m {id} \rfloor \\ = & \sum\limits_{T = 1}^{\min(n,m)} \lfloor \dfrac n T \rfloor \lfloor \dfrac m T \rfloor \sum\limits_{d | T} [f(d) \le p] \mu(\dfrac T d) \end{align*}\]

然后套路地套上一个树状数组,复杂度 \(O(n \log^2 n + T \sqrt n \log n)\)(认为所有 \(n,m\) 同阶)。

代码:

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
#include <cstdio>
#include <algorithm>
#define lowbit(x) ((x) & -(x))
using namespace std;

const int BUFF_SIZE = 1 << 20;
char BUFF[BUFF_SIZE],*BB,*BE;
#define gc() (BB == BE ? (BE = (BB = BUFF) + fread(BUFF,1,BUFF_SIZE,stdin),BB == BE ? EOF : *BB++) : *BB++)
template<class T>
inline void read(T &x)
{
x = 0;
char ch = 0,w = 0;
for(;ch < '0' || ch > '9';w |= ch == '-',ch = gc());
for(;ch >= '0' && ch <= '9';x = (x << 3) + (x << 1) + (ch ^ '0'),ch = gc());
w ? x = -x : x;
}

const int N = 5e5;
const int Q = 5e3;
int q;
struct s_query
{
int n,m,p,id;
inline bool operator<(const s_query &o) const
{
return p < o.p;
}
} qry[Q + 5];
int vis[N + 5],prime[N + 5],cnt,mu[N + 5],f[N + 5];
struct note
{
int id,f;
inline bool operator<(const note &o) const
{
return f < o.f;
}
} a[N + 5];
int c[N + 5];
inline void update(int x,int k)
{
for(;x <= N;x += lowbit(x))
c[x] += k;
}
inline int query(int x)
{
long long ret = 0;
for(;x;x -= lowbit(x))
ret += c[x];
return ret;
}
long long ans[N + 5];
int main()
{
mu[1] = 1;
for(register int i = 2;i <= N;++i)
{
if(!vis[i])
mu[prime[++cnt] = i] = -1,f[i] = 1;
for(register int j = 1;j <= cnt && i * prime[j] <= N;++j)
{
vis[i * prime[j]] = 1,f[i * prime[j]] = f[i] + 1;
if(!(i % prime[j]))
break;
else
mu[i * prime[j]] = -mu[i];
}
}
for(register int i = 1;i <= N;++i)
a[i].id = i,a[i].f = f[i];
sort(a + 1,a + N + 1);
read(q);
for(register int i = 1;i <= q;++i)
read(qry[i].n),read(qry[i].m),read(qry[i].p),qry[i].id = i;
sort(qry + 1,qry + q + 1);
for(register int i = 1,n,m,j = 1;i <= q;++i)
{
for(;a[j].f <= qry[i].p && j <= N;++j)
for(register int k = a[j].id,l = 1;k <= N;k += a[j].id,++l)
update(k,mu[l]);
n = qry[i].n,m = qry[i].m;
if(n > m)
swap(n,m);
for(register int l = 1,r;l <= n;l = r + 1)
{
r = min(n / (n / l),m / (m / l));
ans[qry[i].id] += (long long)(query(r) - query(l - 1)) * (n / l) * (m / l);
}
}
for(register int i = 1;i <= q;++i)
printf("%lld\n",ans[i]);
}