Project Euler 752

发表于 分类于 阅读次数： 本文字数： 4.8k 阅读时长 ≈ 4 分钟

Project Euler 752

题目

Powers of \(1+\sqrt 7\)

When \((1+\sqrt 7)\) is raised to an integral power, \(n\), we always get a number of the form \((a+b\sqrt 7)\).

We write \((1+\sqrt 7)^n = \alpha(n) + \beta(n)\sqrt 7\).

For a given number \(x\) we define \(g(x)\) to be the smallest positive integer \(n\) such that:

\[\begin{align} \alpha(n) &\equiv 1 \pmod x\qquad \text{and }\\ \beta(n) &\equiv 0 \pmod x\end{align} \]

and \(g(x) = 0\) if there is no such value of \(n\). For example, \(g(3) = 0\), \(g(5) = 12\).

Further define

\[ G(N) = \sum_{x=2}^{N} g(x)\]

You are given \(G(10^2) = 28891\) and \(G(10^3) = 13131583\).

Find \(G(10^6)\).

解决方案

本题的核心是把\((1+\sqrt7)^n=\alpha(n)+\beta(n)\sqrt7\)的同余条件转化为一个固定 \(2\times2\) 矩阵在模 \(x\) 意义下的乘法阶问题，随后利用中国剩余定理（CRT）与素幂提升把所有 \(g(x)\) 快速算出，从而在 \(N=10^6\) 范围内线性求和。

令\((1+\sqrt7)^n=\alpha(n)+\beta(n)\sqrt7.\)两边同乘 \((1+\sqrt7)\) 得\((1+\sqrt7)^{n+1}=(\alpha(n)+\beta(n)\sqrt7)(1+\sqrt7)= \bigl(\alpha(n)+7\beta(n)\bigr)+\bigl(\alpha(n)+\beta(n)\bigr)\sqrt7.\)

因此\(\alpha(n+1)=\alpha(n)+7\beta(n),\beta(n+1)=\alpha(n)+\beta(n).\)

把它写成矩阵形式，令列向量\(v(n)=\begin{pmatrix}\alpha(n)\\\beta(n)\end{pmatrix},\)转移矩阵\(T=\begin{pmatrix}1&7\\1&1\end{pmatrix},\)则

\[v(n+1)=Tv(n).\]

又因为 \((1+\sqrt7)^0=1\)，所以\(v(0)=\begin{pmatrix}1\\ 0\end{pmatrix}.\)于是\(v(n)=T^nv(0).\)

更进一步，注意到在二次整数环 \(\mathbb Z[\sqrt7]\) 中，乘以 \(\alpha+\beta\sqrt7\)在基 \((1,\sqrt7)\) 下对应线性变换：\((a+b\sqrt7)(\alpha+\beta\sqrt7)=(a\alpha+7b\beta)+(a\beta+b\alpha)\sqrt7,\)因此对应矩阵为\(M(\alpha,\beta)=\begin{pmatrix}\alpha&7\beta\\ \beta&\alpha\end{pmatrix}.\)

当 \((\alpha,\beta)=(1,1)\) 时正是 \(T\)，而 \((1+\sqrt7)^n=(\alpha(n),\beta(n))\) 的乘法矩阵就是 \(T^n\)，从而得到一个关键恒等式：

\[ T^n=\begin{pmatrix}\alpha(n)&7\beta(n)\\ \beta(n)&\alpha(n)\end{pmatrix}. \]

题目定义 \(g(x)\) 为最小正整数 \(n\)，使得\(\alpha(n)\equiv 1\pmod x, \beta(n)\equiv 0\pmod x.\)结合上面的矩阵表达式，有 \(\alpha(n)\equiv1, \beta(n)\equiv0\pmod x \iff T^n\equiv I_2\pmod x.\)因此：\(g(x)=\operatorname{ord}_x(T),\)即 \(T\) 在模 \(x\) 的可逆矩阵群中对应的乘法阶。

可见，\(\det(T)=1\cdot1-7\cdot1=-6.\)如果 \(T^n\equiv I_2\pmod x\) 对某个 \(n\) 成立，那么 \(T\) 在模 \(x\) 下必然可逆，因而必须有 \(\gcd(\det(T),x)=1 \iff \gcd(6,x)=1.\)于是得到必要条件：\(\gcd(x,6)\neq1 \Rightarrow g(x)=0.\)

