콘텐츠로 이동

연분수

연분수는 실수를 유리수로 이루어진 특정한 수렴 수열로 나타낸 표현입니다. 연분수는 계산하기 쉽고, (분모가 주어진 값을 넘지 않는 모든 수 중에서) 원래 실수에 대한 최선의 유리 근사를 효율적으로 찾는 데 사용할 수 있으므로 경쟁 프로그래밍에서 유용합니다.

그 밖에도 연분수는 유클리드 알고리즘과 밀접한 관련이 있어 여러 정수론 문제에서 유용합니다.

연분수 표현

정의

$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]$로 나타냅니다.

Example

$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} $$

모든 유리수는 정확히 $2$가지 방법으로 연분수로 나타낼 수 있음을 증명할 수 있습니다.

$$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,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_0 = 0$일 때 $\frac{1}{r} = [a_1; a_2, \dots]$라는 것입니다.

정의

위의 정의에서 유리수 $r_0, r_1, r_2, \dots$를 $r$의 근사분수라고 합니다.

이에 따라 개별 $r_k = [a_0; a_1, \dots, a_k] = \frac{p_k}{q_k}$를 $r$의 $k$번째 근사분수라고 합니다.

Example

$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}$로 정의되는 피보나치 수열입니다. 비네 공식으로부터 다음이 알려져 있습니다.

$$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]$를 반수렴자라고 합니다.

일반적으로 $r$보다 큰 (준)수렴자를 위쪽 (준)수렴자라고 하고, $r$보다 작은 (준)수렴자를 아래쪽 (준)수렴자라고 부르겠습니다.

정의

수렴자와 짝을 이루는 개념으로, $s_k = [a_k; a_{k+1}, a_{k+2}, \dots]$를 완전몫(complete quotient)으로 정의합니다.

이에 따라 개별 $s_k$를 $r$의 $k$번째 완전몫이라고 부르겠습니다.

위의 정의로부터 $k \geq 1$에 대하여 $s_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$에 대하여 유일하게 정의됩니다.

구현

코드 예시에서는 주로 유한 연분수를 가정하겠습니다.

$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$이라는 사실도 도출됩니다. 그러므로 수렴자는 항상 기약분수입니다.

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;
}

def fraction(p, q):
    a = []
    while q:
        a.append(p // q)
        p, q = q, p % q
    return a

주요 결과

연분수를 더 깊이 공부할 동기를 부여하기 위해, 이제 몇 가지 핵심적인 사실을 소개합니다.

점화식

수렴자 $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}$입니다.

편차

$r_k = \frac{p_k}{q_k}$와 $r$의 편차는 일반적으로 다음과 같이 추정할 수 있습니다.

$$\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}$는 파란 점선으로 표시되어 있습니다. 홀수 번째 수렴자는 위쪽에서, 짝수 번째 수렴자는 아래쪽에서 이 값에 접근합니다.

격자 껍질

직선 $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}$의 수렴자와 반수렴자(중간의 회색 점)를 볼 수 있습니다.

최선 근사

어떤 $x$에 대해 $q \leq x$라는 조건 아래에서 $\left|r-\frac{p}{q}\right|$를 최소화하는 분수를 $\frac{p}{q}$라고 합시다.

그러면 $\frac{p}{q}$는 $r$의 반수렴자입니다.

이 마지막 사실을 이용하면 $r$의 반수렴자를 확인하여 $r$의 최선 유리 근사를 찾을 수 있습니다.

아래에서는 이러한 사실들에 대한 추가 설명과 함께 약간의 직관과 해석을 제공합니다.

수렴자

앞에서 정의한 수렴자를 좀 더 자세히 살펴봅시다. $r=[a_0, a_1, a_2, \dots]$에 대해, 그 수렴자는 다음과 같습니다.

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

수렴자는 연분수의 핵심 개념이므로, 그 성질을 공부하는 것이 중요합니다.

수 $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라고도 불립니다. 주대각선의 수열을 뒤집어도 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$의 쌍으로 계산합니다.

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);
}

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

연분수의 트리

가능한 모든 연분수를 유용한 트리 구조로 묶는 주요 방법은 두 가지입니다.

Stern-Brocot 트리

Stern-Brocot 트리는 서로 다른 모든 양의 유리수를 포함하는 이진 탐색 트리입니다.

이 트리는 일반적으로 다음과 같은 모습입니다.

이미지 제작: Aaron Rotenberg 라이선스: CC BY-SA 3.0

분수 $\frac{0}{1}$과 $\frac{1}{0}$은 각각 트리의 왼쪽과 오른쪽에 "가상으로" 놓여 있습니다.

그러면 노드에 있는 분수는 그 위에 있는 두 분수 $\frac{a}{b}$와 $\frac{c}{d}$의 중간수 $\frac{a+c}{b+d}$입니다.

점화식 $\frac{p_k}{q_k}=\frac{a_k p_{k-1} + p_{k-2}}{a_k q_{k-1} + q_{k-2}}$은 연분수 표현이 트리에서 $\frac{p_k}{q_k}$까지의 경로를 인코딩한다는 뜻입니다. $[a_0; a_1, \dots, a_{k}, 1]$을 찾으려면 오른쪽으로 $a_0$번, 왼쪽으로 $a_1$번, 오른쪽으로 $a_2$번 이동하는 식으로 $a_k$까지 이동해야 합니다.

그러면 $[a_0; a_1, \dots, a_k,1]$의 부모는 마지막으로 이동한 방향으로 한 단계 되돌아가서 얻는 분수입니다.

