跳转至

連分數

連分數

連分數 是實數作為有理數的特定收斂序列的表示。它們在算法競賽（competitive programming）中很有用，因為它們易於計算，並且可以有效地用於在分母不超過給定值的所有數字中，找到基礎實數（underlying real number）的最佳可能有理近似（best possible rational approximation）。

除此之外，連分數與歐幾里得算法密切相關，這使得它們在一系列數論問題中非常有用。

定義

連分數是一種記號。例如，長為 \(4\) 的連分數：

\[ [a_0,a_1,a_2,a_3]=a_0+\frac{1}{a_1+\frac{1}{a_2+\frac{1}{a_3}}} \]

只是為了形式上簡潔，才記成等號左邊的樣子。這裏的四個變元可以任意取值。

連分數各變元的下標從 \(0\) 開始。

簡單連分數

可以證明，任何有理數都可以精確地以兩種方式表示為連分數：

\[ r = [a_0;a_1,\dots,a_k,1] = [a_0;a_1,\dots,a_k+1] \]

此外，對於 \(r=\frac{p}{q}\)，這種連分數的長度 \(k\) 估計為 \(k = O(\log \min(p, q))\)。

一旦深入研究了連分數構造的細節，這背後的原因就會很清楚。

定義

設 \(a_0, a_1, \dots, a_k \in \mathbb Z\) 和 \(a_1, a_2, \dots, a_k \geq 1\)。然後表達式

\[ r=a_0 + \frac{1}{a_1 + \frac{1}{\dots + \frac{1}{a_k}}} \]

稱為有理數 \(r\) 的 連分數表示，並簡短地表示為 \(r=[a_0;a_1,a_2,\dots,a_k]\)。

示例

設 \(r = \frac{5}{3}\)。有兩種方法可以將其表示為連分數：

\[ \begin{align} r = [1;1,1,1] &= 1+\frac{1}{1+\frac{1}{1+\frac{1}{1}}},\\ r = [1;1,2] &= 1+\frac{1}{1+\frac{1}{2}}. \end{align} \]

對於有限連分數，全體結尾為 \(1\) 的有限連分數和全體結尾不為 \(1\) 的有限連分數一一對應，即同一個連分數有兩種表示：

\([a_0,a_1,a_2,a_3]=[a_0,a_1,a_2,a_3-1,1]\)

簡單連分數：連分數從第 \(1\) 項開始全都是正整數。如果有限，要求最後一項不為 \(1\)。（第 \(0\) 項可以任意）

簡單連分數的值，一定大於偶數的漸進分數，一定小於奇數的漸進分數。無限簡單連分數一定收斂。

仿照一般分數的概念，第 \(0\) 項是 \(0\) 的連分數稱為「真分數」。顯然如果這之後的所有變元都大於等於 \(1\)，那麼得到的真分數一定落在 \(0\) 到 \(1\) 之間。

無限連分數

如果分式無限地寫下去，有無限個變元，就得到無限連分數。無限連分數收斂等價於漸進分數收斂。

有定理：

無限連分數，如果各變元均大於等於 \(1\)，那麼一定收斂。

因為只要各變元為正，無限連分數的偶漸進分數單調遞增（都比它小），奇漸進分數單調遞減（都比它大）。而在均大於等於 \(1\) 時，相鄰（奇偶間）兩個漸進分數之間距離可以給出估計式，趨於 \(0\)，因此收斂。

顯然可以看到，連分數關於下標為偶數的變元單調遞增，關於下標為奇數的變元單調遞減。這無論它有限或無限都成立。

定義

設 \(a_0,a_1,a_2, \dots\) 為整數序列，使得 \(a_1, a_2, \dots \geq 1\)。設 \(r_k = [a_0; a_1, \dots, a_k]\)。然後表達式

\[ r = a_0 + \frac{1}{a_1 + \frac{1}{a_2+\dots}} = \lim\limits_{k \to \infty} r_k \]

稱為無理數 \(r\) 的 連分數表示，並簡短地表示為 \(r = [a_0;a_1,a_2,\dots]\)。

注意，對於 \(r=[a_0;a_1,\dots]\) 和整數 \(k\)，有 \(r+k = [a_0+k; a_1, \dots]\)。

另一個重要的觀察結果是，當 \(a_0>0\) 時，\(\frac{1}{r}=[0;a_0, a_1, \dots]\)；當 \(a_1＝0\) 時，\(\frac{1}{r} = [a_1; a_2, \dots]\)。

漸進分數

定義

在上面的定義中，有理數 \(r_0, r_1, r_2, \dots\) 稱為 \(r\) 的 漸進分數（convergents，意為「收斂」）。

相應地，單個 \(r_k = [a_0; a_1, \dots, a_k] = \frac{p_k}{q_k}\) 稱為 \(r\) 的第 \(k\) 個漸進分數。

示例

考慮 \(r = [1; 1, 1, 1, \dots]\)。可以通過歸納法證明 \(r_k = \frac{F_{k+2}}{F_{k+1}}\)，其中 \(F_k\) 是斐波那契序列，定義為 \(F_0 = 0\)、\(F_1 = 1\) 和 \(F_{k} = F_{k-1} + F_{k-2}\)。從 Binet 公式可知

\[ r_k = \frac{\phi^{k+2} - \psi^{k+2}}{\phi^{k+1} - \psi^{k+1}} \]

其中 \(\phi = \frac{1+\sqrt{5}}{2} \approx 1.618\) 是黃金比率，\(\psi = \frac{1-\sqrt{5}}{2} = -\frac{1}{\phi} \approx -0.618\)。因此

\[ r = 1+\frac{1}{1+\frac{1}{1+\dots}}=\lim\limits_{k \to \infty} r_k = \phi = \frac{1+\sqrt{5}}{2}. \]

請注意，在這種特定情況下，找到 \(r\) 的另一種方法是求解方程

\[ r = 1+\frac{1}{r} \implies r^2 = r + 1 \]

定義

設 \(r_k = [a_0; a_1, \dots, a_{k-1}, a_k]\)。對於 \(1 \leq t \leq a_k\)，\([a_0; a_1, \dots, a_{k-1}, t]\) 稱為 中間分數（semiconvergents，「semi」意為「半」）。

