Project Euler 846

题目

Magic Bracelets

A bracelet is made by connecting at least three numbered beads in a circle. Each bead can only display \(1\), \(2\), or any number of the form \(p^k\) or \(2p^k\) for odd prime \(p\).

In addition a magic bracelet must satisfy the following two conditions:

• no two beads display the same number
• the product of the numbers of any two adjacent beads is of the form \(x^2+1\)

Define the potency of a magic bracelet to be the sum of numbers on its beads.

The example is a magic bracelet with five beads which has a potency of 155.

Let \(F(N)\) be the sum of the potency of each magic bracelet which can be formed using positive integers not exceeding \(N\), where rotations and reflections of an arrangement are considered equivalent. You are given \(F(20)=258\) and \(F(10^2)=538768\).

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

解决方案

若奇素数 \(p\mid (x^2+1)\)，则\(x^2\equiv -1\pmod p,\)即 \(-1\) 是模 \(p\) 的二次剩余。经典结论：对奇素数 \(p\)，\(-1\) 为二次剩余当且仅当 \(p\equiv 1\pmod 4\)。因此：若 \(p\equiv 3\pmod 4\)，则 \(p\nmid (x^2+1)\)，不可能整除任何相邻乘积；所以任何能出现在手链上的数，其奇素因子只能来自 \(p\equiv 1\pmod 4\) 的素数；\(2\) 当然可能出现。

因此，实际相关的数字集合就是\(\mathcal{V}(N)=\{1,2\}\cup\{p^k\le N\}\cup\{2p^k\le N\},\)其中 \(p\) 遍历所有满足 \(p\equiv 1\pmod 4\) 的奇素数。

在高斯整数环 \(\mathbb{Z}[i]\) 中，定义范数\(N(a+bi)=a^2+b^2.\)它满足乘法性：\(N(zw)=N(z)N(w)\)，且 \(\mathbb{Z}[i]\) 是唯一分解整环。

对 \(p\equiv 1\pmod 4\) 的奇素数，\(p\) 在 \(\mathbb{Z}[i]\) 中可分解为\(p=\pi\overline{\pi},N(\pi)=p,\)从而 \(p=a^2+b^2\) 有一组互素整数解 \(\gcd(a,b)=1\)。并且对\(n\in\{1,2,p^k,2p^k\},\)存在（并且在互素、忽略符号与交换 \(a,b\)意义下唯一）的表示\(n=a^2+b^2, \gcd(a,b)=1.\)

其依据是：\(\mathbb{Z}[i]\) 为唯一分解整环，因此 \(n\) 在 \(\mathbb{Z}[i]\) 中的分解（除去单位元与共轭）唯一；而 \(\gcd(a,b)=1\) 的平方和表示 \(n=a^2+b^2\)，恰与把 \(n\) 表示为某个高斯整数范数 \(N(a+bi)\) 一一对应（同样按单位元与共轭同类视为等价）。

于是我们可以把每个可用数字 \(n\) 唯一对应到一个原始向量 \((a,b)\)（约定 \(a\ge b\ge 0\)），并把它看作一个分数\(\dfrac{b}{a}\in[0,1].\)

接下来我们把相邻条件 \(uv=x^2+1\) 转化为：对应既约分数在 Farey 图中相邻（Farey 图是一张以所有既约有理数（并加入 \(\infty=1/0\)）为顶点的无穷图。两点 \(\dfrac{a}{b},\dfrac{c}{d}\) 之间连边当且仅当 \(|ad-bc|=1\)，这正是 Farey 邻接判定）。

令\(u=a^2+b^2=N(a+bi), v=c^2+d^2=N(c+di),\)且 \(\gcd(a,b)=\gcd(c,d)=1\)。注意\(x^2+1 = N(x+i).\)由于范数乘法性，若 \(uv=x^2+1\)，则在 \(\mathbb{Z}[i]\) 中（忽略单位元）必有\((a+bi)(c-di)=\pm(x+i).\)展开得\((a+bi)(c-di)=(ac+bd) + (bc-ad)i.\) 因此虚部必须是 \(\pm 1\)，也就是\(|ad-bc|=1.\)这正是 Farey 图中的相邻判定：两既约分数 \(\dfrac{b}{a},\dfrac{d}{c}\) 相邻当且仅当\(|ad-bc|=1.\)

结论： 把每个可用数字 \(n\) 映射为分数 \(\dfrac{b}{a}\)，则数字相邻当且仅当Farey图上相邻。

