循环连分数

线性分式变换¶

连分数的另一个重要概念是所谓的线性分式变换（Linear fractional transformations）。

定义¶

线性分式变换 是一个函数 \(f : \mathbb R \to \mathbb R\)，使得 \(f(x) = \frac{ax+b}{cx+d}\) 对于一些 \(a,b,c,d \in \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))\) 也是线性分式变换：

线性分式变换的逆，也是线性分式变换：

DMOPC '19 Contest 7 P4 - Bob and Continued Fractions¶

给您一个正整数数组 \(a_1, \dots, a_n\)。您需要回答 \(m\) 查询。每个查询都要计算 \([a_l; a_{l+1}, \dots, a_r]\)。

例题 连分数的线性分式变换¶

设 \(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=[a_0; a_1, \dots, a_n]\)，\(B=[b_0; b_1, \dots, b_m]\)。计算 \(A+B\) 和 \(A \cdot B\) 的连分数表示。

循环连分数¶

仿照循环小数的概念，如果在连分数后面形成了循环，则形成 循环连分数。

循环连分数的最小正周期一般用字母 \(l\) 来表示，循环的部分称为循环节。容易发现，进入循环的充要条件是余项 \(r_k\) 的重复出现。

纯循环连分数¶

定义¶

如果存在 \(k\) 使得 \(x = [a_0; a_1, \dots, a_k, x]\)，则连分数 \(x = [a_0; a_1, \dots]\) 被称为 纯循环（periodic）。

如果 \(x = [a_0; a_1, \dots, a_k, y]\)，其中 \(y\) 是纯循环，则连分数 \(x = [a_0; a_1, \dots]\) 被称为 混循环（eventually periodic）。

例如纯循环连分数：

\([a_0,a_1,a_2,a_3,a_0,a_1,a_2,a_3,…]=[\overline{a_0,a_1,a_2,a_3}]\)

混循环连分数：

\([a_0,a_1,a_2,a_3,a_1,a_2,a_3,…]=[a_0,\overline{a_1,a_2,a_3}]\)

混循环连分数后面循环的部分，最早循环的部分称为它的「纯循环部分」。

纯循环连分数的整数部分 \(a_0\) 直接进入到了循环里面，因此要求 \(a_0\) 必须至少是 1。因此，纯循环连分数大于 1。

对于 \(x = [1; 1, 1, \dots]\)，有 \(x = 1 + \frac{1}{x}\)，因此 \(x^2 = x + 1\)。在循环连分数和二次方程之间存在着一般的联系。考虑以下等式：

一方面，这个方程意味着 \(x\) 的连分数表示是周期为 \(k+1\)。

另一方面，使用渐进分数的公式，这个方程意味着

也就是说，\(x\) 是其自身的线性分数变换。从等式中可以看出，\(x\) 是二次方程式的根：

类似的推理代表了混循环连分数，即 \(x = [a_0; a_1, \dots, a_k, y]\) 代表 \(y=[b_0; b_1, \dots, b_k, y]\)。事实上，从第一个方程导出 \(x = L_0(y)\)，从第二个方程导出 \(y = L_1(y)\)，其中 \(L_0\) 和 \(L_1\) 是线性分数变换。因此

可以进一步证明（由拉格朗日首先完成），对于具有整数系数的任意二次方程 \(ax^2+bx+c=0\)，其解 \(x\) 是一个混循环连分数。

拉格朗日连分数定理¶

关于循环连分数有一个特别重要的基础定理：

拉格朗日连分数定理：实二次有理数与循环连分数一一对应。

该定理最早由拉格朗日（Lagrange）于 1769 年证明。

根据余项的重复出现，证明循环连分数一定是二次有理数非常容易。重点在于证明二次有理数一定是循环连分数。

例题 二次有理数¶

找到 \(\alpha = \frac{x+y\sqrt{n}}{z}\) 的连分数，其中 \(x, y, z, n \in \mathbb Z\) 和 \(n > 0\) 不是完全平方。

# compute the continued fraction of sqrt(n)
def sqrt(n):
    n0 = math.floor(math.sqrt(n))
    x, y, z = 0, 1, 1
    a = []
    def step(x, y, z):
        a.append((x * n0 + y) // z)
        t = y - a[-1]*z
        x, y, z = z*t, -z*y, 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

伽罗瓦连分数定理¶

纯循环连分数有类似于反序定理的定理。记：

则两个 x 互为「倒数负共轭」。即，若记：

则 x 与 y 共轭。

该定理最早由伽罗瓦（Galois）在他 17 岁那年（1828 年）发现并证明。

根据伽罗瓦连分数定理，纯循环连分数的共轭一定落在 \(-1\) 到 \(0\) 之间。考虑它的逆问题，就有下面的定理：

如果二次有理数 \(x\) 大于 \(1\)，它的共轭落在 \(-1\) 到 \(0\) 之间，则它是纯循环连分数。

根号 d 的连分数¶

对于非平方整数 d，有：

这是根号 d 的连分数形式。不仅最后一位是整数部分的两倍，而且中间部分还是对称的。

同样也可以证明，整数部分为半整数的相同结构：

一个实例¶

以 \(\sqrt{74}\) 为例，实际运算一下连分数算法：

各个余项分别是：

根据伽罗瓦连分数定理，对称的余项 \(r_k\) 和 \(r_{l+1-k}\) 循环部分恰好相反，因此互为倒数负共轭。

如果循环节 \(l\) 是奇数，则对称部分最中间的余项与自己互为倒数负共轭。但如果循环节 \(l\) 是偶数，就不会出现这种现象。

在后面的 Pell 方程一节中将指出，在根号 \(d\) 的连分数中，循环节 \(l\) 的奇偶性，将直接决定 Pell 方程中的 \(-1\) 形式是否有解。

Tavrida NU Akai Contest - Continued Fraction¶

你得到了 \(x\) 和 \(k\),\(x\) 不是一个完全平方数。让 \(\sqrt x = [a_0; a_1, \dots]\)，找到 \(\frac{p_k}{q_k}=[a_0; a_1, \dots, a_k]\) 的 \(0 \leq k \leq 10^9\)。

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