Project Euler 780

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

Project Euler 780

题目

Toriangulations

For positive real numbers \(a,b\), an \(a\times b\) torus is a rectangle of width \(a\) and height \(b\), with left and right sides identified, as well as top and bottom sides identified. In other words, when tracing a path on the rectangle, reaching an edge results in “wrapping round” to the corresponding point on the opposite edge.

A tiling of a torus is a way to dissect it into equilateral triangles of edge length \(1\). For example, the following three diagrams illustrate respectively a \(1\times \dfrac{\sqrt{3}}{2}\) torus with two triangles, a \(\sqrt{3}\times 1\) torus with four triangles, and an approximately \(2.8432\times 2.1322\) torus with fourteen triangles:

Two tilings of an \(a\times b\) torus are called equivalent if it is possible to obtain one from the other by continuously moving all triangles so that no gaps appear and no triangles overlap at any stage during the movement. For example, the animation below shows an equivalence between two tilings:

Let \(F(n)\) be the total number of non-equivalent tilings of all possible tori with exactly \(n\) triangles. For example, \(F(6)=8\), with the eight non-equivalent tilings with six triangles listed below:

Let \(G(N)=\sum_{n=1}^N F(n)\). You are given that \(G(6)=14, G(100)=8090\), and \(G(10^5)\equiv 645124048 \pmod{1\,000\,000\,007}\).

Find \(G(10^9)\). Give your answer modulo \(1\,000\,000\,007\).

解决方案

把一个 \(a\times b\) 的 torus 看成矩形基本域，两组对边分别粘合。若它被边长 \(1\) 的正三角形铺满，那么把基本域在平面上按周期复制，就得到一个整张平面的边长 \(1\) 等边三角形铺法（周期性铺法）。

因此，可见任意这样的平面铺法都可以看成由若干条平行条带组成，每条条带宽度固定为 \(\sqrt{3} / 2\)，内部是上下交替的单位等边三角形链，并且一般情况下条带方向是唯一的；但当铺法恰好是三角晶格（六角/蜂窝意义下的那类）时，会有 \(3\) 个等价方向可选。

这立即带来两个结论：每条条带里上下三角形数量相同，所以总三角形数 \(n\) 必须是偶数；我们可以先按条带方向计数，再把三角晶格导致的三重计数修正掉。因此奇数 \(n\) 时 \(F(n)=0\)。

现在假设\(n\)是偶数。设 \(n=2m\)。设 torus 基本域边长为 \(X\)（水平）和 \(Y\)（竖直），面积满足\(XY = \dfrac{\sqrt 3}{4}n = \dfrac{\sqrt 3}{2}m.\)

固定某个条带方向。沿该方向绕 torus 一圈后会闭合，因此该方向相对矩形的斜率必须是有理绕行数：设横向绕行 \(p\) 次、纵向绕行 \(q\) 次（\(p,q\) 互素，允许其中一个为 \(0\)）；若共有 \(k\) 条平行条带（同方向的不同平移类），则每条条带长度恰好是 \(m/k\)（因为每条条带含 \(m/k\) 个上三角与 \(m/k\) 个下三角，总 \(2m/k\) 个三角，沿长度方向每个三角贡献底边长 \(1\)）。于是有长度方程：

\[ (pX)^2 + (qY)^2 = \left(\frac{m}{k}\right)^2. \]

同时面积方程为上面的 \(XY = (\sqrt3/2)m\)。把\(A = pX, B = qY\)代入，则\(A^2+B^2=\left(\dfrac{m}{k}\right)^2, AB = pq\cdot XY = pq\cdot\dfrac{\sqrt3}{2}m.\)因此\((A^2-B^2)^2=(A^2+B^2)^2-4A^2B^2=\left(\dfrac{m}{k}\right)^4-3p^2q^2m^2.\)要有实数解，必须且只需右边非负，即\(\left(\dfrac{m}{k}\right)^4 \ge 3p^2q^2m^2\iff m \ge \sqrt3k^2pq.\)