反过来，若 \(\gcd(x,6)=1\)，则 \(T\in GL(2,\mathbb Z/x\mathbb Z)\)，该群有限，因此 \(T\) 一定有有限阶，从而 \(g(x)>0\)。所以： \(g(x)>0 \iff \gcd(x,6)=1.\)

若 \(\gcd(a,b)=1\)，则 CRT 给出环同构\(\mathbb Z/(ab)\mathbb Z \cong \mathbb Z/a\mathbb Z \times \mathbb Z/b\mathbb Z,\) 从而\(GL(2,\mathbb Z/(ab)\mathbb Z)\cong GL(2,\mathbb Z/a\mathbb Z)\times GL(2,\mathbb Z/b\mathbb Z).\)

在直积群中元素的阶满足\(\operatorname{ord}(u,v)=\mathrm{lcm}(\operatorname{ord}(u),\operatorname{ord}(v)).\) 因此\(\gcd(a,b)=1\Rightarrow g(ab)=\mathrm{lcm}(g(a),g(b)).\)

也就是说，这一步是把计算所有\(x\)变成只要知道素幂即可。

设 \(p\) 为奇素数且 \(p\nmid 6\)，考虑 \(T\) 在域 \(\mathbb F_p\) 上的作用。对任意非零向量 \(v\in\mathbb F_p^2\setminus\{0\}\)，因为 \(T\) 可逆，其轨道\(v,Tv,T^2v,\dots\)都落在有限集合中，必循环。

尤其对 \(v(0)=(1,0)^T\)，其最小循环长度正是满足\(T^n v(0)=v(0)\)的最小正 \(n\)，而由前面矩阵形式这等价于 \(T^n=I_2\)，所以该循环长度就是 \(g(p)\)。

同时，轨道大小一定整除集合大小 \(|\mathbb F_p^2\setminus\{0\}|=p^2-1\)，因此得到：\(g(p)\mid p^2-1 (p\ge5,p\ne7).\)

注意，特殊素数 \(p=7\) 会出现异常。这是因为当 \(p=7\) 时矩阵 \(T\) 在模 \(7\) 下退化为\(T\equiv\begin{pmatrix}1&0\\1&1\end{pmatrix}=I+N, N^2=0,\)从而\(^n\equiv I+nN\pmod7,\)因此 \(T^n\equiv I\) 的最小正解为 \(n=7\)，即\(g(7)=7.\)这也是 \(p=7\) 不满足 \(g(p)\mid(p^2-1)\) 的根本原因。

设 \(p\ge5\) 且 \(p\nmid 6\)，令 \(t=g(p)\)，则\(T^t\equiv I_2\pmod p.\)记\(U=T^t.\)由于 \(U\equiv I_2\pmod p\)，存在整数矩阵 \(A\) 使得\(U=I_2+pA.\)

为了严格描述提升到模 \(p^k\) 时阶如何变化，对整数矩阵 \(M\) 定义\(\displaystyle{v_p(M)=\min_{i,j} \{v_p(m_{ij})\},}\) 即 \(M\) 所有元素的 \(p\)-进指数的最小值。于是存在唯一整数 \(s\ge1\) 满足\(U\equiv I_2\pmod{p^s}\land U\not\equiv I_2\pmod{p^{s+1}},\)也就是\(v_p(U-I_2)=s.\)

接下来给出一个标准的提升引理。设 \(p\) 为奇素数，且\(V=I_2+p^sB, v_p(B)=0.\)则\(V^p\equiv I_2\pmod{p^{s+1}}, V^p\not\equiv I_2\pmod{p^{s+2}}.\)

证明如下：由二项式展开\(\displaystyle{(I_2+p^sB)^p=\sum_{i=0}^{p}\binom{p}{i}(p^sB)^i.}\)其中 \(i=1\) 项为 \(p^{s+1}B\)；而对 \(2\le i\le p\)，因为 \(p\mid\dbinom{p}{i}\)，第 \(i\) 项至少含有因子\(p\cdot p^{si}=p^{si+1}\ge p^{2s+1}\ge p^{s+2}(s\ge1),\)因此\(V^p\equiv I_2+p^{s+1}B\pmod{p^{s+2}}.\)

又因 \(v_p(B)=0\)，故 \(p^{s+1}B\not\equiv0\pmod{p^{s+2}}\)，从而结论成立。