다시 말해, $a_k > 1$일 때는 $[a_0; a_1, \dots, a_k-1,1]$이고 $a_k = 1$일 때는 $[a_0; a_1, \dots, a_{k-1}, 1]$입니다.

따라서 $[a_0; a_1, \dots, a_k, 1]$의 자식은 $[a_0; a_1, \dots, a_k+1, 1]$과 $[a_0; a_1, \dots, a_k, 1, 1]$입니다.

Stern-Brocot 트리에 인덱스를 매겨 봅시다. 루트 정점에는 인덱스 $1$을 부여합니다. 그런 다음 정점 $v$에 대해, 왼쪽 자식의 인덱스는 $v$의 최상위 비트를 $1$에서 $10$으로 바꾸어 부여하고, 오른쪽 자식의 인덱스는 최상위 비트를 $1$에서 $11$로 바꾸어 부여합니다.

이 인덱싱에서 유리수의 연분수 표현은 그 이진 인덱스의 런 렝스 인코딩을 나타냅니다.

$\frac{5}{2} = [2;2] = [2;1,1]$의 경우 인덱스는 $1011_2$이고, 비트를 오름차순으로 볼 때 런 렝스 인코딩은 $[2;1,1]$입니다.

또 다른 예로 $\frac{2}{5} = [0;2,2]=[0;2,1,1]$이 있습니다. 이 수의 인덱스는 $1100_2$이고, 실제로 이 인덱스의 런 렝스 인코딩은 $[0;2,2]$입니다.

Stern-Brocot 트리는 사실 트립이라는 점에 주목할 만합니다. 즉, $\frac{p}{q}$를 기준으로는 이진 탐색 트리이지만, $p$와 $q$ 모두를 기준으로는 힙입니다.

연분수 비교하기

$A=[a_0; a_1, \dots, a_n]$와 $B=[b_0; b_1, \dots, b_m]$가 주어집니다. 어느 분수가 더 작습니까?

풀이

우선 $A$와 $B$가 무리수이고, 이들의 연분수 표현이 Stern-Brocot 트리에서의 무한한 하강을 나타낸다고 가정합시다.

앞서 언급했듯이, 이 표현에서 $a_0$은 하강 과정에서 오른쪽으로 꺾는 횟수를, $a_1$은 그 뒤에 이어서 왼쪽으로 꺾는 횟수를 나타내며, 이후도 마찬가지입니다. 따라서 $a_k$와 $b_k$를 비교할 때 $a_k = b_k$이면 그냥 $a_{k+1}$과 $b_{k+1}$의 비교로 넘어가면 됩니다. 그렇지 않은 경우, $A < B$인지 판단하려면 오른쪽 하강 중이라면 $a_k < b_k$인지를, 왼쪽 하강 중이라면 $a_k > b_k$인지를 확인해야 합니다.

다시 말해, 무리수 $A$와 $B$에 대해 $A < B$일 필요충분조건은 사전순 비교로 $(a_0, -a_1, a_2, -a_3, \dots) < (b_0, -b_1, b_2, -b_3, \dots)$인 것입니다.

이제 형식적으로 $\infty$를 연분수 표현의 원소로 사용하면 무리수 $A-\varepsilon$과 $A+\varepsilon$, 즉 $A$보다 작지만(크지만) 다른 어떤 실수보다도 큰(작은) 원소를 흉내 낼 수 있습니다. 구체적으로, $A=[a_0; a_1, \dots, a_n]$에 대해 이 두 원소 중 하나는 $[a_0; a_1, \dots, a_n, \infty]$로, 다른 하나는 $[a_0; a_1, \dots, a_n - 1, 1, \infty]$로 흉내 낼 수 있습니다.

어느 쪽이 $A-\varepsilon$에 해당하고 어느 쪽이 $A+\varepsilon$에 해당하는지는 $n$의 홀짝성으로 판단하거나, 두 수를 무리수로서 비교하여 판단할 수 있습니다.

# check if a < b assuming that a[-1] = b[-1] = infty and a != b
def less(a, b):
    a = [(-1)**i*a[i] for i in range(len(a))]
    b = [(-1)**i*b[i] for i in range(len(b))]
    return a < b

# [a0; a1, ..., ak] -> [a0, a1, ..., ak-1, 1]
def expand(a):
    if a: # empty a = inf
        a[-1] -= 1
        a.append(1)
    return a

# return a-eps, a+eps
def pm_eps(a):
    b = expand(a.copy())
    a.append(float('inf'))
    b.append(float('inf'))
    return (a, b) if less(a, b) else (b, a)

최적의 내부 점

$\frac{0}{1} \leq \frac{p_0}{q_0} < \frac{p_1}{q_1} \leq \frac{1}{0}$이 주어집니다. $(q; p)$가 사전순으로 가장 작으면서 $\frac{p_0}{q_0} < \frac{p}{q} < \frac{p_1}{q_1}$을 만족하는 유리수 $\frac{p}{q}$를 구하십시오.

풀이

Stern-Brocot 트리의 관점에서 보면, 이는 $\frac{p_0}{q_0}$와 $\frac{p_1}{q_1}$의 LCA를 구해야 한다는 뜻입니다. Stern-Brocot 트리와 연분수 사이의 관계 때문에, 이 LCA는 대략 $\frac{p_0}{q_0}$와 $\frac{p_1}{q_1}$의 연분수 표현의 최장 공통 접두사에 대응합니다.

따라서 $\frac{p_0}{q_0} = [a_0; a_1, \dots, a_{k-1}, a_k, \dots]$와 $\frac{p_1}{q_1} = [a_0; a_1, \dots, a_{k-1}, b_k, \dots]$가 무리수라면, LCA는 $[a_0; a_1, \dots, \min(a_k, b_k)+1]$입니다.

