看到题面想到了莫比乌斯反演模板题，即求(sum_{x = 1}^{N}sum_{y = 1}^{M} [gcd(x, y) = k])，其中(k)为指定常数。
这题将(k)的取值变为了素数，不过还是按照常规套路推一下式子：
[egin{aligned}
ans
&= sum_{p in {P}} sum_{x = 1}^{N} sum_{y = 1}^{M} [gcd(x, y) = p] \
&= sum_{p in {P}} sum_{x = 1}^{lfloor frac{N}{p}
floor} sum_{y = 1}^{lfloor frac{M}{p}
floor} [gcd(x, y) = 1] \
end{aligned}
]
根据(varepsilon = mu ast 1)有：
[egin{aligned}
ans
&= sum_{p in {P}} sum_{x = 1}^{lfloor frac{N}{p}
floor} sum_{y = 1}^{lfloor frac{M}{p}
floor} [gcd(x, y) = 1] \
&= sum_{p in {P}} sum_{x = 1}^{lfloor frac{N}{p}
floor} sum_{y = 1}^{lfloor frac{M}{p}
floor} sum_{d mid gcd(x, y)} mu(d) \
&= sum_{p in {P}} sum_{d} mu(d) sum_{x = 1}^{lfloor frac{N}{p}
floor} [d mid x] sum_{y = 1}^{lfloor frac{M}{p}
floor} [d mid y] \
&= sum_{p in {P}} sum_{d} mu(d) lfloor frac{N}{pd}
floor lfloor frac{M}{pd}
floor
end{aligned}
]
到这一步常规套路就推不下去了，然后根据这个式子枚举(p)去算的话铁定会超时。
令(T = pd)，有
[egin{aligned}
ans
&= sum_{p in {P}} sum_{d} mu(d) lfloor frac{N}{pd}
floor lfloor frac{M}{pd}
floor \
&= sum_{T} lfloor frac{N}{T}
floor lfloor frac{M}{T}
floor sum_{t mid T, t in P} mu(frac{T}{t})
end{aligned}
]
前半部分可以数论分块，后半部分可以线性筛预处理。
现在就可以(O(N) sim O(sqrt{N}))做了。