于是魔法手链等价于：由这些分数组成的一个简单环（顶点互不相同），且环上相邻顶点在 Farey 图里相邻。

Farey 图有一个关键结构：任意既约分数 \(\dfrac{b}{a}\in(0,1)\) 都有且仅有两个更小的相邻点（可视作父节点），记为\(\dfrac{b_1}{a_1}<\dfrac{b}{a}<\dfrac{b_2}{a_2},\)满足三者两两相关的行列式条件，并且\(\dfrac{b}{a}=\dfrac{b_1+b_2}{a_1+a_2},\)也就是 \(\dfrac{b}{a}\) 是其两父节点的mediant。

把这用向量描述就是：给定互素 \((a,b)\)，存在唯一的拆分\((a,b)=(c,d)+(a-c,b-d),\)并且选定一个方向\(ad-bc=1\)后，两父向量的范数更小。

在本题中，我们只保留那些父范数也在 \(\mathcal{V}(N)\) 的边；这样每个点的父节点数最多为 \(2\)。

考虑 Farey 图（或其诱导子图）中的任意一个简单环。设其中范数最大的顶点为 \(v\)，则 \(v\) 在该环上的两个邻点恰为它的两个父节点，且这两个父节点彼此相邻，因此可用它们之间的边替换经过 \(v\) 的两条边，从而将环长度减少 \(1\)。重复此消去过程，环会逐步缩短，直至退化为一条边。

另一方面，任何环都可以从一条 Farey 边 \((U,V)\) 出发，通过在某些边之间插入媒介点（再递归）逐步长出来。

因此，对每一条 Farey 边 \((U,V)\)（在诱导子图中存在），所有可形成的环与从 \(U\) 到 \(V\) 的路径一一对应：路径允许在区间内递归插入媒介点；取一条非平凡路径（至少包含一个内部点），再加上直接边 \((U,V)\)，就形成一个长度 \(\ge 3\) 的简单环；这样的底边在该环的归约过程中是唯一的，所以不会重复计数。

于是问题变为：对每条边 \((U,V)\)，统计从 \(U\) 到 \(V\) 的所有路径的数量与路径点权和（点权即数字本身），再把非平凡路径对应的环的 potency 累计。

把每条有向边记为子 \(\to\) 父。对一条边 \((i\to q)\)，维护两个量：

• \(C(i,q)\)表示从 \(q\) 到 \(i\) 的非平凡路径条数（内部至少一个点）；
• \(S(i,q)\)表示这些非平凡路径的内部点权和总和（不含端点 \(q,i\)）。

初始时：\(C(i,q)=0, S(i,q)=0\)对所有边都成立；

当我们要把这些路径闭合成环时：每条非平凡路径加上底边 \((q,i)\) 就得到一个 bracelet，其 potency 为\(q+i+\sum_{w\in I} w\)，其中 \(I\) 表示该路径的内部顶点集合。

因此，边 \((i,q)\) 对答案的累积贡献为\(C(i,q)(q+i)+S(i,q).\)

随后，把直接边路径（无内部点）并入状态：\(C(i,q)\leftarrow C(i,q)+1,\)对应地 \(S(i,q)\) 不变（因为该路径内部点和为 \(0\)）。

若点 \(i\) 有两个父节点 \(r<s\)，构成三角形 \((r,i,s)\)。则所有经过 \(i\) 的 \(r\to s\) 路径由一条 \(r\to i\) 路径与一条 \(i\to s\) 路径拼接而成，所以直接做如下累积更新：

\[ C(s,r)\leftarrow C(i,r)\cdot C(i,s),S(s,r)\leftarrow C(i,r)C(i,s)\cdot i + C(i,r)\cdot S(i,s)+C(i,s)\cdot S(i,r). \]

这里 \(C(i,r),C(i,s)\) 已经包含了平凡路径（直接边）。按数值从大到小处理 \(i\)，即可保证更新时所需子状态都已就绪。

最终，我们枚举所有互素对 \((a,b)\)（\(1\le b\le a,a^2+b^2\le N\)），若 \(a^2+b^2\) 是素数且 \(\equiv 1\pmod 4\)，记录该素数\(p\)对应的 \((a,b)\)。随后对每个这样的素数 \(p\)，用高斯整数乘法生成 \((a+bi)^k\) 的表示，从而得到 \(p^k\) 与 \(2p^k\) 的平方和表示，标记为可用点。那么对每个可用点 \(n=a^2+b^2\)，通过解\(a d - b c = 1\)得到两父向量 \((c,d)\) 与 \((a-c,b-d)\)，进而得到父节点范数（并检查是否也在可用集合中）。最终用DP得到答案。