유리수 $r_0$과 $r_1$의 경우, 둘 중 하나가 LCA 자체일 수 있으므로 경우를 나누어 처리해야 합니다. 유리수 $r_0$과 $r_1$에 대한 풀이를 단순화하기 위해, 앞의 문제에서 유도한 $r_0 + \varepsilon$과 $r_1 - \varepsilon$의 연분수 표현을 사용할 수 있습니다.

# finds lexicographically smallest (q, p)
# such that p0/q0 < p/q < p1/q1
def middle(p0, q0, p1, q1):
    a0 = pm_eps(fraction(p0, q0))[1]
    a1 = pm_eps(fraction(p1, q1))[0]
    a = []
    for i in range(min(len(a0), len(a1))):
        a.append(min(a0[i], a1[i]))
        if a0[i] != a1[i]:
            break
    a[-1] += 1
    p, q = convergents(a)
    return p[-1], q[-1]

GCJ 2019, Round 2 - New Elements: Part 2

양의 정수 쌍 $(C_i, J_i)$가 $N$개 주어집니다. $C_i x + J_i y$가 순증가 수열이 되도록 하는 양의 정수 쌍 $(x, y)$를 구해야 합니다.

그러한 쌍 중에서 사전순으로 가장 작은 것을 구하십시오.

풀이

문제를 바꾸어 말하면, 모든 $i$에 대해 $A_i x + B_i y$가 양수여야 합니다. 여기서 $A_i = C_i - C_{i-1}$이고 $B_i = J_i - J_{i-1}$입니다.

$A_i x + B_i y > 0$ 꼴의 이러한 방정식들은 네 가지 중요한 그룹으로 나눌 수 있습니다:

  1. $x, y > 0$인 해를 찾고 있으므로 $A_i, B_i > 0$인 경우는 무시할 수 있습니다.
  2. $A_i, B_i \leq 0$인 경우에는 답이 "IMPOSSIBLE"이 됩니다.
  3. $A_i > 0$, $B_i \leq 0$인 경우입니다. 이러한 제약 조건은 $\frac{y}{x} < \frac{A_i}{-B_i}$와 동치입니다.
  4. $A_i \leq 0$, $B_i > 0$인 경우입니다. 이러한 제약 조건은 $\frac{y}{x} > \frac{-A_i}{B_i}$와 동치입니다.

네 번째 그룹의 $\frac{-A_i}{B_i}$ 중 가장 큰 값을 $\frac{p_0}{q_0}$, 세 번째 그룹의 $\frac{A_i}{-B_i}$ 중 가장 작은 값을 $\frac{p_1}{q_1}$이라고 합시다.

이제 문제는 $\frac{p_0}{q_0} < \frac{p_1}{q_1}$이 주어졌을 때, $(q;p)$가 사전순으로 가장 작으면서 $\frac{p_0}{q_0} < \frac{p}{q} < \frac{p_1}{q_1}$을 만족하는 분수 $\frac{p}{q}$를 찾는 것입니다.

    def solve():
    n = int(input())
    C = [0] * n
    J = [0] * n
    # p0/q0 < y/x < p1/q1
    p0, q0 = 0, 1
    p1, q1 = 1, 0
    fail = False
    for i in range(n):
        C[i], J[i] = map(int, input().split())
        if i > 0:
            A = C[i] - C[i-1]
            B = J[i] - J[i-1]
            if A <= 0 and B <= 0:
                fail = True
            elif B > 0 and A < 0: # y/x > (-A)/B if B > 0
                if (-A)*q0 > p0*B:
                    p0, q0 = -A, B
            elif B < 0 and A > 0: # y/x < A/(-B) if B < 0
                if A*q1 < p1*(-B):
                    p1, q1 = A, -B
    if p0*q1 >= p1*q0 or fail:
        return 'IMPOSSIBLE'

    p, q = middle(p0, q0, p1, q1)
    return str(q) + ' ' + str(p)

Calkin-Wilf 트리

연분수를 이진 트리로 구성하는 조금 더 간단한 방법으로 Calkin-Wilf 트리가 있습니다.

이 트리는 대체로 다음과 같은 모습입니다:

이미지 작성자: Olli Niemitalo, Proz 라이선스: CC0 1.0

트리의 루트에는 수 $\frac{1}{1}$이 위치합니다. 그리고 수 $\frac{p}{q}$를 가진 정점의 자식은 $\frac{p}{p+q}$와 $\frac{p+q}{q}$입니다.

Stern-Brocot 트리와 달리 Calkin-Wilf 트리는 이진 탐색 트리가 아니므로, 유리수 이진 탐색을 수행하는 데 사용할 수 없습니다.

Calkin-Wilf 트리에서 분수 $\frac{p}{q}$의 직접 부모는 $p>q$일 때 $\frac{p-q}{q}$이고, 그렇지 않으면 $\frac{p}{q-p}$입니다.

Stern-Brocot 트리에서는 수렴 분수에 대한 점화식을 사용했습니다. 연분수와 Calkin-Wilf 트리 사이의 관계를 파악하려면 완전몫에 대한 점화식을 떠올려야 합니다. $s_k = \frac{p}{q}$이면 $s_{k+1} = \frac{q}{p \mod q} = \frac{q}{p-\lfloor p/q \rfloor \cdot q}$입니다.

