Project Euler 312

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

Project Euler 312

题目

Cyclic paths on Sierpiński graphs

  • A Sierpiński graph of order-\(1\) \((S_1)\) is an equilateral triangle.

  • \(S_{n+1}\) is obtained from \(S_n\) by positioning three copies of \(S_n\) so that every pair of copies has one common corner.

Let \(C(n)\) be the number of cycles that pass exactly once through all the vertices of \(S_n\).

For example, \(C(3) = 8\) because eight such cycles can be drawn on \(S_3\), as shown below:

It can also be verified that :

\(\begin{aligned} &C(1) = C(2) = 1\\ &C(5) = 71328803586048\\ &C(10 000) \bmod 10^8 = 37652224\\ &C(10 000) \bmod 13^8 = 617720485 \end{aligned}\)

Find \(C(C(C(10 000))) \bmod 13^8\).

解决方案

令\(N=10000\)。把第 \(n\) 级的 Sierpiński gasket graph 记作 \(S_n\)。关键自相似性是：\(S_{n+1}\) 由三个拷贝 \(S_n\)（上、左下、右下）在角点处粘接得到；

三个外角点在整体结构里起到“边界端点”的作用，使得哈密顿回路/路径在粘接处的走法受到强约束。这让计数可以写成乘法型的递推。

令：

  • \(C(n)\)：\(S_n\) 上无向哈密顿回路的条数；
  • \(f(n)\)：\(S_n\) 上的“哈密顿路径开口态”条数：一条路径经过 \(S_n\) 的所有顶点恰好一次，并且路径两端点是三个外角点中的某特定的两点（把三种端点对称情况一起计数，因为对称性）。这也意味着，三个外角点只有一个在哈密顿路径内部。

直觉上：每条回路可以在外角的三个位置之一“剪开”成一条开口态路径，所以会出现一个乘 \(3\) 的关系。

好的，我把上条回答里的公式全部用 $...$ / $$...$$ 包起来（其余文字基本不变）。

当 \(n\) 足够大进入稳定递推区间时，把 \(S_n\) 看成由三个 \(S_{n-1}\) 组成。任何一条 \(S_n\) 上的哈密顿回路，在每个 \(S_{n-1}\) 子块内部都必须表现为“从该子块的一个粘接点进来、从另一个粘接点出去”的哈密顿路径开口态。

三块彼此独立选择，于是得到：\(C(n)=f^3(n-1).\)

把“端点是哪一对外角”也算进去时，每条回路对应 \(3\) 种对称的剪开方式，因此\(f(n)=3C(n).\)

把上述两条式子合并，消去 \(f\)，得到：

\[ C(n)=\bigl(C(n-1)\bigr)^3,n\ge 4, \]

两边取对数。令 \(x_n=\ln C(n)\)，则\(x_n = 3x_{n-1}+3\ln 3\)。

解这个线性递推：\(\displaystyle{ x_n =3^{n-3}\ln 8+\left(\sum_{k=0}^{n-4}3^k\right)\cdot 3\ln 3 =3^{n-3}\ln 8+\frac{3^{n-3}-1}{3-1}\cdot 3\ln 3 =3^{n-3}\ln 8+\frac{3^{n-2}-3}{2}\ln 3. }\)

重新指数化得到：\(C(n)=8^{3^{n-3}}\cdot 3^{(3^{n-2}-3)/2}.\)

化简后，即得到最简形式：

\[ C(n)=8\cdot 12^{\frac{3^{,n-2}-3}{2}},n\ge 3. \]

如果要算\(C(n)=8\cdot 12^{E},E=\dfrac{3^{n-2}-3}{2}.\)设模数是 \(13^k\)，则\(\varphi(13^k)=13^k-13^{k-1}=12\cdot 13^{k-1}.\)

且 \(\gcd(12,13^k)=1\)，因此由欧拉定理，得：\(12^{x+\varphi(13^k)}\equiv 12^x\pmod{13^k}.\)