通常將大於 \(r\) 的分數稱為 上（upper）漸進分數或中間分數，將小於 \(r\) 者稱為 下（lower）漸進分數或中間分數。

餘項和部分商

定義

作為漸進分數的補充，定義 餘項（完全商，complete quotients）為 \(s_k = [a_k; a_{k+1}, a_{k+2}, \dots]\)。

相應地，將單個 \(s_k\) 稱為 \(r\) 的第 \(k\) 個完全商。

相應地，取整得到的 \(a_k\) 稱為部分商。

對於各項均為正數的連分數，所有的餘項也都是正數。

根據以上定義，可以得出 \(s_k \geq 1\) 代表 \(k \geq 1\) 的結論。

將 \([a_0; a_1, \dots, a_k]\) 視為一個形式代數表達式，並允許任意實數代替 \(a_i\)，得到

\[ r = [a_0; a_1, \dots, a_{k-1}, s_k] \]

特別地，\(r = [s_0] = s_0\)。另一方面，可以將 \(s_k\) 表示為

\[ s_k = [a_k; s_{k+1}] = a_k + \frac{1}{s_{k+1}} \]

這意味着可以從 \(s_k\) 計算 \(a_k = \lfloor s_k \rfloor\) 和 \(s_{k+1} = (s_k - a_k)^{-1}\)。

序列 \(a_0, a_1, \dots\) 定義良好，除非 \(s_k=a_k\)，這僅在 \(r\) 為有理數時發生。

因此，對於任何無理數 \(r\)，連分數表示都是唯一的。

求解簡單連分數表示

在代碼片段中主要假設有限的連分數。

如果要求 \((0,1)\) 區間內某個數的簡單連分數表示（第 \(0\) 項為 \(0\)），只需：

  • 取倒數，得到的餘項大於 \(1\)。
  • 取整得到整數部分為部分商，小數部分在 \(0\) 到 \(1\) 之間。
  • 對小數部分重複上述操作。

這樣就得到了相應的表示。

從 \(s_k\) 到 \(s_{k+1}\) 的轉換如下

\[ s_k =\left\lfloor s_k \right\rfloor + \frac{1}{s_{k+1}} \]

從這個表達式中，下一個完全商 \(s_{k+1}\) 如下

\[ s_{k+1} = \left(s_k-\left\lfloor s_k\right\rfloor\right)^{-1} \]

對於 \(s_k=\frac{p}{q}\)，這意味着

\[ s_{k+1} = \left(\frac{p}{q}-\left\lfloor \frac{p}{q} \right\rfloor\right)^{-1} = \frac{q}{p-q\cdot \lfloor \frac{p}{q} \rfloor} = \frac{q}{p \bmod q}. \]

因此，\(r=\frac{p}{q}\) 的連分數表示的計算遵循 \(p\) 和 \(q\) 的歐幾里得算法的步驟。

由此，\(\frac{p_k}{q_k} = [a_0; a_1, \dots, a_k]\) 的 \(\gcd(p_k, q_k) = 1\)。因此，漸進分數總是不可約的。

1
2
3
4
5
6
7
8
auto fraction(int p, int q) {
  vector<int> a;
  while (q) {
    a.push_back(p / q);
    tie(p, q) = make_pair(q, p % q);
  }
  return a;
}