由引理可归纳推出更一般的结论：对任意 \(j\ge0\)，都有\(V^{p^j}\equiv I_2\pmod{p^{s+j}}, V^{p^j}\not\equiv I_2\pmod{p^{s+j+1}}.\)

将其应用到 \(U\)：把 \(U\) 写成\(U=I_2+p^sB, v_p(B)=0.\)我们需要最小的 \(m\) 使得\(U^m\equiv I_2\pmod{p^k}.\)那么

  • 若 \(k\le s\)，则 \(U\equiv I_2\pmod{p^k}\) 已经成立，因此最小 \(m=1\)。
  • 若 \(k>s\)，由上面的归纳结论知最小 \(m=p^{k-s}\)。

于是得到统一公式\(g(p^k)=t\cdot p^{\max(0,k-s)}\).

特别地，在最常见的情形\(s=v_p(T^t-I_2)=1\Longleftrightarrow T^t\equiv I_2\pmod p\land T^t\not\equiv I_2\pmod{p^2},\)上式化简为

\[g(p^k)=t\cdot p^{k-1}.\]

因此，每提升一层模 \(p\) 的精度，阶会乘一个 \(p\)这一现象；而若出现 \(s\ge2\)，则意味着 \(T^t\) 在模 \(p^2\) 下已经更接近单位矩阵，使得阶的增长会相应变慢，但仍由上面的公式完全刻画。

最终，把 \(x\) 分解得到：\(\displaystyle{x=\prod_{i=1}^{t} p_i^{e_i}, \gcd(x,6)=1.}\)由素幂公式与互素 lcm 公式，得到

\[g(x)=\mathrm{lcm}\Bigl(g(p_1)p_1^{e_1-1},g(p_2)p_2^{e_2-1},\dots\Bigr).\]

若 \(\gcd(x,6)\ne1\)，则 \(g(x)=0\)。

最终我们需要回到如何求\(g(p)\)的问题上，已知 \(g(p)\mid p^2-1\)，令\(M=p^2-1=(p-1)(p+1).\)先分解 \(M\) 的质因子集合 \(\{q\}\)，再用标准降阶算法：

从 \(O=M\) 开始，对每个质因子 \(q\)：只要\(O\equiv0\pmod q\land T^{O/q}\equiv I_2\pmod p,\) 就更新 \(O\leftarrow O/q\)。最终得到的 \(O\) 就是最小的 \(g(p)\)。

因为 \(p\le10^6\)，而 \(M=(p-1)(p+1)\) 的分解只需分解两个不超过 \(10^6+1\) 的数，用最小质因子筛即可完成得非常快。

代码

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
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
from math import gcd
from tools import *

N = 10 ** 6

def mat_mul(A, B, mod):
    a, b, c, d = A
    e, f, g, h = B
    return ((a * e + b * g) % mod,
            (a * f + b * h) % mod,
            (c * e + d * g) % mod,
            (c * f + d * h) % mod)

def mat_pow(A, n, mod):
    R = (1, 0, 0, 1)
    while n:
        if n & 1:
            R = mat_mul(R, A, mod)
        A = mat_mul(A, A, mod)
        n >>= 1
    return R

def compute_g_prime(N: int, factorizer):
    T = (1, 7, 1, 1)
    I = (1, 0, 0, 1)
    gp = [0] * (N + 1)
    for p in factorizer.primes:
        if p > N:
            break
        if p == 7:
            gp[p] = 7
            continue
        if gcd(p, 6) != 1:
            gp[p] = 0
            continue
        m = p * p - 1
        facs = [v[0] for v in factorizer.factor_small(p - 1)] + [v[0] for v in factorizer.factor_small(p + 1)]
        ordv = m
        for q in facs:
            while ordv % q == 0:
                if mat_pow(T, ordv // q, p) == I:
                    ordv //= q
                else:
                    break
        gp[p] = ordv
    return gp

def G(N: int):
    factorizer = Factorizer(N)
    gp = compute_g_prime(N, factorizer)
    ans = 0
    for x in range(2, N + 1):
        if gcd(x, 6) != 1:
            continue
        cur = 1
        for p, e in factorizer.factor_small(x):
            cur = lcm(cur, gp[p] * p ** (e - 1))
        ans += cur
    return ans

print(G(N))