한편, $p > q$인 동안 $s_k = \frac{p}{q}$에서 Calkin-Wilf 트리의 부모로 반복해서 이동하면 결국 $\frac{p \mod q}{q} = \frac{1}{s_{k+1}}$에 도달합니다. 이 과정을 계속하면 $s_{k+2}$에 도달하고, 이어서 $\frac{1}{s_{k+3}}$에 도달하는 식으로 진행됩니다. 이로부터 다음을 유도할 수 있습니다:

  1. $a_0> 0$일 때, Calkin-Wilf 트리에서 $[a_0; a_1, \dots, a_k]$의 직접 부모는 $\frac{p-q}{q}=[a_0 - 1; a_1, \dots, a_k]$입니다.
  2. $a_0 = 0$이고 $a_1 > 1$일 때, 그 직접 부모는 $\frac{p}{q-p} = [0; a_1 - 1, a_2, \dots, a_k]$입니다.
  3. 그리고 $a_0 = 0$이고 $a_1 = 1$일 때, 직접적인 부모는 $\frac{p}{q-p} = [a_2; a_3, \dots, a_k]$입니다.

이에 대응하여, $\frac{p}{q} = [a_0; a_1, \dots, a_k]$의 자식은 다음과 같습니다.

  1. $\frac{p+q}{q}=1+\frac{p}{q}$, 즉 $[a_0+1; a_1, \dots, a_k]$,
  2. $\frac{p}{p+q} = \frac{1}{1+\frac{q}{p}}$, 즉 $a_0 > 0$일 때는 $[0, 1, a_0, a_1, \dots, a_k]$이고 $a_0=0$일 때는 $[0, a_1+1, a_2, \dots, a_k]$입니다.

주목할 만한 점은, Calkin-Wilf 트리의 정점에 너비 우선 탐색 순서로 번호를 매기면(즉, 루트의 번호가 $1$이고, 정점 $v$의 자식들의 인덱스가 각각 $2v$와 $2v+1$이면) Calkin-Wilf 트리에서 유리수의 인덱스가 Stern-Brocot 트리에서의 인덱스와 같아진다는 것입니다.

따라서 Stern-Brocot 트리와 Calkin-Wilf 트리의 같은 레벨에 있는 수들은 동일하지만, 그 순서는 비트 반전 순열에 의해 달라집니다.

수렴

수 $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}$의 분자의 부호를 반전한 것입니다. 이 값은 다시 다음의 경우 $1$과 같습니다.

$$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=\phi = \frac{1+\sqrt{5}}{2}=[1;1,1,\dots]$의 근사분수들과 $r$로부터의 거리.

이 그림으로부터 다음을 알 수 있습니다.

$$|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$가 주어집니다. $Ax + By = C$를 만족하는 $x, y \in \mathbb Z$를 찾으십시오.

풀이

이 문제는 보통 확장 유클리드 알고리즘으로 풀지만, 연분수를 이용한 간단하고 직관적인 풀이도 있습니다.

$\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}$입니다.

# 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]

일차 분수 변환

연분수와 관련된 또 다른 중요한 개념은 이른바 일차 분수 변환입니다.

정의

일차 분수 변환은 어떤 $a,b,c,d \in \mathbb R$에 대해 $f(x) = \frac{ax+b}{cx+d}$를 만족하는 함수 $f : \mathbb R \to \mathbb R$입니다.

일차 분수 변환 $L_0(x)=\frac{a_0 x + b_0}{c_0 x + d_0}$와 $L_1(x)=\frac{a_1 x + b_1}{c_1 x + d_1}$의 합성 $(L_0 \circ L_1)(x) = L_0(L_1(x))$ 역시 그 자체로 일차 분수 변환입니다:

$$\frac{a_0\frac{a_1 x + b_1}{c_1 x + d_1} + b_0}{c_0 \frac{a_1 x + b_1}{c_1 x + d_1} + d_0} = \frac{a_0(a_1 x + b_1) + b_0 (c_1 x + d_1)}{c_0 (a_1 x + b_1) + d_0 (c_1 x + d_1)} = \frac{(a_0 a_1 + b_0 c_1) x + (a_0 b_1 + b_0 d_1)}{(c_0 a_1 + d_0 c_1) x + (c_0 b_1 + d_0 d_1)}.$$

일차 분수 변환의 역변환 또한 일차 분수 변환입니다:

$$y = \frac{ax+b}{cx+d} \iff y(cx+d) = ax + b \iff x = -\frac{dy-b}{cy-a}.$$

DMOPC '19 Contest 7 P4 - Bob and Continued Fractions

양의 정수로 이루어진 배열 $a_1, \dots, a_n$이 주어집니다. $m$개의 쿼리에 답해야 합니다. 각 쿼리는 $[a_l; a_{l+1}, \dots, a_r]$을 계산하는 것입니다.

풀이

연분수를 이어 붙일 수 있다면 이 문제를 세그먼트 트리로 풀 수 있습니다.

일반적으로 $[a_0; a_1, \dots, a_k, b_0, b_1, \dots, b_k] = [a_0; a_1, \dots, a_k, [b_1; b_2, \dots, b_k]]$가 성립합니다.

$L_{k}(x) = [a_k; x] = a_k + \frac{1}{x} = \frac{a_k\cdot x+1}{1\cdot x + 0}$이라고 표기합시다. $L_k(\infty) = a_k$임에 유의하십시오. 이 표기법에서 다음이 성립합니다.

$$[a_0; a_1, \dots, a_k, x] = [a_0; [a_1; [\dots; [a_k; x]]]] = (L_0 \circ L_1 \circ \dots \circ L_k)(x) = \frac{p_k x + p_{k-1}}{q_k x + q_{k-1}}.$$

따라서 문제는 다음을 계산하는 것으로 귀결됩니다.