代码

from array import array
from tools import Factorizer, int_sqrt, gcd

N = 10 ** 6

def F(N):
    factorizer = Factorizer(N)
    primes = factorizer.primes
    is_prime = [i >= 2 and factorizer.spf[i] == i for i in range(N+1)]
    lim = int_sqrt(N)

    base_a = array('H', [0]) * (N + 1)
    base_b = array('H', [0]) * (N + 1)

    for b in range(1, lim + 1):
        b2 = b * b
        a_max = int_sqrt(N - b2)
        for a in range(b, a_max + 1):
            if gcd(a, b) != 1:
                continue
            s = a * a + b2
            if s <= N and is_prime[s] and (s & 3) == 1 and base_a[s] == 0:
                base_a[s] = a
                base_b[s] = b

    rep_x = array('H', [0]) * (N + 1)
    rep_y = array('H', [0]) * (N + 1)
    has = bytearray(N + 1)

    def normalize(x, y):
        return sorted([abs(x), abs(y)])

    def gmul(x1, y1, x2, y2):
        return x1 * x2 - y1 * y2, x1 * y2 + y1 * x2

    has[1] = 1
    rep_x[1] = 1
    rep_y[1] = 0
    if N >= 2:
        has[2] = 1
        rep_x[2] = 1
        rep_y[2] = 1

    for p in primes:
        if p == 2:
            continue
        if base_a[p] == 0:
            continue
        a = int(base_a[p])
        b = int(base_b[p])
        curx, cury = a, b
        q = p
        while q <= N:
            nx, ny = normalize(curx, cury)
            has[q] = 1
            rep_x[q] = nx
            rep_y[q] = ny

            qq = q + q
            if qq <= N:
                tx, ty = gmul(curx, cury, 1, 1)
                tx, ty = normalize(tx, ty)
                has[qq] = 1
                rep_x[qq] = tx
                rep_y[qq] = ty

            if q > N // p:
                break
            curx, cury = gmul(curx, cury, a, b)
            q *= p

    par_len = bytearray(N + 1)
    par0 = array('I', [0]) * (N + 1)
    par1 = array('I', [0]) * (N + 1)

    if N >= 2:
        par_len[2] = 1
        par0[2] = 1

    for n in range(3, N + 1):
        if not has[n]:
            continue
        x = int(rep_x[n])
        y = int(rep_y[n])

        if y == 1:
            d = 1
            c = x - 1
        else:
            d = pow(x, -1, y)
            c = (x * d - 1) // y

        n1 = c * c + d * d
        c2 = x - c
        d2 = y - d
        n2 = c2 * c2 + d2 * d2

        if n1 > n2:
            n1, n2 = n2, n1

        cnt = 0
        if n1 <= N and has[n1]:
            par0[n] = n1
            cnt = 1
        if n2 <= N and has[n2] and n2 != n1:
            if cnt == 0:
                par0[n] = n2
                cnt = 1
            else:
                par1[n] = n2
                cnt = 2
        par_len[n] = cnt

    dp_c0 = array('Q', [0]) * (N + 1)
    dp_s0 = array('Q', [0]) * (N + 1)
    dp_c1 = array('Q', [0]) * (N + 1)
    dp_s1 = array('Q', [0]) * (N + 1)

    ans = 0

    for i in range(N, 0, -1):
        l = par_len[i]
        if l == 0:
            continue

        p = par0[i]
        c = dp_c0[i]
        s = dp_s0[i]
        ans += c * (i + p) + s
        dp_c0[i] = c + 1

        if l == 2:
            p = par1[i]
            c = dp_c1[i]
            s = dp_s1[i]
            ans += c * (i + p) + s
            dp_c1[i] = c + 1

            u = par0[i]
            v = par1[i]

            w0 = dp_c0[i]
            w1 = dp_c1[i]
            s0 = dp_s0[i]
            s1 = dp_s1[i]

            if par0[v] == u:
                dp_c0[v] += w0 * w1
                dp_s0[v] += w0 * w1 * i + w0 * s1 + w1 * s0
            else:
                dp_c1[v] += w0 * w1
                dp_s1[v] += w0 * w1 * i + w0 * s1 + w1 * s0

    return ans

print(F(N))