因此，若 \(p=0\) 或 \(q=0\)（条带平行于坐标轴），给定 \((k,p,q)\) 时只有 \(1\) 个正解 \((X,Y)\)；若 \(p,q>0\)，由 \(A^2+B^2\) 与 \(|A^2-B^2|\) 可得两种解（交换大/小分配），所以有 \(2\) 个正解。此外还要计入正负斜率对 torus 铺法计数的区别，再乘 \(2\)。

于是轴向条带（\(p=0\) 或 \(q=0\)）每组贡献 \(2\)；非轴向条带（\(p,q>0\)）每组贡献 \(4\)。

令\(u=kp, v=kq.\)则 \(\gcd(u,v)=k\)，且 \(p,q\) 互素。非轴向时 \(u,v\ge1\)。条件\(m \ge \sqrt3k^2pq\)变成\(m \ge \sqrt3uv.\)但 \(m\) 还必须是 \(k=\gcd(u,v)\) 的倍数（因为 \(m=kt\)）。所以在 \(1\le n\le N\)（即 \(1\le m\le \lfloor N/2\rfloor\)）范围内，对固定 \((u,v)\) 的可选 \(m\) 数量是\(\left\lfloor \dfrac{N}{2\gcd(u,v)} \right\rfloor-\left\lfloor \sqrt3\cdot \dfrac{uv}{\gcd(u,v)} \right\rfloor.\)

设\(M=\left\lfloor\dfrac N2\right\rfloor,L=\left\lfloor \dfrac{M}{\sqrt3}\right\rfloor=\left\lfloor\dfrac{N}{2\sqrt3}\right\rfloor.\)则所有条带方向计数（暂未修正三重计数）为

\[ G_{\text{strip}}(N) =2D(M)+ 4\sum_{uv\le L} \left( \left\lfloor \frac{N}{2\gcd(u,v)} \right\rfloor- \left\lfloor \sqrt3\cdot \frac{uv}{\gcd(u,v)} \right\rfloor \right), \]

其中\(\displaystyle D(x)=\sum_{t=1}^{x}\left\lfloor \frac{x}{t}\right\rfloor.\)上式第一项 \(2D(M)\) 就是轴向条带的贡献，第二项是非轴向条带的贡献。

不过上面的主计数 \(G_{\text{strip}}(N)\) 是按条带方向来数铺法的。对于一般铺法，这样计数没有问题；但当铺法的所有顶点都落在三角晶格点上时，会出现一个特殊现象：同一个铺法可以沿三角晶格的三组轴方向中的任意一组切成条带，因此在 \(G_{\text{strip}}(N)\) 中会被计数 \(3\) 次，而在 \(F(n)\) 中它本应只计 \(1\) 次。

因此，这一类三角晶格铺法需要做统一修正：每个铺法多算了 \(2\) 次，所以最终要减去这类铺法总数的 \(2\) 倍。下面我们直接参数化并计数这类铺法。

为了计算这部分修正项，下面采用等价的晶格矩形原始块参数化。

三角晶格上的（按旋转归类后的）原始矩形可由整数对 \((u,v)\) 参数化，并对应到基本面积参数\(Q(u,v)=u^2+uv+v^2.\)满足原始性（即去掉可进一步分解、以及会导致三方向重复退化的情形）的条件可写成：\(u>v,\gcd(u,v)=1,2u+v \not\equiv 0 \pmod 3\)

此外，还要单独加上最基本的轴向原始块（对应 \(Q=1\) 的那一族）。

令\(X=\left\lfloor\frac{N}{4}\right\rfloor.\)定义 \(H(X)\) 为这些原始晶格矩形所生成的 torus 铺法前缀计数（并且把横向/纵向堆叠产生的所有倍数都计入）：

