loj #509. 「LibreOJ NOI Round #1」动态几何问题「反演」
#509. 「LibreOJ NOI Round #1」动态几何问题
题意:
给了一个初中几何题,使用三角函数、相似、勾股定理等技巧,最后变成了:
给定 \(n,m \leq 1.5\times 10^{16}\),求多少对有序二元组 \((x,y)\) 满足 \(x\in[1,n],y\in[1,m],xy\) 是完全平方数。
题解:
省略亿点点内容,直接跳到题解的算法 5
算法 5
不妨设 \(n\leq m\)。
容易发现令 \(x = da,y=d'b\),其中 \(a,b\) 分别是 \(x,y\) 最大的平方因子,也就意味着 \(d,d'\) 是形如 \(\prod p_i\) 的数。由于 \(xy=dd'ab\) 是完全平方数,能推出 \(d=d'\)。于是我们想到可以枚举无平方因子数 \(d\),然后计算多少个 \(x,y\) 是 \(d\) 的倍数,且 \(\frac{x}{d},\frac{y}{d}\) 是完全平方数。答案为:
稍微解释一下,\(\mu^2(i)\) 表示的是 \(i\) 是否含平方因子;\(1\) 到 \(x\) 中的完全平方数有 \(\lfloor\sqrt x\rfloor\) 个。
有一个重要的式子不得不提,即 \(\mu^2\) 的前缀和:
证明其实就是容斥,等价于证明 \(\mu^2(n) = \sum_{d^2 \mid n} \mu(d)\),若 \(n\) 有 \(k\) 个素因子的次数至少为 \(2\),\(n\) 被计算到的次数为 \(\sum_{i = 0}^k {k\choose i}(-1)^k=[k=0]\)。
于是我们可以 \(\mathcal O(\sqrt n)\) 计算 \(\mu^2\) 的前缀和了。
注意到 \(\left\lfloor \sqrt{\dfrac{n}{i}}\right\rfloor,\left\lfloor \sqrt{\dfrac{m}{i}}\right\rfloor\) 分别有 \(\mathcal O( n^{\frac{1}{3}})\) 和 \(\mathcal O(\min(n, m^{\frac{1}{3}}))\) 种取值,所以我们可以分段求:对于 \(i\leq m^{\frac{1}{3}}\) 求,再对于 \(\left\lfloor \sqrt{\dfrac{m}{i}}\right\rfloor \leq m^{\frac{1}{3}}\) 求。时间复杂度是多少?我们可以通过积分求,最后大概得到 \(n = m\) 时是 \(\mathcal O(n^{0.5} \log n)\)。
算法 6
如果我们处理了 \(\sqrt n\) 以内 \(\mu,\mu^2\) 的前缀和,那么 \(\mu\) 的前缀和 \(S(n)\) 就可以 \(\mathcal O(n^{\frac{1}{3}})\) 求了。
用积分精细分析一下复杂度大约是 \(\mathcal O(\sqrt n)\)。
算法 7
如果我们预处理 \(S = \min(n,m^{\frac{3}{7}})\) 内的 \(\mu,\mu^2\) 前缀和,使用杜教筛求 \(\mu\) 就可以在 \(\mathcal O(\max(n,m)^{\frac37})\) 时间内计算,具体推导见题解(太毒瘤了)
求解有很多技巧。
对于外层有个结论:
对于筛 \(\mu^2\),也有个技巧。前面说我们对于 \(i\leq m^{\frac{1}{3}}\) 求,再对于 \(\left\lfloor \sqrt{\dfrac{m}{i}}\right\rfloor \leq m^{\frac{1}{3}}\) 求。具体如何实现?先记录最大的 \(i,i^3 \leq m^{\frac{1}{3}}\),求出 \(\sum_{j\leq i} \mu(j)\left\lfloor {\dfrac{n}{j^2}}\right\rfloor\),并记录 \(t=\sum_{j\leq i} \mu(j)\) 。然后依次求 \(\left\lfloor {\dfrac{n}{j^2}}\right\rfloor = \left\lfloor {\dfrac{n}{(i+1)^2}}\right\rfloor...1\) 的答案,根据上面的结论右端点可以求出来记为 \(r\),我们只需给答案加上 \(-t+\sum_{j\leq r} \mu(j)\)。你会发现每个 \(\mu(i)\) 恰好被前缀和覆盖了 \(\left\lfloor {\dfrac{n}{i^2}}\right\rfloor\) 次,也就神仙地完成了计算。
#include <unordered_map>
#include <algorithm>
#include <cstdio>
#include <cmath>
using namespace std;
typedef long long ll;
const int N = 2e7 + 30;
int p[N >> 3], mu[N], sn, tot;
ll mu1[N], mu2[N];
bool tag[N];
void sieve() {
mu[1] = 1;
for(int i = 2; i <= sn; i++) {
if(!tag[i]) p[++ tot] = i, mu[i] = -1;
for(int j = 1; j <= tot && i * p[j] <= sn; j++) {
tag[i * p[j]] = 1;
if (i % p[j] == 0) break ;
mu[i * p[j]] = -mu[i];
}
}
for(int i = 1; i <= sn; i ++) {
mu1[i] = mu1[i - 1] + mu[i];
mu2[i] = mu2[i - 1] + (mu[i] != 0);
}
}
ll n, m;
unordered_map<int, int> Map;
ll Mu(int n) {
if(n <= sn) return mu1[n];
if(Map.count(n)) return Map[n];
int res = 1;
for(int i = 2, j; i <= n; i = j + 1) {
j = n / (n / i); res -= Mu(n / i) * (j - i + 1);
}
Map[n] = res;
return res;
}
ll Mu2(ll n) {
if (n <= sn) return mu2[n];
ll res = 0, i;
for(i = 1; i * i * i <= n; i ++) {
res += mu[i] * (n / (i * i));
}
ll tmp = Mu(i - 1);
for(ll j = n / (i * i); j >= 1; j --) {
res += Mu(sqrt(n / j)) - tmp;
}
return res;
}
int main() {
scanf("%lld%lld", &n, &m);
if(n > m) swap(n, m);
sn = min((ll) (pow(m, 3 / 7.0)), n); sieve();
ll ans = 0, la = 0, cur = 0;
for(ll l = 1, r; l <= n; l = r + 1) {
ll a = sqrt(n / l), b = sqrt(m / l);
r = min(n / (a * a), m / (b * b));
ans += ((cur = Mu2(r)) - la) * a * b;
la = cur;
}
printf("%lld\n", ans);
return 0;
}