Problem
給定 n 和 m,請你統計有序對 (a,b) 的個數,其中 1≤a≤n,1≤b≤m 且 a 與 b 互質。
a 與 b 互質的定義是:a 與 b 的最大公約數等於 1。> count coprime pair
## Sample Input ##
|
|
Sample Output
|
|
Solution
前言
柳柳州出數論基礎-「莫比烏斯反演 (Möbius inversion formula)」
之前放棄看完的內容又要撿回來看,百般痛苦,之前沒看懂啊。但能知道的是莫比烏斯反演類似排容原理,就像錯排依樣。先不管莫比烏斯,看到「線性時間求出所有乘法反元素」真的可以嗎?那也是件很有趣的事情。
參考資源
基礎定義
了解莫比烏斯反演之前,要先介紹積性函數
積性函數
- 何謂積性函數 (Multiplicative function)?
對於定義域 N+ 的函數 f 而言,任兩個互質的 gcd(a,b)=1 正整數 a,b 滿足 f(ab)=f(a)f(b)。 - 何謂完全積性函數?
f 是一個積性函數,同時滿足 f(pn)=f(p)n。 - 若 f(n),g(n) 是積性函數,則 h(n)=f(n)g(n) 也是積性函數。
- 若 f(n) 為積性函數,則函數 g(n)=∑d|nf(d) 也是積性函數。
歐拉函數
回顧歐拉函數 (Euler’s totient function) ϕ
- 定義 ϕ(n) 為 1⋯n 中與 n 互質的個數。
ϕ(n) 是一個積性函數,但不是完全積性函數。根據中國餘數定理 或者 基本算術可以證明之。
歐拉定理 a^{\phi(n)} \equiv 1 \mod n, \text{ when } gcd(a, n) = 1。
特性 1:\sum_{d|n} \phi(d) = n。當作把互質個數 \phi(n) 相加不互質個數 s,由於 d 是 n 的因數,則與 d 互質 x 個數 \phi(d),把那些 x' = x \times d/n 就會補上那些不互質的個數 s = |set(x')|。
特性 2:1 \cdots n 與 n 互質的數和為 n\phi(n)/2。原因很簡單,若 gcd(x, n) = 1,則會滿足 gcd(n-x, n) = 1,看起來就是一個對稱總和。
莫比烏斯
根據積性函數的性質,我們得到莫比烏斯反演的基礎:
莫比烏斯反演公式 f(n)
f(n) = \sum_{d|n} \mu(d) g(\frac{n}{d})莫比烏斯函數 \mu(n)
\mu(n) = \left\{\begin{matrix} 1 && n = 1\\ (-1)^k && n = p_1 p_2 \cdots p_k \\ 0 && \text{otherwise} \end{matrix}\right.- 特性 1:\sum_{d|n} \mu(d) = [n = 1]
在此簡單說一次,莫比烏斯反演公式就像排容原理,而莫比烏斯函數 \mu(n) 就像一加一減,之所以在 n = p_1 p_2 \cdots p_k 為 \mu(n) = (-1)^k,也就是說當 n 只全由質數相乘 (不允許冪次方大於 1,否則為 0),相當於一般認知,奇數次要扣,偶數次要補回來的排容口訣。而特定自然不必多說,其中數學表達 [n = 1] 相當於程式中的 n == 1 ? 1 : 0
。
簡單展示歐拉函數 \phi(n)
- f(n) = \sum_{d|n} \phi(d) = n,由歐拉函數特性 2 得知。
*\phi(n) = \sum_{d|n} \mu(d) f(\frac{n}{d}) = \sum_{d|n} \frac{\mu(d) n}{d},套用莫比烏斯公式後整理。
關於公式\phi(n) = \sum_{d|n} \frac{\mu(d) n}{d} 可以這麼理解。
- 求出 gcd(n, k) = 1 的個數,其餘要捨棄掉。
- 計算 gcd(n, k) = d 的個數,利用排容原理即莫比烏斯函數,顯而易見排容方案只當 d 為 n 個質因數組合而成。
應用此題
利用莫比烏斯函數 特性 1:\sum_{d|n} \mu(d) = [n = 1]
\begin{align*} & \sum_{a = 1}^{N} \sum_{b = 1}^{M} [gcd(a, b) = 1] \\ & = \sum_{a = 1}^{N} \sum_{b = 1}^{M} \sum_{d|gcd(a, b)} \mu(d) \\ & = \sum \mu(d) \left ( \sum_{1 \le a \le N \text{ and } d | a} \left ( \sum_{1 \le b \le M \text{ and } d | b} 1 \right ) \right )\\ & = \sum \mu(d) \left \lfloor \frac{n}{d} \right \rfloor \left \lfloor \frac{m}{d} \right \rfloor \end{align*}第二行就是抽換莫比烏斯,第三行則是交換順序,第四行則是可以快速找到 d|a 的總數為 \lfloor n/d \rfloor 所導致。
若直接窮舉 d 時間複雜度為 O(min(n, m)),可以利用塊狀的概念優化到 O(\sqrt{n} + \sqrt{m})。因為 \lfloor n/d \rfloor 的值只有 2\lfloor \sqrt{n} \rfloor 種可能,同理 \lfloor m/d \rfloor 也是,那麼對於同一個 d 的計數會是一個連續區間不斷地移動。
以下是一個簡單的 \lfloor n/d \rfloor 示意圖,用 \lfloor n/d \rfloor 不同當作劃分。
|
|
程式只要處理打 X
的位置即可,時間複雜度 O(\sqrt{n} + \sqrt{m}),可以參照一般解。代碼短,但除法次數多。
加速一般解由 liouzhou_101 提供。開一次根號 sqrt()
,省下 O(\sqrt{n}) 次的除法,賺了 400 ms,代碼長了 800 B。加速 25% 的效能,可謂除法的可怕,根據研究差異後才發現到這一點。
一般解
|
|
加速運算解
|
|