\[ H(X)=D(X)+2\sum_{\substack{u>v\ge1\\ \gcd(u,v)=1\\ 2u+v\not\equiv0\pmod 3\\ Q(u,v)\le X}} D\left(\left\lfloor\frac{X}{Q(u,v)}\right\rfloor\right). \]

在这里，\(D(X)\) 是轴向原始块族的贡献，以及非轴向部分前面的系数 \(2\) 来自旋转/对称带来的两种取向。

由此可得，三角晶格型 torus 铺法的总前缀数为 \(2H(X)\)（每个原始矩形对应两种取向）。而这部分在 \(G_{\text{strip}}(N)\) 中被数了 \(3\) 次，所以相对于正确计数需要额外减去\(2\times (2H(X)) = 4H(X).\)

因此最终修正后的结果为

\[ G(N)=G_{\text{strip}}(N)-4H\left(\left\lfloor\frac N4\right\rfloor\right). \]

直接双重枚举 \(uv\le L\) 会很慢，我们需要两个关键优化。

在非轴向条带项中，固定某个参数（例如固定 \(v\)）后，会出现形如“统计与某个整数互素的个数”的子问题。

设需要统计\(\#\{t\le T:\gcd(t,m)=1\},\)则可用 Möbius 反演写成\(\displaystyle\#\{t\le T:\gcd(t,m)=1\}=\sum_{q\mid m}\mu(q)\left\lfloor\frac{T}{q}\right\rfloor.\)

这样一来，原本带有 \(\gcd\) 约束的计数，就转化为对 \(m\) 的因子（或无平方因子）的求和。在实现中，只需预处理：Möbius 函数 \(\mu(\cdot)\)；每个整数的因子表（或无平方因子表）；就可以把互素计数压缩到 \(O\left(\sum \tau(v)\right)\) 量级处理（这里 \(\tau\) 是约数个数函数）。

这一步主要用于条带主项中由 \(\gcd\) 导致的计数拆分，是后续 Beatty 求和的前提。

可见，条带主项的另一部分会反复出现形如\(\sum \left\lfloor c\sqrt3\cdot t \right\rfloor\)的求和（其中 \(c\) 为整数或由当前枚举参数确定的常数）。这正是 Beatty 型前缀和。定义\(\displaystyle T(n)=\frac{n(n+1)}2,B(\alpha,n)=\sum_{k=1}^{n}\lfloor k\alpha\rfloor.\)对无理数 \(\alpha>1\)，递归分两种情况。大致过程如下：

当 \(\alpha\ge 2\) 时，先拆整数部分。设\(\alpha=r+\beta, r=\lfloor \alpha\rfloor, 0<\beta<1.\)则\(\lfloor k\alpha\rfloor=\lfloor kr+k\beta\rfloor=kr+\lfloor k\beta\rfloor,\)所以\(\displaystyle B(\alpha,n)=rT(n)+\sum_{k=1}^{n}\lfloor k\beta\rfloor.\)

为了继续使用参数 \(>1\) 的递归形式，通常把 \(\beta\) 改写成 \(1+(\beta-1)\) 的处理方式。因此，更直接地写成（取 \(i=\lfloor\alpha\rfloor-1\)，使 \(\alpha-i\in(1,2)\)）：

\[B(\alpha,n)=iT(n)+B(\alpha-i,n).\]

当 \(1<\alpha<2\) 时，用 Beatty 互补递归。设\(m=\lfloor (\alpha-1)n\rfloor, t=\lfloor \alpha n\rfloor,\)并令\(\alpha'=\dfrac{\alpha}{\alpha-1}.\)则有递归

\[ B(\alpha,n)=T(t)-B(\alpha',m). \]

其中 \(t=\lfloor \alpha n\rfloor\)，\(m=\lfloor(\alpha-1)n\rfloor=t-n\)。

这一步会把参数从 \((1,2)\) 映射到一个新的大于 \(2\) 的无理数，再配合上一条拆整数部分继续递归。由此我们方便地计算出了这类无理数向下取整前缀和。

其余部分（包括三角晶格修正项）在实现中按前述参数化直接计算即可.

代码

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
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
from typing import List, Tuple

from tools import Factorizer, gcd, int_sqrt
N = 10**9
MOD = 1_000_000_007

class Fact2(Factorizer):
    def __init__(self, lim: int):
        super().__init__(max(lim, 1))
        self.all_divs: List[List[int]] = [[] for _ in range(self.maxN + 1)]
        self.sf_mu_divs: List[List[Tuple[int, int]]] = [[] for _ in range(self.maxN + 1)]
        self.all_divs[0] = []
        self.sf_mu_divs[0] = []
        self.all_divs[1] = [1]
        self.sf_mu_divs[1] = [(1, 1)]
        for x in range(2, self.maxN + 1):
            fac = self.factor_small(x)
            self.all_divs[x] = self._build_divs_from_factor(fac)
            self.sf_mu_divs[x] = self._build_sqfree_mu_divs_from_factor(fac)

    def _build_divs_from_factor(self, fac: List[Tuple[int, int]]) -> List[int]:
        divs = [1]
        for p, e in fac:
            nxt = []
            pe = 1
            for _ in range(e + 1):
                for d in divs:
                    nxt.append(d * pe)
                pe *= p
            divs = nxt
        return divs

    def _build_sqfree_mu_divs_from_factor(self, fac: List[Tuple[int, int]]) -> List[Tuple[int, int]]:
        res = [(1, 1)]
        for p, _ in fac:
            res += [(d * p, -sgn) for d, sgn in res]
        return res

    def divisors(self, x: int) -> List[int]:
        return self.all_divs[x]

    def squarefree_mu_divs(self, x: int) -> List[Tuple[int, int]]:
        return self.sf_mu_divs[x]


def floor_sqrt3_mul(n: int) -> int:
    return int_sqrt(3 * n * n)


def floor_div_sqrt3(n: int, m: int) -> int:
    nn = n * n
    den = 3 * m * m
    t = int_sqrt(nn // den)
    while den * (t + 1) * (t + 1) <= nn:
        t += 1
    while den * t * t > nn:
        t -= 1
    return t


def divisor_summatory(n: int) -> int:
    if n <= 0:
        return 0
    res = 0
    i = 1
    while i <= n:
        q = n // i
        j = n // q
        res += q * (j - i + 1)
        i = j + 1
    return res


def norm_abc(a: int, b: int, c: int):
    if c < 0:
        a, b, c = -a, -b, -c
    g = gcd(gcd(abs(a), abs(b)), abs(c))
    if g > 1:
        a //= g
        b //= g
        c //= g
    return a, b, c


def floor_qsqrt3(a: int, b: int, c: int) -> int:
    if b == 0:
        return a // c
    x = (a + floor_sqrt3_mul(abs(b))) // c
    bb3 = 3 * b * b
    while True:
        y = (x + 1) * c - a
        if y <= 0 or y * y <= bb3:
            x += 1
        else:
            break
    while True:
        y = x * c - a
        if y <= 0 or y * y <= bb3:
            break
        x -= 1
    return x


def alpha_floor(alpha):
    a, b, c = alpha
    return floor_qsqrt3(a, b, c)


def alpha_mul_floor(alpha, n: int):
    a, b, c = alpha
    return floor_qsqrt3(a * n, b * n, c)


def alpha_sub_int(alpha, k: int):
    a, b, c = alpha
    return norm_abc(a - k * c, b, c)


def alpha_div_alpha_minus1(alpha):
    a, b, c = alpha
    ac = a - c
    A = a * ac - 3 * b * b
    B = -b * c
    D = ac * ac - 3 * b * b
    if D < 0:
        A, B, D = -A, -B, -D
    if B < 0:
        A, B = -A, -B
    return norm_abc(A, B, D)


def tri(n: int) -> int:
    return n * (n + 1) // 2


def beatty_sum(alpha, n: int) -> int:
    res = 0
    sign = 1
    while n > 0:
        f = alpha_floor(alpha)
        if f > 1:
            res += sign * (f - 1) * tri(n)
            alpha = alpha_sub_int(alpha, f - 1)
        m = alpha_mul_floor(alpha, n)
        res += sign * tri(m)
        n = m - n
        if n <= 0:
            break
        alpha = alpha_div_alpha_minus1(alpha)
        sign = -sign
    return res


def beatty_sqrt3(c: int, n: int) -> int:
    if n <= 0:
        return 0
    return beatty_sum((0, c, 1), n)


def strip_sums(N: int, V: int, L: int, nt: Fact2):
    s1 = 0
    s2 = 0
    for v in range(1, V + 1):
        Umax = L // v
        for d in nt.divisors(v):
            m = v // d
            hi = Umax // d
            if hi < m:
                continue
            w = N // (2 * d)
            lo = m - 1
            cnt = 0
            sm = 0
            for q, muq in nt.squarefree_mu_divs(m):
                cnt += muq * (hi // q - lo // q)
                hiq = hi // q
                loq = lo // q
                if hiq:
                    sm += muq * (beatty_sqrt3(v * q, hiq) - beatty_sqrt3(v * q, loq))
            s1 += w * cnt
            s2 += sm
    d1 = 0
    d2 = 0
    for v in range(1, V + 1):
        d1 += N // (2 * v)
        d2 += floor_sqrt3_mul(v)
    return 2 * s1 - d1, 2 * s2 - d2


def count_mod3(lo: int, hi: int, r: int) -> int:
    if hi < lo:
        return 0
    first = lo + (r - lo % 3) % 3
    if first > hi:
        return 0
    return (hi - first) // 3 + 1


def hex_hsum(X: int, nt: Fact2):
    if X <= 0:
        return 0
    dcache = {}

    def D(n: int) -> int:
        if n <= 0:
            return 0
        if n not in dcache:
            dcache[n] = divisor_summatory(n)
        return dcache[n]

    V = int_sqrt(X)
    extra = 0
    for v in range(1, V + 1):
        disc0 = 4 * X - 3 * v * v
        if disc0 <= 0:
            break
        Umax = (-v + int_sqrt(disc0)) // 2
        u = v + 1
        sfv = nt.squarefree_mu_divs(v)
        vmod = v % 3
        while u <= Umax:
            t = u * u + u * v + v * v
            q = X // t
            if q == 0:
                break
            T = X // q
            uhi = (-v + int_sqrt(4 * T - 3 * v * v)) // 2
            if uhi > Umax:
                uhi = Umax
            lo = u - 1
            total = 0
            for d, mud in sfv:
                total += mud * (uhi // d - lo // d)
            if vmod:
                bad = 0
                r = vmod
                for d, mud in sfv:
                    dm3 = d % 3
                    inv = 1 if dm3 == 1 else 2
                    tlo = (u + d - 1) // d
                    thi = uhi // d
                    rr = (r * inv) % 3
                    bad += mud * count_mod3(tlo, thi, rr)
                total -= bad
            if total:
                extra += total * D(q)
            u = uhi + 1
    return D(X) + 2 * extra


def solve(N: int) -> int:
    if N <= 0:
        return 0
    M = N // 2
    L = floor_div_sqrt3(M, 1) if M > 0 else 0
    V1 = int_sqrt(L) if L > 0 else 0
    X = N // 4
    V2 = int_sqrt(X) if X > 0 else 0
    nt = Fact2(max(V1, V2))
    base = 2 * divisor_summatory(M)
    s1, s2 = strip_sums(N, V1, L, nt)
    h = hex_hsum(X, nt)
    return (base + 4 * (s1 - s2) - 4 * h) % MOD


print(solve(N))