$$(L_l \circ L_{l+1} \circ \dots \circ L_r)(\infty).$$

변환의 합성은 결합법칙을 만족하므로, 세그먼트 트리의 각 노드에서 그 서브트리에 속한 변환들의 합성을 계산할 수 있습니다.

연분수의 일차 분수 변환

$L(x) = \frac{ax+b}{cx+d}$라고 합시다. $A=[a_0; a_1, \dots, a_n]$에 대해 $L(A)$의 연분수 표현 $[b_0; b_1, \dots, b_m]$을 계산하십시오.

이를 이용하면 임의의 $\frac{p}{q}$에 대해 $A + \frac{p}{q} = \frac{qA + p}{q}$와 $A \cdot \frac{p}{q} = \frac{p A}{q}$를 계산할 수 있습니다.

풀이

앞에서 언급했듯이 $[a_0; a_1, \dots, a_k] = (L_{a_0} \circ L_{a_1} \circ \dots \circ L_{a_k})(\infty)$이므로, $L([a_0; a_1, \dots, a_k]) = (L \circ L_{a_0} \circ L_{a_1} \circ \dots L_{a_k})(\infty)$입니다.

따라서 $L_{a_0}$, $L_{a_1}$ 등을 차례로 덧붙여 나가면 다음을 계산할 수 있습니다.

$$(L \circ L_{a_0} \circ \dots \circ L_{a_k})(x) = L\left(\frac{p_k x + p_{k-1}}{q_k x + q_{k-1}}\right)=\frac{a_k x + b_k}{c_k x + d_k}.$$

$L(x)$는 가역이므로 $x$에 대해 단조이기도 합니다. 따라서 임의의 $x \geq 0$에 대해 $L(\frac{p_k x + p_{k-1}}{q_k x + q_{k-1}})$은 $L(\frac{p_k}{q_k}) = \frac{a_k}{c_k}$와 $L(\frac{p_{k-1}}{q_{k-1}}) = \frac{b_k}{d_k}$ 사이에 있습니다.

또한 $x=[a_{k+1}; \dots, a_n]$일 때 이 값은 $L(A)$와 같습니다. 그러므로 $b_0 = \lfloor L(A) \rfloor$는 $\lfloor L(\frac{p_k}{q_k}) \rfloor$와 $\lfloor L(\frac{p_{k-1}}{q_{k-1}}) \rfloor$ 사이에 있습니다. 이 두 값이 같으면 둘 다 $b_0$와도 같습니다.

$L(A) = (L_{b_0} \circ L_{b_1} \circ \dots \circ L_{b_m})(\infty)$임에 유의하십시오. $b_0$를 알면 현재 변환에 $L_{b_0}^{-1}$을 합성하고 $L_{a_{k+1}}$, $L_{a_{k+2}}$ 등을 계속 추가하면서 새로운 바닥값들이 일치하는 때를 찾을 수 있으며, 이를 통해 $b_1$ 등을 차례로 알아내어 $[b_0; b_1, \dots, b_m]$의 모든 값을 복원할 수 있습니다.

연분수 산술

$A=[a_0; a_1, \dots, a_n]$이고 $B=[b_0; b_1, \dots, b_m]$이라고 합시다. $A+B$와 $A \cdot B$의 연분수 표현을 계산하십시오.

풀이

기본 아이디어는 앞의 문제와 비슷하지만, $L(x) = \frac{ax+b}{cx+d}$ 대신 쌍선형 분수 변환 $L(x, y) = \frac{axy+bx+cy+d}{exy+fx+gy+h}$를 고려해야 합니다.

$L(x) \mapsto L(L_{a_k}(x))$ 대신 현재 변환을 $L(x, y) \mapsto L(L_{a_k}(x), y)$ 또는 $L(x, y) \mapsto L(x, L_{b_k}(y))$로 바꿉니다.

그런 다음 $\lfloor \frac{a}{e} \rfloor = \lfloor \frac{b}{f} \rfloor = \lfloor \frac{c}{g} \rfloor = \lfloor \frac{d}{h} \rfloor$인지 확인하고, 모두 일치하면 이 값을 결과 분수의 $c_k$로 사용하며 변환을 다음과 같이 바꿉니다.

$$L(x, y) \mapsto \frac{1}{L(x, y) - c_k}.$$

정의

어떤 $k$에 대해 $x = [a_0; a_1, \dots, a_k, x]$이면 연분수 $x = [a_0; a_1, \dots]$를 주기적이라고 합니다.

$y$가 주기적일 때 $x = [a_0; a_1, \dots, a_k, y]$이면 연분수 $x = [a_0; a_1, \dots]$를 결국 주기적이라고 합니다.

$x = [1; 1, 1, \dots]$에 대해 $x = 1 + \frac{1}{x}$가 성립하므로 $x^2 = x + 1$입니다. 주기적 연분수와 이차방정식 사이에는 일반적인 연관성이 있습니다. 다음 방정식을 생각해 봅시다:

$$ x = [a_0; a_1, \dots, a_k, x].$$

한편으로 이 방정식은 $x$의 연분수 표현이 주기 $k+1$로 주기적이라는 것을 의미합니다.

다른 한편으로 근사분수 공식을 이용하면 이 방정식은 다음을 의미합니다.

$$x = \frac{p_k x + p_{k-1}}{q_k x + q_{k-1}}.$$

즉, $x$는 자기 자신의 일차 분수 변환입니다. 이 방정식으로부터 $x$가 다음 이차방정식의 근임을 알 수 있습니다:

$$q_k x^2 + (q_{k-1}-p_k)x - p_{k-1} = 0.$$

