Project Euler 354

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

Project Euler 354

题目

Distances in a bee’s honeycomb

Consider a honey bee’s honeycomb where each cell is a perfect regular hexagon with side length \(1\).

One particular cell is occupied by the queen bee.

For a positive real number \(L\), let \(\text{B}(L)\) count the cells with distance \(L\) from the queen bee cell (all distances are measured from centre to centre); you may assume that the honeycomb is large enough to accommodate for any distance we wish to consider.

For example, \(\text{B}(\sqrt 3)=6, \text{B}(\sqrt {21}) = 12\) and \(\text{B}(111\,111\,111) = 54\).

Find the number of \(L \le 5 \times 10^{11}\) such that \(\text{B}(L) = 450\).

解决方案

令\(N= 5 \times 10^{11}\)蜂巢每个格子是边长为 \(1\) 的正六边形。相邻六边形的中心距为\(d=\sqrt{3}.\)把蜂巢中心点都放到平面上，它们构成一个三角晶格（等价于 Eisenstein 整数晶格的缩放）。

取复平面，令\(\omega=e^{i\pi/3}=\dfrac12+i\dfrac{\sqrt3}{2}.\)三角晶格可表示为\(\Lambda=\{a+b\omega: a,b\in\mathbb Z\}.\)蜂巢中心点与 \(\Lambda\) 只差一个整体缩放：由于相邻中心距是 \(\sqrt3\)，可取中心点集合为\(\sqrt3\Lambda.\)

于是从原点到点 \(\sqrt3(a+b\omega)\) 的距离平方为\(L^2 = 3|a+b\omega|^2.\)而\(|a+b\omega|^2=(a+b\omega)(a+b\bar\omega)=a^2+ab+b^2.\)令\(n=\dfrac{L^2}{3},\)那么距离为 \(L\) 等价于存在整数 \(a,b\) 使\(a^2+ab+b^2=n.\)

因此 \(\mathrm B(L)\) 就是该方程的整数解 \((a,b)\) 的个数。

为了使用更标准的形式，把它化到 \(x^2+3y^2\) 上。注意恒等变换：\(4(a^2+ab+b^2)=(2a+b)^2+3b^2.\)令$ x=2a+b,y=b,\(得到等价计数：\)x2+3y2=4n.$

所以 \(\mathrm B(L)\) 等于 \(4n\) 被二次型 \(x^2+3y^2\) 表示的方式数（并且每个 \((x,y)\) 对应唯一的 \((a,b)\)）。

在 Eisenstein 整数环 \(\mathbb Z[\omega]\) 里，范数为\(N(a+b\omega)=a^2+ab+b^2.\)

该环是唯一分解整环。整数素数在 \(\mathbb Z[\omega]\) 中的分解规律如下（只用到结论本身即可）：

  • \(p=3\) 是特殊素数：在 \(\mathbb Z[\omega]\) 中平方分裂（本质上不影响计数，因为\(\mathbb Z[\omega]\)有6个单位元\(\{\pm 1,\pm \omega,\pm \omega^2\},\)因此\(3 = (1-\omega)(1-\omega^2)= -\omega^2(1-\omega)^2.\)）。
  • 若普通素数 \(p\equiv 2\pmod 3\)（对奇素数即 \(p\equiv 5\pmod 6\)），则 \(p\) 在 \(\mathbb Z[\omega]\) 中仍为素（不分裂）。
  • 若普通素数 \(p\equiv 1\pmod 3\)（对奇素数即 \(p\equiv 1\pmod 6\)），则 \(p\) 在 \(\mathbb Z[\omega]\) 中分裂成一对共轭素因子。

把“解的个数”理解为：把整数 \(n\) 的素因子在 \(\mathbb Z[\omega]\) 里分配到 \(\alpha\) 与其共轭 \(\bar\alpha\) 之间，使 \(N(\alpha)=n\)，然后再乘上 \(6\) 个单位元带来的旋转对称（蜂巢六方向）。