所以只需要知道\(E \bmod \varphi(13^k).\)

关键难点：\(E\) 里有除以 \(2\)。做法是先在模 \(2\varphi(13^k)\) 下算出 \(3^{n-2}\)，再减 \(3\)，此时必为偶数，就能安全地“除以 2”：

令\(u\equiv 3^{n-2}\pmod{2\varphi(13^k)}.\)

则\(t\equiv (u-3)\pmod{2\varphi(13^k)}\)是偶数（因为 \(u\) 是奇数，\(3\) 也是奇数），于是\(E\equiv \dfrac{t}{2}\pmod{\Phi}.\)

最终：

\[ C(n)\equiv 8\cdot 12^{E\bmod \varphi(13^k)}\pmod{13^k}. \]

现在题目要求：\(C(C(C(N)))\bmod 13^8.\)

如果直接算内层 \(C(10000)\) 是天文数。因此我们用先取模再当作下一层的输入的方式把它压下来：

设\(m_0=13^k,m_1=\varphi(m_0)=12\cdot 13^{k-1},m_2=\varphi(m_1)=48\cdot 13^{k-2}.\)

然后计算

  • \(x_1=C(N)\bmod m_2,\)
  • \(x_2=C(x_1)\bmod m_1,\)
  • \(\text{ans}=C(x_2)\bmod m_0.\)

在算 \(C(\cdot)\bmod 13^8\)时，我们只会遇到形如 \(3^{n-2}\bmod(2\varphi(13^8))\) 这样的东西；而要算这个又会用到 \(n\) 在更小周期里的信息，于是出现再往下的 \(\varphi(\varphi(\cdot))\)。

可以注意到，当 \(n\ge 4\) 时，\(E\ge 3\)，所以 \(12^E\) 含有至少一个因子 \(12\)，从而\(C(n)=8\cdot 12^E\equiv 0\pmod{12},\)

甚至也有\(C(n)\equiv 0\pmod{48}.\)

因此像 \(m_1=12\cdot 13^7\)、\(m_2=48\cdot 13^6\) 这种模数，我们只要先算出模 \(13^7\) 或者是 \(13^6\) 的部分，再用 CRT 拼回去即可。

代码

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
from tools import CRT


def C_mod_13power(n, k):
    """
    Compute C(n) mod 13^k using:
        C(n) = 8 * 12^((3^(n-2) - 3)/2)  for n >= 3
        C(1)=C(2)=1
    """
    mod = 13 ** k
    if n <= 2:
        return 1 % mod
    if n == 3:
        return 8 % mod

    phi = 12 * (13 ** (k - 1))       # φ(13^k)
    u = pow(3, n - 2, 2 * phi)        # 3^(n-2) mod (2*phi)
    t = (u - 3) % (2 * phi)           # even
    e_mod = (t // 2) % phi            # ((3^(n-2)-3)/2) mod phi

    return (8 * pow(12, e_mod, mod)) % mod


def C_mod_s_times_13power(n, s, k):
    """
    Compute C(n) mod (s * 13^k), with gcd(s,13)=1 (here s is 12 or 48).
    Use:
      - r13 = C(n) mod 13^k
      - rs  = C(n) mod s
      - CRT([s, 13^k], [rs, r13])
    """
    mod13 = 13 ** k
    r13 = C_mod_13power(n, k)

    if n <= 2:
        rs = 1 % s
    elif n == 3:
        rs = 8 % s
    else:
        # n >= 4 => C(n) divisible by 12 and by 48
        rs = 0

    return CRT([s, mod13], [rs, r13])


if __name__ == "__main__":
    x1 = C_mod_s_times_13power(10000, 48, 6)  # mod m2
    x2 = C_mod_s_times_13power(x1,    12, 7)  # mod m1
    ans = C_mod_13power(x2, 8)                # mod m0
    print(ans)  # expected: 324681947