1
2
3
4
5
6
def fraction(p, q):
    a = []
    while q:
        a.append(p // q)
        p, q = q, p % q
    return a

如果規定第 \(0\) 項是該數的取整，那麼全體實數都有「唯一的簡單連分數表示」。其中：

如果兩個無限簡單連分數的值相等，必然逐項相等。

如果兩個有限簡單連分數的值相等，不僅要逐項相等，而且必然項數也相同。

無限簡單連分數不能與有限簡單連分數值相等。有理數與有限簡單連分數具有一一對應關係，因此無限簡單連分數全都是無理數。

性質

為了給連分數的進一步研究提供一些動力，現在給出一些性質。

遞推

遞推

對於漸進分數 \(r_k = \frac{p_k}{q_k}\)，以下遞推公式適用於它們的快速計算：

\[ \frac{p_k}{q_k}=\frac{a_k p_{k-1} + p_{k-2}}{a_k q_{k-1} + q_{k-2}} \]

其中 \(\frac{p_{-1}}{q_{-1}}=\frac{1}{0}\) 並且 \(\frac{p_{-2}}{q_{-2}}=\frac{0}{1}\)。

漸進分數分子和分母具有完全相同的遞推關係：

\[ p_k=a_kp_{k-1}+p_{k-2} \]
\[ q_k=a_kq_{k-1}+q_{k-2} \]

這裏和 Farey 數列的遞推關係很像。

形式上記初項：

\[ p_{-1}=1 \quad p_0=a_0 \]
\[ q_{-1}=0 \quad q_0=1 \]

只是形式上成立。第 \(-1\) 項漸進分數是 1/0，沒有實際意義。

證明

可以注意到，\(p_k\) 與 \(q_k\) 對於 \(a_k\) 和 \(b_k\) 都是線性函數。這是因為，\(a_k\) 和 \(b_k\) 都只出現了一次，無論如何通分也不會有另一個 \(a_k\) 或 \(b_k\) 乘上去。於是通過待定係數，即可解得這個遞推關係。

反序定理

\[ \frac{q_k}{q_{k-1}}=[a_k,a_{k-1},\ldots,a_1] \]

如果 \(a_0\neq 0\)：

\[ \frac{p_k}{p_{k-1}}=[a_k,a_{k-1},\ldots,a_0] \]

如果 \(a_0=0\)：

\[ \frac{p_k}{p_{k-1}}=[a_k,a_{k-1},\ldots,a_2] \]
證明

對遞推關係稍加改造，有：

\[ \frac{p_k}{p_{k-1}}=a_k+\frac{1}{\frac{p_{k-1}}{p_{k-2}}} \]
\[ \frac{q_k}{q_{k-1}}=a_k+\frac{1}{\frac{q_{k-1}}{q_{k-2}}} \]

又利用初值，即可證明反序定理。

漸進分數的差分

計算相鄰兩項漸進分數的差，需要通分。通分後的分子代入遞推關係：

\[ p_{k+1}q_k-q_{k+1}p_k=\left(a_{k+1}p_k+p_{k-1}\right)q_k-\left(a_{k+1}q_k+q_{k-1}\right)p_k=-\left(p_kq_{k-1}-q_kp_{k-1}\right) \]

代入初值就有漸進分數的差分：

\[ p_{k+1}q_k-q_{k+1}p_k=(-1)^k \]
\[ \frac{p_{k+1}}{q_{k+1}}-\frac{p_k}{q_k}=\frac{(-1)^k}{q_{k+1}q_k} \]
注

可以觀察到，式 \(p_{k+1}q_k-q_{k+1}p_k\) 特別像一個行列式，完全可以按「行列式」理解。

漸進分數的遞推關係很像行列式的列變換。行列式一列加到另一列上不改變它的值，兩列交換則反號。

根據遞推式，如果連分數各項均為整數，則漸進分數分子分母總是互素。

對於有理數的簡單連分數展開，常用漸進分數差分的等式，求解一次線性不定方程（參見 擴展歐幾里得算法）：

\[ p_{k+1}q_k-q_{k+1}p_k=(-1)^k \]
\[ ax-by=1 \]

因為 a 與 b 互素，\(\frac{a}{b}\) 就是最簡的有理數，也就是它本身的最後一個漸進分數。那麼，它的前一個漸進分數就是所求的解。

倒數定理

由於實數與簡單連分數一一對應，稱實數的簡單連分數的漸進分數，就是實數的漸進分數。於是就有倒數定理：

對於大於 1 的實數 x，x 的漸進分數的倒數恰好是 \(\frac{1}{x}\) 的漸進分數。顯然，該定理也應該對於 0 到 1 之間的實數 x 成立。

證明
\[ x=[a_0,a_1,a_2,\ldots] \]
\[ \frac{1}{x}=\frac{1}{[a_0,a_1,a_2,\ldots]}=[0,a_0,a_1,a_2,\ldots] \]

於是根據新的初值與遞推就能發現倒數關係成立。

最佳逼近

最佳逼近

讓 \(\frac{p}{q}\) 是最小化 \(\left|r-\frac{p}{q}\right|\) 的分數，對於某些 \(x\)，該分數服從 \(q \leq x\)。

那麼 \(\frac{p}{q}\) 是 \(r\) 的中間分數。

因此允許通過檢查 \(r\) 是否中間分數來找到其最佳有理逼近。

下面會對這些性質建立一些直覺並做出進一步解釋。

連行列式

繼續看前面定義的漸近分數。對於 \(r=[a_0, a_1, a_2, \dots]\)，其漸近分數為

\[ \begin{gathered} r_0=[a_0] \\ r_1=[a_0, a_1] \\ \dots \\ r_k=[a_0, a_1, \dots, a_k] \end{gathered} \]

漸近分數是連分數的核心概念，因此研究它們的性質很重要。

對於數字 \(r\)，其第 \(k\) 個漸近分數 \(r_k = \frac{p_k}{q_k}\) 可以計算為

\[ r_k = \frac{P_k(a_0,a_1,\dots,a_k)}{P_{k-1}(a_1,\dots,a_k)} = \frac{a_k p_{k-1} + p_{k-2}}{a_k q_{k-1} + q_{k-2}} \]

其中 \(P_k(a_0,\dots,a_k)\) 是連行列式（continuant），定義為

\[ P_k(x_0,x_1,\dots,x_k) = \det \begin{bmatrix} x_k & 1 & 0 & \dots & 0 \\ -1 & x_{k-1} & 1 & \dots & 0 \\ 0 & -1 & x_2 & . & \vdots \\ \vdots & \vdots & . & \ddots & 1 \\ 0 & 0 & \dots & -1 & x_0 \end{bmatrix}_{\textstyle .} \]

因此，\(r_k\) 是 \(r_{k-1}\) 和 \(r_{k-2}\) 的加權中間值（mediant）。

為了一致性，定義了兩個額外的漸近分數 \(r_{-1} = \frac{1}{0}\) 和 \(r_{-2} = \frac{0}{1}\)。

詳細説明

漸進分數 \(r_k\) 的分子和分母可以看作 \(a_0, a_1, \dots, a_k\) 的多元多項式：

\[ r_k = \frac{P_k(a_0, a_1, \dots, a_k)}{Q_k(a_0,a_1, \dots, a_k)} \]

根據漸進分數的定義，

\[ r_k = a_0 + \frac{1}{[a_1;a_2,\dots, a_k]}= a_0 + \frac{Q_{k-1}(a_1, \dots, a_k)}{P_{k-1}(a_1, \dots, a_k)} = \frac{a_0 P_{k-1}(a_1, \dots, a_k) + Q_{k-1}(a_1, \dots, a_k)}{P_{k-1}(a_1, \dots, a_k)} \]

由此得出 \(Q_k(a_0, \dots, a_k) = P_{k-1}(a_1, \dots, a_k)\)。這產生了關係

\[ P_k(a_0, \dots, a_k) = a_0 P_{k-1}(a_1, \dots, a_k) + P_{k-2}(a_2, \dots, a_k) \]

最初，\(r_0 = \frac{a_0}{1}\) 和 \(r_1 = \frac{a_0 a_1 + 1}{a_1}\)，因此

\[ \begin{align}P_0(a_0)&=a_0,\\ P_1(a_0, a_1) &= a_0 a_1 + 1.\end{align} \]

為了保持一致性，可以方便地定義 \(P_{-1} = 1\) 和 \(P_{-2}=0\)，並正式表示 \(r_{-1} = \frac{1}{0}\) 和 \(r_{-2}=\frac{0}{1}\)。

從數值分析可知，任意三對角矩陣的行列式

\[ T_k = \det \begin{bmatrix} a_0 & b_0 & 0 & \dots & 0 \\ c_0 & a_1 & b_1 & \dots & 0 \\ 0 & c_1 & a_2 & . & \vdots \\ \vdots & \vdots & . & \ddots & c_{k-1} \\ 0 & 0 & \dots & b_{k-1} & a_k \end{bmatrix} \]

可以遞歸地計算為 \(T_k = a_k T_{k-1} - b_{k-1} c_{k-1} T_{k-2}\)。將其與 \(P_k\) 進行比較，得到一個直接表達式

\[ P_k = \det \begin{bmatrix} x_k & 1 & 0 & \dots & 0 \\ -1 & x_{k-1} & 1 & \dots & 0 \\ 0 & -1 & x_2 & . & \vdots \\ \vdots & \vdots & . & \ddots & 1 \\ 0 & 0 & \dots & -1 & x_0 \end{bmatrix}_{\textstyle .} \]

這個多項式也被稱為連行列式（continuant），由於其與連分數的密切關係。如果主對角線上的順序顛倒，則連行列式不會改變。這產生了一個計算公式：

\[ P_k(a_0, \dots, a_k) = a_k P_{k-1}(a_0, \dots, a_{k-1}) + P_{k-2}(a_0, \dots, a_{k-2}) \]

實現

把漸進分數計算為一對序列 \(p_{-2}, p_{-1}, p_0, p_1, \dots, p_k\) 和 \(q_{-2}, q_{-1}, q_0, q_1, \dots, q_k\)：

1
2
3
4
5
6
7
8
9
auto convergents(vector<int> a) {
  vector<int> p = {0, 1};
  vector<int> q = {1, 0};
  for (auto it : a) {
    p.push_back(p[p.size() - 1] * it + p[p.size() - 2]);
    q.push_back(q[q.size() - 1] * it + q[q.size() - 2]);
  }
  return make_pair(p, q);
}

1
2
3
4
5
6
7
def convergents(a):
    p = [0, 1]
    q = [1, 0]
    for it in a:
        p.append(p[-1]*it + p[-2])
        q.append(q[-1]*it + q[-2])
    return p, q

誤差和餘項的估計

誤差

漸進分數 \(r_k = \frac{p_k}{q_k}\) 與 \(r\) 的誤差（deviation）通常可以估計為

\[ \left|\frac{p_k}{q_k}-r\right| \leq \frac{1}{q_k q_{k+1}} \leq \frac{1}{q_k^2} \]

將兩邊乘以 \(q_k\)，得到另一個估計：

\[ |p_k - q_k r| \leq \frac{1}{q_{k+1}} \]

從上面的循環可以看出，\(q_k\) 的增長速度至少與斐波那契數一樣快。

在下圖中可以看到收斂 \(r_k\) 漸進分數 \(r=\frac{1+\sqrt 5}{2}\) 的可視化：

無理數 \(r=\frac{1+\sqrt 5}{2}\) 由藍色虛線表示。奇數漸進分數從上面接近它，偶數漸進分數從下面接近它。

實數 x 也可以寫成：

\[ x=[a_0,a_1,\ldots,a_k,r_{k+1}] \]

最後一項漸近分數就是 x 本身。於是根據漸進分數的遞推式，就有：

\[ x=\frac{r_{k+1}p_k+p_{k-1}}{r_{k+1}q_k+q_{k-1}} \]

於是可以估計漸進分數的誤差：

\[ x-\frac{p_k}{q_k}=\frac{(-1)^k}{q_k\left(r_{k+1}q_k+q_{k-1}\right)} \]

分別對 k 取奇數偶數就得到，x 總小於其奇數階漸近分數，大於其偶數階漸近分數。

對於數字 \(r\) 及其第 \(k\) 個漸進分數 \(r_k=\frac{p_k}{q_k}\)，以下公式成立：

\[ r_k = a_0 + \sum\limits_{i=1}^k \frac{(-1)^{i-1}}{q_i q_{i-1}} \]

特別地，這意味着

\[ r_k - r_{k-1} = \frac{(-1)^{k-1}}{q_k q_{k-1}} \]

並且

\[ p_k q_{k-1} - p_{k-1} q_k = (-1)^{k-1} \]

由此可以得出結論

\[ \left| r-\frac{p_k}{q_k} \right| \leq \frac{1}{q_{k+1}q_k} \leq \frac{1}{q_k^2} \]

後一種不等式是由於 \(r_k\) 和 \(r_{k+1}\) 通常位於 \(r\) 的不同側面，因此

\[ |r-r_k| = |r_k-r_{k+1}|-|r-r_{k+1}| \leq |r_k - r_{k+1}| \]
詳細説明

為了估計 \(|r-r_k|\)，首先估計相鄰漸進分數之間的差異。根據定義，

\[ \frac{p_k}{q_k} - \frac{p_{k-1}}{q_{k-1}} = \frac{p_k q_{k-1} - p_{k-1} q_k}{q_k q_{k-1}} \]

將分子中的 \(p_k\) 和 \(q_k\) 替換為它們的循環，得到

\[ \begin{align} p_k q_{k-1} - p_{k-1} q_k &= (a_k p_{k-1} + p_{k-2}) q_{k-1} - p_{k-1} (a_k q_{k-1} + q_{k-2}) \\&= p_{k-2} q_{k-1} - p_{k-1} q_{k-2},\end{align} \]

因此，\(r_k - r_{k-1}\) 的分子總是與 \(r_{k-1} - r_{k-2}\) 相反。反過來，它等於

\[ r_1 - r_0=\left(a_0+\frac{1}{a_1}\right)-a_0=\frac{1}{a_1} \]

因此

\[ r_k - r_{k-1} = \frac{(-1)^{k-1}}{q_k q_{k-1}} \]

這產生了 \(r_k\) 作為無限級數的部分和的替代表示：

\[ r_k = (r_k - r_{k-1}) + \dots + (r_1 - r_0) + r_0 = a_0 + \sum\limits_{i=1}^k \frac{(-1)^{i-1}}{q_i q_{i-1}} \]

根據遞歸關係，\(q_k\) 單調增加的速度至少與斐波那契數一樣快，因此

\[ r = \lim\limits_{k \to \infty} r_k = a_0 + \sum\limits_{i=1}^\infty \frac{(-1)^{i-1}}{q_i q_{i-1}} \]

總是定義明確的，因為基礎系列總是收斂的。值得注意的是，剩餘系列

\[ r-r_k = \sum\limits_{i=k+1}^\infty \frac{(-1)^{i-1}}{q_i q_{i-1}} \]

由於 \(q_i q_{i-1}\) 下降的速度，與 \((-1)^k\) 具有相同的符號。因此，偶數索引的 \(r_k\) 從下面接近 \(r\)，而奇數索引的 \(r_k\) 從上面接近：

從這張圖中可以看到

\[ |r-r_k| = |r_k - r_{k+1}| - |r-r_{k+1}| \leq |r_k - r_{k+1}| \]

因此，\(r\) 和 \(r_k\) 之間的距離永遠不會大於 \(r_k\) 和 \(r_{k+1}\) 之間的間距：

\[ \left|r-\frac{p_k}{q_k}\right| \leq \frac{1}{q_k q_{k+1}} \leq \frac{1}{q_k^2} \]

例題 擴展歐幾里得

您將獲得 \(A, B, C \in \mathbb Z\)。查找 \(x, y \in \mathbb Z\)，使 \(Ax + By = C\).

解答

雖然這個問題通常是用擴展歐幾里得算法解決的，但有一個簡單而直接的連分數的解決方案。

設 \(\frac{A}{B}=[a_0; a_1, \dots, a_k]\)。上面證明了 \(p_k q_{k-1} - p_{k-1} q_k = (-1)^{k-1}\)。將 \(p_k\) 和 \(q_k\) 替換為 \(A\) 和 \(B\)，得到

\[ Aq_{k-1} - Bp_{k-1} = (-1)^{k-1} g \]

其中 \(g = \gcd(A, B)\)。如果 \(C\) 可被 \(g\) 整除，則解為 \(x = (-1)^{k-1}\frac{C}{g} q_{k-1}\) 和 \(y = (-1)^{k}\frac{C}{g} p_{k-1}\)。

1
2
3
4
5
6
7
# return (x, y) such that Ax+By=C
# assumes that such (x, y) exists
def dio(A, B, C):
    p, q = convergents(fraction(A, B))
    C //= A // p[-1] # divide by gcd(A, B)
    t = (-1) if len(p) % 2 else 1
    return t*C*q[-2], -t*C*p[-2]

幾何解釋

格點

考慮線 \(y=rx\) 上方和下方的點的凸包。

奇數漸進分數 \((q_k;p_k)\) 是上殼的頂點，而偶數漸進分數 \((q_k;p_k)\) 則是下殼的頂點。

外殼上的所有整數頂點都作為 \((q;p)\) 獲得，這樣

\[ \frac{p}{q} = \frac{tp_{k-1} + p_{k-2}}{tq_{k-1} + q_{k-2}} \]

對於整數 \(0 \leq t \leq a_k\)。換句話説，外殼上的格點集對應於中間分數。

在下圖中可以看到 \(r=\frac{9}{7}\) 的漸進分數和中間分數（灰點）。

對於漸進分數 \(r_k = \frac{p_k}{q_k}\)，設 \(\vec r_k = (q_k;p_k)\)。然後，以下重複出現：

\[ \vec r_k = a_k \vec r_{k-1} + \vec r_{k-2} \]

設 \(\vec r = (1;r)\)。然後，每個向量 \((x;y)\) 對應於等於其斜率係數 \(\frac{y}{x}\) 的數字。

利用外積 \((x_1;y_1) \times (x_2;y_2) = x_1 y_2 - x_2 y_1\) 的概念，可以看出（參見下面的解釋）

\[ s_k = -\frac{\vec r_{k-2} \times \vec r}{\vec r_{k-1} \times \vec r} = \left|\frac{\vec r_{k-2} \times \vec r}{\vec r_{k-1} \times \vec r}\right| \]

最後一個等式是由於 \(r_{k-1}\) 和 \(r_{k-2}\) 位於 \(r\) 的不同側，因此 \(\vec r_{k-1}\) 和 \(\vec r_{k-2}\) 與 \(\vec r\) 的外積具有不同的符號。考慮到 \(a_k = \lfloor s_k \rfloor\)，\(\vec r_k\) 的公式如下

\[ \vec r_k = \vec r_{k-2} + \left\lfloor \left| \frac{\vec r \times \vec r_{k-2}}{\vec r \times \vec r_{k-1}}\right|\right\rfloor \vec r_{k-1} \]

注意到 \(\vec r_k \times r = (q;p) \times (1;r) = qr - p\)，因此

\[ a_k = \left\lfloor \left| \frac{q_{k-1}r-p_{k-1}}{q_{k-2}r-p_{k-2}} \right| \right\rfloor \]
解釋

正如已經注意到的，\(a_k = \lfloor s_k \rfloor\)，其中 \(s_k = [a_k; a_{k+1}, a_{k+2}, \dots]\)。另一方面，從漸進分數的遞推，可以得出

\[ r = [a_0; a_1, \dots, a_{k-1}, s_k] = \frac{s_k p_{k-1} + p_{k-2}}{s_k q_{k-1} + q_{k-2}} \]

在向量形式中，它重寫為

\[ \vec r \parallel s_k \vec r_{k-1} + \vec r_{k-2} \]

這意味着 \(\vec r\) 和 \(s_k \vec r_{k-1} + \vec r_{k-2}\) 共線（即具有相同的斜率係數）。用 \(\vec r\) 計算兩個部分的外積，可以得到

\[ 0 = s_k (\vec r_{k-1} \times \vec r) + (\vec r_{k-2} \times \vec r) \]

得出最終公式

\[ s_k = -\frac{\vec r_{k-2} \times \vec r}{\vec r_{k-1} \times \vec r} \]

例題 鼻子拉伸算法

每次將 \(\vec r_{k-1}\) 添加到向量 \(\vec p\) 時，\(\vec p \times \vec r\) 的值都會增加 \(\vec r_{k-1} \times \vec r\)。

因此，\(a_k=\lfloor s_k \rfloor\) 是 \(\vec r_{k-1}\) 向量的最大整數，可以將其添加到 \(\vec r_{k-2}\)，而無需更改與 \(\vec r\) 的外積的符號。

換句話説，\(a_k\) 是您可以將 \(\vec r_{k-1}\) 添加到 \(\vec r_{k-2}\) 的最大整數次數，而無需跨越 \(\vec r\) 定義的線：

在上面的圖片中，\(\vec r_2 = (4;3)\) 是通過將 \(\vec r_1 = (1;1)\) 重複添加到 \(\vec r_0 = (1;0)\) 而獲得的。

當不可能在不跨越 \(y=rx\) 線的情況下將 \(\vec r_1\) 進一步添加到 \(\vec r_0\) 時，轉到另一側，重複將 \(\vec r_2\) 添加到 \(\vec r_1\) 以獲得 \(\vec r_3 = (9;7)\)。

此過程生成接近直線的指數較長的向量。

對於這一特性，Boris Delaunay 將生成結果收斂向量的過程稱為 鼻子拉伸算法（Nose stretching algorithm）。

如果觀察在點 \(\vec r_{k-2}\)、\(\vec r_{k}\) 和 \(\vec 0\) 上繪製的三角形，會注意到它的加倍面積是

\[ |\vec r_{k-2} \times \vec r_k| = |\vec r_{k-2} \times (\vec r_{k-2} + a_k \vec r_{k-1})| = a_k |\vec r_{k-2} \times \vec r_{k-1}| = a_k \]

結合 Pick 定理，這意味着三角形內部沒有嚴格的格點，其邊界上的唯一格點是 \(\vec 0\) 和 \(\vec r_{k-2} + t \cdot \vec r_{k-1}\)，對於所有整數 \(t\)，使得 \(0 \leq t \leq a_k\)。當連接所有可能的 \(k\) 時，這意味着在由偶數索引和奇數索引收斂向量形成的多邊形之間的空間中沒有整數點。

這反過來意味着，具有奇數係數的 \(\vec r_k\) 形成了線 \(y=rx\) 上方 \(x \geq 0\) 的格點凸包，而具有偶數係數的 \(\vec r_k\) 形成線 \(y=rx\) 下方 \(x > 0\) 的格點凸包。

定義

這些多邊形也被稱為 克萊因多邊形（Klein polygons），以費利克斯·克萊因（Felix Klein）的名字命名，他首次提出了對連續分數的幾何解釋。

例題

既然已經介紹了最重要的事實和概念，那麼是時候深入研究具體的例題了。

線下凸包

找到格點 \((x;y)\) 的凸包，使得 \(r=[a_0;a_1,\dots,a_k]=\frac{p_k}{q_k}\) 的 \(0 \leq x \leq N\) 和 \(0 \leq y \leq rx\)。

解答

如果我們考慮無界集合 \(0 \leq x\)，則上凸包將由線 \(y=rx\) 本身給出。

然而，在附加約束 \(x \leq N\) 的情況下，最終需要偏離直線以保持適當的凸包。

設 \(t = \lfloor \frac{N}{q_k}\rfloor\)，則對於整數 \(1 \leq \alpha \leq t\)，在 \((0;0)\) 之後的外殼上的第一個 \(t\) 格點是 \(\alpha \cdot (q_k; p_k)\)。

然而，\((t+1)(q_k; p_k)\) 不能是下一個格點，因為 \((t+1)q_k\) 大於 \(N\)。

為了到達外殼中的下一個格點，應該到達點 \((x;y)\)，該點與 \(y=rx\) 相差最小，同時保持 \(x \leq N\)。

設 \((x; y)\) 為凸包中的最後一個當前點。然後，下一點 \((x'; y')\) 是這樣的：\(x' \leq N\) 和 \((x'; y') - (x; y) = (\Delta x; \Delta y)\) 儘可能接近線 \(y=rx\)。換句話説，\((\Delta x; \Delta y)\) 根據 \(\Delta x \leq N - x\) 和 \(\Delta y \leq r \Delta x\) 最大化 \(r \Delta x - \Delta y\)。

這樣的點位於 \(y=rx\) 以下的格點的凸包上。換句話説，\((\Delta x; \Delta y)\) 必須是 \(r\) 的下中間分數。

也就是説，對於某些奇數 \(i\) 和 \(0 \leq t < a_i\)，\((\Delta x; \Delta y)\) 的形式為 \((q_{i-1}; p_{i-1}) + t \cdot (q_i; p_i)\)。

要找到這樣的 \(i\)，可以遍歷所有可能的 \(i\)，從最大的一個開始，並對 \(i\) 使用 \(t = \lfloor \frac{N-x-q_{i-1}}{q_i} \rfloor\)，這樣 \(N-x-q_{i-1} \geq 0\)。

當 \((\Delta x; \Delta y) = (q_{i-1}; p_{i-1}) + t \cdot (q_i; p_i)\) 時，條件 \(\Delta y \leq r \Delta x\) 由中間分數的性質保持。

並且 \(t < a_i\) 成立，因為已經耗盡了從 \(i+2\) 獲得的半收斂，因此 \(x + q_{i-1} + a_i q_i = x+q_{i+1}\) 大於 \(N\)。

現在，可以將 \((\Delta x; \Delta y)\) 添加到 \((x;y)\) 中 \(k = \lfloor \frac{N-x}{\Delta x} \rfloor\) 次，然後再超過 \(N\)，之後將嘗試下一個中間分數。

 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
// returns [ah, ph, qh] such that points r[i]=(ph[i], qh[i]) constitute upper
// convex hull of lattice points on 0 <= x <= N and 0 <= y <= r * x, where r =
// [a0; a1, a2, ...] and there are ah[i]-1 integer points on the segment between
// r[i] and r[i+1]
auto hull(auto a, int N) {
  auto [p, q] = convergents(a);
  int t = N / q.back();
  vector ah = {t};
  vector ph = {0, t * p.back()};
  vector qh = {0, t * q.back()};

  for (int i = q.size() - 1; i >= 0; i--) {
    if (i % 2) {
      while (qh.back() + q[i - 1] <= N) {
        t = (N - qh.back() - q[i - 1]) / q[i];
        int dp = p[i - 1] + t * p[i];
        int dq = q[i - 1] + t * q[i];
        int k = (N - qh.back()) / dq;
        ah.push_back(k);
        ph.push_back(ph.back() + k * dp);
        qh.push_back(qh.back() + k * dq);
      }
    }
  }
  return make_tuple(ah, ph, qh);
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
# returns [ah, ph, qh] such that points r[i]=(ph[i], qh[i]) constitute upper convex hull
# of lattice points on 0 <= x <= N and 0 <= y <= r * x, where r = [a0; a1, a2, ...]
# and there are ah[i]-1 integer points on the segment between r[i] and r[i+1]
def hull(a, N):
    p, q = convergents(a)
    t = N // q[-1]
    ah = [t]
    ph = [0, t*p[-1]]
    qh = [0, t*q[-1]]
    for i in reversed(range(len(q))):
        if i % 2 == 1:
            while qh[-1] + q[i-1] <= N:
                t = (N - qh[-1] - q[i-1]) // q[i]
                dp = p[i-1] + t*p[i]
                dq = q[i-1] + t*q[i]
                k = (N - qh[-1]) // dq
                ah.append(k)
                ph.append(ph[-1] + k * dp)
                qh.append(qh[-1] + k * dq)
    return ah, ph, qh

Timus - Crime and Punishment

您將得到整數 \(A\)、\(B\) 和 \(N\)。查找 \(x \geq 0\) 和 \(y \geq 0\)，使 \(Ax + By \leq N\) 和 \(Ax + By\) 達到最大值。

解答

在這個問題中有 \(1 \leq A, B, N \leq 2 \cdot 10^9\)，因此可以用 \(O(\sqrt N)\) 來解決。但是，有一個 \(O(\log N)\) 解決方案包含連分數。

為了方便起見，通過替換 \(x \mapsto \lfloor \frac{N}{A}\rfloor - x\) 來反轉 \(x\) 的方向，因此需要找到點 \((x; y)\)，使得 \(0 \leq x \leq \lfloor \frac{N}{A} \rfloor\)、\(By - Ax \leq N \;\bmod\; A\) 和 \(By - Ax\) 是可能的最大值。每個 \(x\) 的最佳 \(y\) 值為 \(\lfloor \frac{Ax + (N \bmod A)}{B} \rfloor\)。

為了更一般地對待它，編寫一個函數，該函數在 \(0 \leq x \leq N\) 和 \(y = \lfloor \frac{Ax+B}{C} \rfloor\) 上找到最佳點。

這個問題的核心解決方案思想基本上重複了前面的問題，但不是使用下中間分數來偏離直線，而是使用上中間分數來接近直線，而不跨越直線，也不違反 \(x \leq N\)。不幸的是，與前一個問題不同，您需要確保在靠近 \(y=\frac{Ax+B}{C}\) 線時不會越過該線，因此在計算中間分數的係數 \(t\) 時應牢記這一點。

 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
# (x, y) such that y = (A*x+B) // C,
# Cy - Ax is max and 0 <= x <= N.
def closest(A, B, C, N):
    # y <= (A*x + B)/C <=> diff(x, y) <= B
    def diff(x, y):
        return C*y-A*x
    a = fraction(A, C)
    p, q = convergents(a)
    ph = [B // C]
    qh = [0]
    for i in range(2, len(q) - 1):
        if i % 2 == 0:
            while diff(qh[-1] + q[i+1], ph[-1] + p[i+1]) <= B:
                t = 1 + (diff(qh[-1] + q[i-1], ph[-1] + p[i-1]) - B - 1) // abs(diff(q[i], p[i]))
                dp = p[i-1] + t*p[i]
                dq = q[i-1] + t*q[i]
                k = (N - qh[-1]) // dq
                if k == 0:
                    return qh[-1], ph[-1]
                if diff(dq, dp) != 0:
                    k = min(k, (B - diff(qh[-1], ph[-1])) // diff(dq, dp))
                qh.append(qh[-1] + k*dq)
                ph.append(ph[-1] + k*dp)
    return qh[-1], ph[-1]

def solve(A, B, N):
    x, y = closest(A, N % A, B, N // A)
    return N // A - x, y

June Challenge 2017 - Euler Sum

計算 \(\sum\limits_{x=1}^N \lfloor \mathrm{e}x \rfloor\)，其中 \(\mathrm{e} = [2; 1, 2, 1, 1, 4, 1, 1, 6, 1, \dots, 1, 2n, 1, \dots]\) 是自然對數的底，\(N \leq 10^{4000}\)。

解答

此和等於格點 \((x;y)\) 的數量，使得 \(1 \leq x \leq N\) 和 \(1 \leq y \leq \mathrm{e}x\)。

在構造了 \(y=\mathrm{e}x\) 以下的點的凸包之後，可以使用 Pick 定理計算這個數：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
// sum floor(k * x) for k in [1, N] and x = [a0; a1, a2, ...]
int sum_floor(auto a, int N) {
  N++;
  auto [ah, ph, qh] = hull(a, N);

  // The number of lattice points within a vertical right trapezoid
  // on points (0; 0) - (0; y1) - (dx; y2) - (dx; 0) that has
  // a+1 integer points on the segment (0; y1) - (dx; y2).
  auto picks = [](int y1, int y2, int dx, int a) {
    int b = y1 + y2 + a + dx;
    int A = (y1 + y2) * dx;
    return (A - b + 2) / 2 + b - (y2 + 1);
  };

  int ans = 0;
  for (size_t i = 1; i < qh.size(); i++) {
    ans += picks(ph[i - 1], ph[i], qh[i] - qh[i - 1], ah[i - 1]);
  }
  return ans - N;
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
# sum floor(k * x) for k in [1, N] and x = [a0; a1, a2, ...]
def sum_floor(a, N):
    N += 1
    ah, ph, qh = hull(a, N)

    # The number of lattice points within a vertical right trapezoid
    # on points (0; 0) - (0; y1) - (dx; y2) - (dx; 0) that has
    # a+1 integer points on the segment (0; y1) - (dx; y2).
    def picks(y1, y2, dx, a):
        b = y1 + y2 + a + dx
        A = (y1 + y2) * dx
        return (A - b + 2) // 2 + b - (y2 + 1)

    ans = 0
    for i in range(1, len(qh)):
        ans += picks(ph[i-1], ph[i], qh[i]-qh[i-1], ah[i-1])
    return ans - N

NAIPC 2019 - It's a Mod, Mod, Mod, Mod World

給定 \(p\)、\(q\) 和 \(n\)，計算 \(\sum\limits_{i=1}^n [p \cdot i \bmod q]\)。

解答

如果您注意到 \(a \bmod b = a - \lfloor \frac{a}{b} \rfloor b\)，則此問題會減少到上一個問題。有了這個事實，總數減少到

\[ \sum\limits_{i=1}^n \left(p \cdot i - \left\lfloor \frac{p \cdot i}{q} \right\rfloor q\right) = \frac{pn(n+1)}{2}-q\sum\limits_{i=1}^n \left\lfloor \frac{p \cdot i}{q}\right\rfloor \]

然而，將 \(x\) 從 \(1\) 到 \(N\) 的 \(\lfloor rx \rfloor\) 相加，是我們能夠從上一個問題中得出的結果。

1
2
3
void solve(int p, int q, int N) {
  cout << p * N * (N + 1) / 2 - q * sum_floor(fraction(p, q), N) << "\n";
}

1
2
def solve(p, q, N):
    return p * N * (N + 1) // 2 - q * sum_floor(fraction(p, q), N)

Library Checker - Sum of Floor of Linear

給定 \(N\)、\(M\)、\(A\) 和 \(B\)，計算 \(\sum\limits_{i=0}^{N-1} \lfloor \frac{A \cdot i + B}{M} \rfloor\)。

解答

這是迄今為止技術上最麻煩的問題。

可以使用相同的方法來構造線 \(y = \frac{Ax+B}{M}\) 以下的點的全凸包。

已經知道如何解決 \(B = 0\) 的問題。此外，已經知道如何構造這個凸包，直到 \([0, N-1]\) 段上的這條線的最近格點（這在上面的「罪與罰」問題中完成）。

現在應該注意到，一旦到達了離直線最近的點，就可以假設直線實際上通過了最近的點。因為在實際直線和稍微向下移動以通過最近點的直線之間，\([0, N-1]\) 上沒有其他格點。

也就是説，要在 \([0, N-1]\) 上的線 \(y=\frac{Ax+B}{M}\) 下方構造全凸包，可以將其構造到與 \([0, N-1]\) 的線最近的點，然後繼續，就像該線通過該點一樣，重用用於構造 \(B=0\) 的凸包的算法：

 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
# hull of lattice (x, y) such that C*y <= A*x+B
def hull(A, B, C, N):
    def diff(x, y):
        return C*y-A*x
    a = fraction(A, C)
    p, q = convergents(a)
    ah = []
    ph = [B // C]
    qh = [0]

    def insert(dq, dp):
        k = (N - qh[-1]) // dq
        if diff(dq, dp) > 0:
            k = min(k, (B - diff(qh[-1], ph[-1])) // diff(dq, dp))
        ah.append(k)
        qh.append(qh[-1] + k*dq)
        ph.append(ph[-1] + k*dp)

    for i in range(1, len(q) - 1):
        if i % 2 == 0:
            while diff(qh[-1] + q[i+1], ph[-1] + p[i+1]) <= B:
                t = (B - diff(qh[-1] + q[i+1], ph[-1] + p[i+1])) // abs(diff(q[i], p[i]))
                dp = p[i+1] - t*p[i]
                dq = q[i+1] - t*q[i]
                if dq < 0 or qh[-1] + dq > N:
                    break
                insert(dq, dp)

    insert(q[-1], p[-1])

    for i in reversed(range(len(q))):
        if i % 2 == 1:
            while qh[-1] + q[i-1] <= N:
                t = (N - qh[-1] - q[i-1]) // q[i]
                dp = p[i-1] + t*p[i]
                dq = q[i-1] + t*q[i]
                insert(dq, dp)
    return ah, ph, qh

OKC 2 - From Modular to Rational

有一個有理數 \(\frac{p}{q}\)，即 \(1 \leq p, q \leq 10^9\)。您可以詢問幾個素數 \(m\) 的 \(p q^{-1}\) 模 \(m \sim 10^9\) 的值。恢復 \(\frac{p}{q}\)。

這個問題等價於：查找 \(1 \leq x \leq N\) 中，使 \(Ax \bmod M\) 最小的 \(x\)。

解答

根據中國剩餘定理，要求結果模化幾個素數與要求其模化其乘積是相同的。因此，在不喪失一般性的情況下，假設知道餘數模足夠大的數 \(m\)。

對於給定的餘數 \(r\)，可能有幾種可能的解決方案 \((p, q)\) 到 \(p \equiv qr \pmod m\)。然而，如果 \((p_1, q_1)\) 和 \((p_2, q_2)\) 都是解，那麼它也認為 \(p_1 q_2 \equiv p_2 q_1 \pmod m\)。假設 \(\frac{p_1}{q_1} \neq \frac{p_2}{q_2}\)，則意味着 \(|p_1 q_2 - p_2 q_1|\) 至少為 \(m\)。

題面有 \(1 \leq p, q \leq 10^9\)，因此，如果 \(p_1, q_1\) 和 \(p_2, q_2\) 最多都是 \(10^9\) 的話，那麼差額最多為 \(10^{18}\)。對於 \(m > 10^{18}\)，這意味着具有 \(\frac{p}{q}\) 的解 \(1 \leq p, q \leq 10^9\) 作為有理數是唯一的。

因此，問題歸結為，給定 \(r\) 模 \(m\)，找到任何 \(q\)，使得 \(1 \leq q \leq 10^9\) 和 \(qr \;\bmod\; m \leq 10^9\)。

這實際上與找到 \(1 \leq q \leq 10^9\) 的 \(q\) 是相同的，該 \(q\) 提供了可能的最小 \(qr \bmod m\)。

對於 \(qr = km + b\)，這意味着需要找到一對 \((q, m)\)，使得 \(1 \leq q \leq 10^9\) 和 \(qr - km \geq 0\) 是可能的最小值。

由於 \(m\) 是常量，可以除以它，並進一步將其重新表述為求解 \(q\)，這樣 \(1 \leq q \leq 10^9\) 和 \(\frac{r}{m} q - k \geq 0\) 是可能的最小值。

就連分數而言，這意味着 \(\frac{k}{q}\) 是 \(\frac{r}{m}\) 的最佳丟番圖近似值，並且僅檢查 \(\frac{r}{m}\) 的下中間分數就足夠了。

1
2
3
4
5
6
7
8
# find Q that minimizes Q*r mod m for 1 <= k <= n < m 
def mod_min(r, n, m):
    a = fraction(r, m)
    p, q = convergents(a)
    for i in range(2, len(q)):
        if i % 2 == 1 and (i + 1 == len(q) or q[i+1] > n):
            t = (n - q[i-1]) // q[i]
            return q[i-1] + t*q[i]

習題

本頁面主要譯自博文 Continued fractions，版權協議為 CC-BY-SA 4.0。

本页面最近更新：，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用