由此得到经典结论：

  • 若 \(n\) 的分解中存在某个 \(p\equiv 2\pmod 3\)（即 \(p\equiv 5\pmod 6\) 的奇素数，或 \(p=2\)）的指数为奇数，则无法在共轭之间平均分配，故\(\mathrm B(\sqrt{3n})=0.\)
  • 否则（所有 \(p\equiv 2\pmod 3\) 的指数都为偶数），这些素因子在 \(\alpha\) 与 \(\bar\alpha\) 的分配方式是唯一的，不贡献乘数；真正贡献组合数的是分裂素数 \(p\equiv 1\pmod 3\)，即奇素数 \(p\equiv 1\pmod 6\)。

设\(\displaystyle{n = 3^t\cdot \Big(\prod_{p\equiv 1(\bmod 6)} p^{e_p}\Big)\cdot \Big(\prod_{q\equiv 5(\bmod 6)\lor q=2} q^{2f_q}\Big),}\)则\(\mathrm B(\sqrt{3n}) = 6\prod_{p\equiv 1(6)}(e_p+1).\)也可以把\(\displaystyle{x=\prod_{p\equiv 1(\bmod 6)} p^{e_p}}\) 看作 \(n\) 去掉所有“非 \(1\bmod 6\)”素因子后的部分，那么\(\mathrm B(\sqrt{3n}) = 6\tau(x),\)其中 \(\tau\) 是约数个数函数。

可见，\(\mathrm B(L)=450 \iff 6,\tau(x)=450 \iff \tau(x)=75.\)而\(75=3\cdot 5^2.\)若\(\displaystyle{x=\prod p_i^{a_i},\tau(x)=\prod(a_i+1)=75,}\)则满足的指数组合只有三类（忽略排列）：

  1. \((a_1,a_2)=(24,2):x=p^{24}q^2.\)
  2. \((a_1,a_2)=(14,4):x=p^{14}q^4.\)
  3. \((a_1,a_2,a_3)=(4,4,2):x=p^4q^4r^2.\)

这里 \(p,q,r\) 都必须是互不相同且 \(\equiv 1\pmod 6\) 的素数。

注意这三类里指数全是偶数，所以 \(x\) 是完全平方数；再加上其他素因子必须是偶次，得到：\(n=L^2/3\) 一定是 “平方数 \(\times 3^t\)”。

于是\(L^2=3n = 3^{t+1}\cdot w^2,\) 从而 \(L\) 只能落在两类：

  • \(t\) 偶数：\(L=\sqrt3\cdot m\)（\(m\) 为整数）
  • \(t\) 奇数：\(L=3\cdot m\)（\(m\) 为整数）

记\(R_1=\left\lfloor \dfrac{N}{\sqrt3}\right\rfloor, R_2=\left\lfloor \dfrac{N}{3}\right\rfloor.\)我们要数：

  • \(L=\sqrt3m\le 5\times 10^{11}\Rightarrow m\le R_1\)
  • \(L=3m\le 5\times 10^{11}\Rightarrow m\le R_2\)

对给定 \(m\)，有

  • 若 \(L=\sqrt3 m\)，则 \(n=m^2\)
  • 若 \(L=3m\)，则 \(n=3m^2\)

两者的 \(x\)（也就是所有 \(p\equiv 1\pmod6\) 的部分）都来自 \(m^2\) 的那部分，因此完全一样。

设\(m = u\cdot v,\) 其中

  • \(u\) 由 \(p\equiv 1\pmod6\) 的素因子构成（并且指数被上面三种形状严格限定）
  • \(v\) 不含任何 \(p\equiv 1\pmod6\) 的素因子（否则 \(x\) 会多出新的 \(1\bmod6\) 素数，\(\tau(x)\) 就不可能还是 \(75\)）

由于 \(n=m^2\) 或 \(3m^2\)，所有 \(v\) 中的素因子在 \(n\) 里都自动变成偶次，因此不会触发 “\(p\equiv 5\bmod6\) 奇次指数导致 \(\mathrm B=0\)” 的禁忌。

三种骨架（对 \(m\) 的形式）

把 \(x\) 的指数减半（因为 \(m^2\)）：

  1. \(x=p^{24}q^2 \Rightarrow m\) 含 \(u=p^{12}q\)
  2. \(x=p^{14}q^4 \Rightarrow m\) 含 \(u=p^7 q^2\)
  3. \(x=p^4q^4r^2 \Rightarrow m\) 含 \(u=p^2 q^2 r\)

其中 \(p,q,r\) 两两不同且都 \(\equiv 1\pmod 6\) 的素数。