결국 주기적인 연분수, 즉 $y=[b_0; b_1, \dots, b_k, y]$에 대해 $x = [a_0; a_1, \dots, a_k, y]$인 경우에도 비슷한 논리가 성립합니다. 실제로 첫 번째 방정식에서 $x = L_0(y)$를, 두 번째 방정식에서 $y = L_1(y)$를 유도할 수 있으며, 여기서 $L_0$와 $L_1$은 일차 분수 변환입니다. 따라서,

$$x = (L_0 \circ L_1)(y) = (L_0 \circ L_1 \circ L_0^{-1})(x).$$

더 나아가 정수 계수를 갖는 임의의 이차방정식 $ax^2+bx+c=0$에 대해 그 해 $x$가 결국 주기적인 연분수라는 것도 증명할 수 있습니다(이는 Lagrange가 처음 증명했습니다).

이차 무리수

$x, y, z, n \in \mathbb Z$이고 $n > 0$이 완전제곱수가 아닐 때 $\alpha = \frac{x+y\sqrt{n}}{z}$의 연분수를 구하십시오.

풀이

이 수의 $k$번째 완전몫 $s_k$에 대해 일반적으로 다음이 성립합니다.

$$\alpha = [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}}.$$

따라서,

$$s_k = -\frac{\alpha q_{k-1} - p_{k-1}}{\alpha q_k - p_k} = -\frac{q_{k-1} y \sqrt n + (x q_{k-1} - z p_{k-1})}{q_k y \sqrt n + (xq_k-zp_k)}.$$

분자와 분모에 $(xq_k - zp_k) - q_k y \sqrt n$을 곱하면 분모에서 $\sqrt n$을 없앨 수 있으므로, 완전몫은 다음과 같은 꼴이 됩니다.

$$s_k = \frac{x_k + y_k \sqrt n}{z_k}.$$

$s_k$를 알고 있다고 가정하고 $s_{k+1}$을 구해 봅시다.

우선 $a_k = \lfloor s_k \rfloor = \left\lfloor \frac{x_k + y_k \lfloor \sqrt n \rfloor}{z_k} \right\rfloor$입니다. 그러면,

$$s_{k+1} = \frac{1}{s_k-a_k} = \frac{z_k}{(x_k - z_k a_k) + y_k \sqrt n} = \frac{z_k (x_k - y_k a_k) - y_k z_k \sqrt n}{(x_k - y_k a_k)^2 - y_k^2 n}.$$

따라서 $t_k = x_k - y_k a_k$로 나타내면 다음이 성립합니다.

$$\begin{align}x_{k+1} &=& z_k t_k, \\ y_{k+1} &=& -y_k z_k, \\ z_{k+1} &=& t_k^2 - y_k^2 n.\end{align}$$

이러한 표현의 좋은 점은 $x_{k+1}, y_{k+1}, z_{k+1}$을 이들의 최대공약수로 나누어 약분하면 그 결과가 유일하다는 것입니다. 따라서 이를 이용하여 현재 상태가 이미 반복되었는지 확인할 수 있고, 이 상태를 가졌던 이전 인덱스가 어디였는지도 확인할 수 있습니다.

아래는 $\alpha = \sqrt n$의 연분수 표현을 계산하는 코드입니다.

