1. 程式人生 > 實用技巧 >loj #509. 「LibreOJ NOI Round #1」動態幾何問題「反演」

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}\) 是完全平方數。答案為:

\[\sum_{d = 1}^{n} \mu^2(i) \left\lfloor \sqrt{\left\lfloor\dfrac{n}{i}\right\rfloor}\right\rfloor \left\lfloor \sqrt{\left\lfloor\dfrac{m}{i}\right\rfloor}\right\rfloor \]

稍微解釋一下,\(\mu^2(i)\) 表示的是 \(i\) 是否含平方因子;\(1\)\(x\) 中的完全平方數有 \(\lfloor\sqrt x\rfloor\) 個。

有一個重要的式子不得不提,即 \(\mu^2\) 的字首和:

\[\sum_{i = 1}^n \mu^2(i) = \sum_{i = 1}^{\lfloor\sqrt n\rfloor} \mu(i)\left\lfloor {\dfrac{n}{i^2}}\right\rfloor \]

證明其實就是容斥,等價於證明 \(\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})\) 時間內計算,具體推導見題解(太毒瘤了)

求解有很多技巧。

對於外層有個結論:

\[令\left\lfloor \sqrt{\left\lfloor\dfrac{n}{i}\right\rfloor}\right\rfloor=a,則這一段右端點是 \left\lfloor\dfrac{n}{a^2}\right\rfloor \]

對於篩 \(\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;
}