于是，对任意上界 \(R\)（取 \(R_1\) 或 \(R_2\)），答案就是\(\displaystyle{\sum_{u\in \mathcal U} \#\{v\le R/u:(\forall k\in\mathbb{Z}: 6k+1\nmid v)\}.}\)

只要我们预处理一个前缀计数函数\(F(t)=\#\{1\le v\le t: (\forall k\in\mathbb{Z}: 6k+1\nmid v)\},\) 那么每个骨架 \(u\) 的贡献就是 \(F(\lfloor R/u\rfloor)\)。

最终只需要把这两块统计结果相加即可。

代码

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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
from bisect import bisect_right
from array import array
from tools import *

L_MAX = 500_000_000_000

def build_free_prefix(max_v: int, primes_1mod6) -> array:
    """
    F(t) 前缀计数：1..t 中“不含任何 p≡1(mod6) 素因子”的数的个数。
    用筛法标记所有含该类素因子的整数。
    """
    bad = bytearray(max_v + 1)  # bad[i]=1 表示 i 含有某个 1(mod6) 素因子
    for p in primes_1mod6:
        if p > max_v:
            break
        bad[p:max_v + 1:p] = b'\x01' * (((max_v - p) // p) + 1)

    pref = array('I', [0]) * (max_v + 1)
    c = 0
    for i in range(1, max_v + 1):
        if bad[i] == 0:
            c += 1
        pref[i] = c
    return pref

def count_R(R: int, primes_1mod6, F: array) -> int:
    """
    统计 m<=R，使得 B(sqrt(3)*m)=450（等价于 tau(x)=75）。
    三类骨架：
      1) p^12 * q
      2) p^7  * q^2
      3) p^2  * q^2 * r
    其中 p,q,r 为两两不同的 1(mod6) 素数；自由因子 v 计数用 F。
    """
    pr = primes_1mod6
    bsr = bisect_right
    total = 0

    # Case 1: p^12 * q
    for p in pr:
        p12 = p ** 12
        if p12 > R:
            break
        lim_q = R // p12
        idx_q = bsr(pr, lim_q)
        for t in range(idx_q):
            q = pr[t]
            if q == p:
                continue
            base = p12 * q
            total += F[R // base]

    # Case 2: p^7 * q^2
    for p in pr:
        p7 = p ** 7
        if p7 > R:
            break
        lim_q = int(math.isqrt(R // p7))
        idx_q = bsr(pr, lim_q)
        for t in range(idx_q):
            q = pr[t]
            if q == p:
                continue
            base = p7 * q * q
            total += F[R // base]

    # Case 3: p^2 * q^2 * r
    Lpr = len(pr)
    for i in range(Lpr - 1):
        p = pr[i]
        p2 = p * p
        q0 = pr[i + 1]
        base_min = p2 * q0 * q0

        # 最小 r 取 7/13/19 中第一个不等于 p,q0 的（足够用于剪枝）
        if p != 7 and q0 != 7:
            minr = 7
        elif p != 13 and q0 != 13:
            minr = 13
        else:
            minr = 19

        if base_min * minr > R:
            break

        for j in range(i + 1, Lpr):
            q = pr[j]
            base2 = p2 * q * q

            if p != 7 and q != 7:
                minr = 7
            elif p != 13 and q != 13:
                minr = 13
            else:
                minr = 19

            if base2 * minr > R:
                break

            lim_r = R // base2
            idx_r = bsr(pr, lim_r)
            for k in range(idx_r):
                r = pr[k]
                if r == p or r == q:
                    continue
                total += F[R // (base2 * r)]

    return total


R1 = int_sqrt(L_MAX ** 2 // 3)
R2 = L_MAX // 3

min_base = 7 * 7 * 13 * 13
max_prime = R1 // min_base

# primes1 = primes_1mod6_upto(max_prime)
primes1 = [p for p in Factorizer(max_prime).primes if p % 6 == 1]

# max_v 取 case3 的最小骨架 7^2*13^2*19 的商（R1 最大）
max_v = R1 // (7 * 7 * 13 * 13 * 19)
F = build_free_prefix(max_v, primes1)

ans = count_R(R1, primes1, F) + count_R(R2, primes1, F)
print(ans)