# compute the continued fraction of sqrt(n)
def sqrt(n):
    n0 = math.floor(math.sqrt(n))
    x, y, z = 1, 0, 1
    a = []
    def step(x, y, z):
        a.append((x * n0 + y) // z)
        t = y - a[-1]*z
        x, y, z = -z*x, z*t, t**2 - n*x**2
        g = math.gcd(x, math.gcd(y, z))
        return x // g, y // g, z // g

    used = dict()
    for i in range(n):
        used[x, y, z] = i
        x, y, z = step(x, y, z)
        if (x, y, z) in used:
            return a

같은 step 함수를 사용하되 초기 $x$, $y$, $z$를 다르게 설정하면 임의의 $\frac{x+y \sqrt{n}}{z}$에 대해서도 이를 계산할 수 있습니다.

Tavrida NU Akai Contest - Continued Fraction

$x$와 $k$가 주어지며, $x$는 완전제곱수가 아닙니다. $\sqrt x = [a_0; a_1, \dots]$라 할 때, $0 \leq k \leq 10^9$에 대해 $\frac{p_k}{q_k}=[a_0; a_1, \dots, a_k]$를 구하십시오.

풀이

$\sqrt x$의 주기를 계산한 후에는, 연분수 표현으로부터 유도되는 일차 분수 변환에 이진 거듭제곱을 적용하여 $a_k$를 계산할 수 있습니다. 최종 변환을 구하려면 크기가 $T$인 주기를 하나의 변환으로 압축하여 $\lfloor \frac{k-1}{T}\rfloor$번 반복한 다음, 나머지 변환들과 직접 결합하면 됩니다.

x, k = map(int, input().split())

mod = 10**9+7

# compose (A[0]*x + A[1]) / (A[2]*x + A[3]) and (B[0]*x + B[1]) / (B[2]*x + B[3])
def combine(A, B):
    return [t % mod for t in [A[0]*B[0]+A[1]*B[2], A[0]*B[1]+A[1]*B[3], A[2]*B[0]+A[3]*B[2], A[2]*B[1]+A[3]*B[3]]]

A = [1, 0, 0, 1] # (x + 0) / (0*x + 1) = x

a = sqrt(x)

T = len(a) - 1 # period of a

# apply ak + 1/x = (ak*x+1)/(1x+0) to (Ax + B) / (Cx + D)
for i in reversed(range(1, len(a))):
    A = combine([a[i], 1, 1, 0], A)

def bpow(A, n):
    return [1, 0, 0, 1] if not n else combine(A, bpow(A, n-1)) if n % 2 else bpow(combine(A, A), n // 2)


C = (0, 1, 0, 0) # = 1 / 0
while k % T:
    i = k % T
    C = combine([a[i], 1, 1, 0], C)
    k -= 1

C = combine(bpow(A, k // T), C)
C = combine([a[0], 1, 1, 0], C)
print(str(C[1]) + '/' + str(C[3]))

기하학적 해석

수렴값 $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$의 의사 스칼라곱과 $\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 p$에 $\vec r_{k-1}$을 더할 때마다 $\vec p \times \vec r$의 값은 $\vec r_{k-1} \times \vec r$만큼 증가합니다.

따라서 $a_k=\lfloor s_k \rfloor$는 $\vec r$과의 외적의 부호를 바꾸지 않으면서 $\vec r_{k-2}$에 더할 수 있는 $\vec r_{k-1}$ 벡터의 최대 정수 개수입니다.

다시 말해, $a_k$는 $\vec r$이 정의하는 직선을 넘지 않으면서 $\vec r_{k-2}$에 $\vec r_{k-1}$을 더할 수 있는 최대 정수 횟수입니다:

$r=\frac{7}{9}=[0;1,3,2]$의 근사분수들입니다. 반수렴자는 회색 화살표 사이의 중간 점에 해당합니다.

위 그림에서 $\vec r_2 = (4;3)$은 $\vec r_0 = (1;0)$에 $\vec r_1 = (1;1)$을 반복해서 더하여 얻습니다.

$y=rx$ 직선을 넘지 않고는 더 이상 $\vec r_0$에 $\vec r_1$을 더할 수 없게 되면, 반대편으로 넘어가 $\vec r_1$에 $\vec r_2$를 반복해서 더하여 $\vec r_3 = (9;7)$을 얻습니다.

이 과정은 직선에 가까워지면서 지수적으로 길어지는 벡터들을 생성합니다.

이러한 성질 때문에 Boris Delaunay는 연속된 근사분수 벡터를 생성하는 이 과정을 코 늘이기 알고리즘이라고 불렀습니다.

점 $\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.$$

픽의 정리와 결합하면, 이는 삼각형의 내부(경계 제외)에는 격자점이 없으며 경계 위의 격자점은 $\vec 0$과 $0 \leq t \leq a_k$를 만족하는 모든 정수 $t$에 대한 $\vec r_{k-2} + t \cdot \vec r_{k-1}$뿐임을 의미합니다. 가능한 모든 $k$에 대해 이를 합치면, 짝수 번째 근사분수 벡터들과 홀수 번째 근사분수 벡터들이 이루는 다각형 사이의 공간에는 정수 점이 없음을 의미합니다.

이는 다시, 홀수 계수를 가진 $\vec r_k$들은 직선 $y=rx$ 위쪽에 있는 $x \geq 0$인 격자점들의 볼록 껍질을 이루고, 짝수 계수를 가진 $\vec r_k$들은 직선 $y=rx$ 아래쪽에 있는 $x > 0$인 격자점들의 볼록 껍질을 이룬다는 것을 의미합니다.

정의

이 다각형들은 클라인 다각형이라고도 하며, 연분수에 대한 이 기하학적 해석을 처음 제안한 Felix Klein의 이름을 딴 것입니다.

문제 예시

가장 중요한 사실과 개념을 소개했으니, 이제 구체적인 문제 예시를 살펴볼 차례입니다.

직선 아래의 볼록 껍질

$r=[a_0;a_1,\dots,a_k]=\frac{p_k}{q_k}$에 대해 $0 \leq x \leq N$이고 $0 \leq y \leq rx$인 격자점 $(x;y)$들의 볼록 껍질을 구하십시오.

풀이

무한 집합 $0 \leq x$를 고려한다면, 위쪽 볼록 껍질은 직선 $y=rx$ 자체로 주어질 것입니다.

하지만 $x \leq N$이라는 추가 제약이 있으면, 올바른 볼록 껍질을 유지하기 위해 결국 직선에서 벗어나야 합니다.

$t = \lfloor \frac{N}{q_k}\rfloor$라 하면, 껍질 위에서 $(0;0)$ 다음에 오는 처음 $t$개의 격자점은 정수 $1 \leq \alpha \leq t$에 대한 $\alpha \cdot (q_k; p_k)$입니다.

하지만 $(t+1)q_k$가 $N$보다 크므로 $(t+1)(q_k; p_k)$는 다음 격자점이 될 수 없습니다.

껍질의 다음 격자점으로 가려면, $x \leq N$을 유지하면서 $y=rx$에서 가장 작은 차이로 벗어나는 점 $(x;y)$로 가야 합니다.

$0 \leq x \leq 19$에서 $y=\frac{4}{7}x$ 아래에 있는 격자점들의 볼록 껍질은 점 $(0;0), (7;4), (14;8), (16;9), (18;10), (19;10)$으로 이루어집니다.

볼록 껍질에서 현재 마지막 점을 $(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$의 아래쪽 반수렴자여야 합니다.

따라서 $(\Delta x; \Delta y)$는 어떤 홀수 $i$와 $0 \leq t < a_i$에 대해 $(q_{i-1}; p_{i-1}) + t \cdot (q_i; p_i)$ 꼴입니다.

이러한 $i$를 찾기 위해, 가능한 모든 $i$를 가장 큰 것부터 차례로 살펴보면서 $N-x-q_{i-1} \geq 0$을 만족하는 $i$에 대해 $t = \lfloor \frac{N-x-q_{i-1}}{q_i} \rfloor$를 사용할 수 있습니다.

$(\Delta x; \Delta y) = (q_{i-1}; p_{i-1}) + t \cdot (q_i; p_i)$이면, 반수렴자의 성질에 의해 조건 $\Delta y \leq r \Delta x$가 유지됩니다.

또한 $i+2$에서 얻은 반수렴자를 이미 모두 사용했으므로 $x + q_{i-1} + a_i q_i = x+q_{i+1}$이 $N$보다 크고, 따라서 $t < a_i$가 성립합니다.

이제 $N$을 넘기 전까지 $(x;y)$에 $(\Delta x; \Delta y)$를 $k = \lfloor \frac{N-x}{\Delta x} \rfloor$번 더할 수 있으며, 그 후에는 다음 반수렴자를 시도합니다.

// 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);
}

# 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$이 주어집니다. $Ax + By \leq N$을 만족하면서 $Ax + By$가 최대가 되는 $x \geq 0$과 $y \geq 0$을 찾으십시오.

풀이

이 문제에서는 $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$의 방향을 뒤집겠습니다. 그러면 이제 $0 \leq x \leq \lfloor \frac{N}{A} \rfloor$, $By - Ax \leq N \;\bmod\; A$를 만족하면서 $By - Ax$가 최대가 되는 점 $(x; y)$를 찾아야 합니다. 각 $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$를 계산할 때 이 점을 염두에 두어야 합니다.

# (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 ex \rfloor$를 계산하십시오. 여기서 $e = [2; 1, 2, 1, 1, 4, 1, 1, 6, 1, \dots, 1, 2n, 1, \dots]$는 오일러 수이고 $N \leq 10^{4000}$입니다.

풀이

이 합은 $1 \leq x \leq N$과 $1 \leq y \leq ex$를 만족하는 격자점 $(x;y)$의 개수와 같습니다.

$y=ex$ 아래에 있는 점들의 볼록 껍질을 구성한 후, 이 개수는 픽의 정리를 이용해 계산할 수 있습니다:

// 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;
}

# 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$를 모두 더하는 것은 이전 문제에서 이미 할 수 있게 된 일입니다.

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

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]$ 구간에서 이 직선에 가장 가까운 격자점까지 이 볼록 껍질을 구성하는 방법도 이미 알고 있습니다(이는 위의 "Crime and Punishment" 문제에서 다루었습니다.

이제 직선에 가장 가까운 점에 도달하고 나면 직선이 실제로 그 가장 가까운 점을 지난다고 가정해도 된다는 점에 주목해야 합니다. 실제 직선과, 가장 가까운 점을 지나도록 약간 아래로 옮긴 직선 사이에는 $[0, N-1]$ 위의 다른 격자점이 없기 때문입니다.

따라서 $[0, N-1]$에서 직선 $y=\frac{Ax+B}{M}$ 아래의 볼록 껍질 전체를 구성하려면, $[0, N-1]$에서 직선에 가장 가까운 점까지 볼록 껍질을 구성한 다음 $B=0$일 때 볼록 껍질을 구성하는 알고리즘을 재사용하여 직선이 이 점을 지나는 것처럼 계속 진행하면 됩니다.

# 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

$1 \leq p, q \leq 10^9$를 만족하는 유리수 $\frac{p}{q}$가 있습니다. 여러 소수 $m$에 대해 $m \sim 10^9$을 법으로 하는 $p q^{-1}$의 값을 질의할 수 있습니다. $\frac{p}{q}$를 복원하십시오.

동등한 표현: $1 \leq x \leq N$에서 $Ax \;\bmod\; M$의 최솟값을 주는 $x$를 찾으십시오.

풀이

중국인의 나머지 정리에 의해, 여러 소수를 법으로 하는 결과를 질의하는 것은 그 소수들의 곱을 법으로 하는 결과를 질의하는 것과 같습니다. 이 때문에 일반성을 잃지 않고 충분히 큰 수 $m$을 법으로 하는 나머지를 알고 있다고 가정하겠습니다.

주어진 나머지 $r$에 대해 $p \equiv qr \pmod m$의 해 $(p, q)$는 여러 개 있을 수 있습니다. 하지만 $(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}$이면 $1 \leq p, q \leq 10^9$인 해 $\frac{p}{q}$는 유리수로서 유일합니다.

따라서 문제는 $m$을 법으로 하는 $r$이 주어졌을 때 $1 \leq q \leq 10^9$이고 $qr \;\bmod\; m \leq 10^9$를 만족하는 임의의 $q$를 찾는 것으로 귀결됩니다.

이는 사실상 $1 \leq q \leq 10^9$에서 가능한 가장 작은 $qr \bmod m$을 주는 $q$를 찾는 것과 같습니다.

$qr = km + b$일 때, 이는 $1 \leq q \leq 10^9$이면서 $qr - km \geq 0$이 가능한 한 최소가 되는 쌍 $(q, m)$을 찾아야 한다는 뜻입니다.

$m$은 상수이므로 $m$으로 나눌 수 있으며, 문제를 $1 \leq q \leq 10^9$이면서 $\frac{r}{m} q - k \geq 0$이 가능한 한 최소가 되는 $q$를 찾는 문제로 다시 쓸 수 있습니다.

연분수의 관점에서 이는 $\frac{k}{q}$가 $\frac{r}{m}$에 대한 최량 디오판토스 근사라는 뜻이며, $\frac{r}{m}$의 아래쪽 반수렴자만 확인하면 충분합니다.

# 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]

연습 문제