Project Euler 295

发表于 更新于 分类于 阅读次数： 本文字数： 3.8k 阅读时长 ≈ 3 分钟

Project Euler 295

题目

Lenticular holes

We call the convex area enclosed by two circles a lenticular hole if:

  • The centres of both circles are on lattice points.
  • The two circles intersect at two distinct lattice points.
  • The interior of the convex area enclosed by both circles does not contain any lattice points.

Consider the circles:

\(\begin{aligned} &C_0: x^2+y^2=25\\ &C_1: (x+4)^2+(y-4)^2=1\\ &C_2: (x-12)^2+(y-4)^2=65 \end{aligned}\)

The circles \(C_0, C_1\) and \(C_2\) are drawn in the picture below.

\(C_0\) and \(C_1\) form a lenticular hole, as well as \(C_0\) and \(C_2\).

We call an ordered pair of positive real numbers \((r_1, r_2)\) a lenticular pair if there exist two circles with radii \(r_1\) and \(r_2\) that form a lenticular hole.

We can verify that \((1, 5)\) and \((5, \sqrt{65})\) are the lenticular pairs of the example above.

Let \(L(N)\) be the number of distinct lenticular pairs \((r_1, r_2)\) for which \(0 < r_1 \le r_2 \le N\).

We can verify that \(L(10) = 30\) and \(L(100) = 3442\).

Find \(L(100 000)\).

解决方案

两圆相交于两格点 \(A,B\)，透镜是两圆盘交集形成的凸“月牙”区域。平移不改变格点结构与“内部是否含格点”，故不妨令

  • \(A=(0,0)\)
  • \(B=(x,y)\)，其中 \(x,y\in\mathbb Z\) 且不全为 \(0\)。

透镜内部不能含格点，尤其线段 \(AB\) 上不能出现除端点外的格点；否则这些点在两圆盘交集中，必在透镜内部，违背定义。因此必须有：\(\gcd(x,y)=1.\)

此外，过 \(A,B\) 的圆心必须在 \(AB\) 的垂直平分线上。\(AB\) 的中点为 \(\left(\dfrac x2,\dfrac y2\right)\)，垂直方向向量可取 \((y,-x)\)，所以垂直平分线为

\[M(t)=\left(\frac x2,\frac y2\right)+t(y,-x), t\in\mathbb R.\]

如果要求圆心 \(M(t)\) 也是格点。由于 \(\gcd(x,y)=1\)，那么 \(x,y\) 不可能同为偶数；要让 \(\left(\dfrac {x}{2},\dfrac {y}{2}\right)\) 与 \(t(y,-x)\) 合起来变成整数坐标，唯一可行的是 \(x,y\) 同为奇数。原因在于，如果\(x\)是偶数，那么\(y\)必须是奇数，由于\(2\mid ty\)，因此\(t\)必须是偶数，但是\(\dfrac{y}{2}-tx\) 此时不可能是整数。因此，\(x,y\) 同为奇数。

因此我们只需枚举 互素的奇数对 \((x,y)\)（利用对称性可取 \(1\le x<y\)，并单独处理 \((1,1)\)）。

当 \(x,y\) 为奇数时，令 \(t=k+\dfrac{1}{2}\)（\(k\in\mathbb Z\)），则

\[ M_k=\left(\frac{x}{2},\frac{y}{2}\right)+\left(k+\dfrac{1}{2}\right)(y,-x)=\left(\frac{x+y}{2}+ky,\frac{y-x}{2}-kx\right) \]

给出了垂直平分线上的全部格点圆心。

因为圆过原点 \(A=(0,0)\)，半径平方就是

\[r_k^2 = |M_k|^2=\left(\frac{x+y}{2}+ky\right)^2+\left(\frac{y-x}{2}-kx\right)^2=\frac{x^2+y^2}{2}(2k^2+2k+1)\]

记\(b=x^2+y^2\)，注意到 \(x,y\) 为奇数时 \(b\equiv 2\pmod 4\)，所以 \(\dfrac{b}{2}\in\mathbb Z\)，从而\(r_k^2=\dfrac{b}{2}(2k^2+2k+1)\in\mathbb Z\)，于是本题可用“半径平方”这个整数来表示半径，去重也在整数域完成。

对于固定的 \((x,y)\)，当 \(|k|\) 增大时，圆心 \(M_k\) 远离弦 \(AB\)，圆弧会更贴近弦，透镜那一侧的“圆弓形”单调变薄；因此存在最小 \(k_{\min}\)，使得当 \(k\ge k_{\min}\) 时，该圆弓形不再包含任何格点，并且所有更大的 \(k\) 都安全。

接下来的任务是求 \(k_{\min}\)。为了保证\(AB\)下方的每个格点都不在内部里，因此对每个整数 \(u\in\{0,1,\dots,x-1\}\)，令\(v(u)=\left\lfloor \dfrac{uy}{x}\right\rfloor+1,P_u=(u,v(u)).\)那么\(P_u\)是固定 \(u\) 的那一列里，它是离直线 \(AB\) 最近的格点，因为在\(P_u\)正下方的格点更早退出。

圆心为 \(M_k=(M_x,M_y)\) 且过原点的圆满足

\[ (X-M_x)^2+(Y-M_y)^2 = M_x^2+M_y^2\iff X^2+Y^2-2(M_xX+M_yY)=0. \]

对格点 \(P=(u,v)\)，定义\(F_k(P)=u^2+v^2-2(M_xu+M_yv).\)则 \(F_k(P)<0\) 表示 \(P\) 在圆内，\(F_k(P)=0\) 在圆上，\(F_k(P)>0\) 在圆外。把 \(M_k\) 代入并整理可得：

\[ F_k(P)=(u^2+v^2-u(x+y)+v(x-y))+2k(vx-uy). \]

由于\(P\)不在圆内（允许在圆上），即 \(F_k(P)\ge 0\)，并且可见，\(P\)在\(\overrightarrow{AB}\)左侧，因此\(\overrightarrow{AB}\times \overrightarrow{AP}=vx-uy>0\)，因此这时有\(\displaystyle{k\ge \frac{-(u^2+v^2+v(x-y)-u(x+y))}{2(vx-uy)}}\)

因此，对于任意\(u\in\{0,\dots,x-1\}\)，有\(k_{\min}(u)=\left\lceil\dfrac{u(x+y)-v(u)(x-y)-u^2-v(u)^2}{2(vx-uy)}\right\rceil.\)

最终有：

\[ k_{\min}=\max_{u=0}^{x-1} \{k_{\min}(u)\}. \]

此外，由于\(r_k^2=\dfrac {b}{2}(2k^2+2k+1)\le N^2\).因此\(k_{\max}\)只需要满足\(k_{\max}^2+2k_{\max}+1 \le \left\lfloor\dfrac{2N^2}{b}\right\rfloor\)即可.

对每个合法 \((x,y)\)，可用的 \(k\) 区间为\(k_{\min}\le k\le k_{\max}\)。设区间长度\(f = k_{\max}-k_{\min}+1.\)

则该 \((x,y)\) 产生的半径平方集合大小为 \(f\)，任取两圆半径（允许相等，表示两圆同半径但圆心在弦两侧）得到的有序对数为\(\dbinom{f+1}{2}\)。把所有 \((x,y)\) 的该值累加得到一个“带重复”的总数。

记每个\((x,y)\)对应的半径平方集合为 \(S_{x,y}\)，我们要数的是所有\(S_{x,y}=\{(a,b): (a,b)\in S_{x,y}, a\le b\}\)的并集大小。重复只会由“某些半径平方同时出现在多个 \(S_{x,y}\)”引起。

做法：

  1. 第一遍枚举所有 \((x,y,k)\) 生成所有 \(r^2\)，找出出现超过一次的 \(r^2\)。
  2. 第二遍再次枚举，把每个重复 \(r^2\) 关联到它出现过的 \((x,y)\) 列表 \(I(r^2)\)。
  3. 对每个 \(I(r^2)\)，把它的所有子集 \(T\)（\(|T|\ge 2\)）的计数 \(c_T\) 加一，其中 \(c_T\) 表示\(\displaystyle{c_T = \left|\bigcap_{(x,y)\in T} S_{x,y}\right|.}\)
  4. 用容斥原理修正：\(\displaystyle{\sum_{|T|\ge 2}(-1)^{|T|+1}\binom{c_T+1}{2}}\)加回到初始累加值即可得到去重后的 \(L(N)\)。

代码

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
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
from tools import int_sqrt, gcd
from collections import defaultdict


def ceil_div(a: int, b: int) -> int:
    """
    向上取整的整除：ceil(a / b)，要求 b > 0。
    这里用 -((-a)//b) 的技巧保证 a 为负数时也正确。

    例：
      ceil_div(5,2)=3
      ceil_div(-5,2)=-2   (因为 -2 <= -2.5 < -1)
    """
    return -((-a) // b)


def rr_from_b_k(b: int, k: int) -> int:
    """
    给定 chord 的 b = x^2 + y^2，以及整数 k，
    计算对应圆的半径平方 r^2（注意 r 可能是无理数，但 r^2 是整数）：

        r^2 = (b/2) * (2k^2 + 2k + 1)

    在本题中 x,y 必为奇数 => b ≡ 2 (mod 4)，因此 b/2 为整数，
    故 r^2 是整数，可直接用整型存储、去重与比较。
    """
    return (b * (2 * k * k + 2 * k + 1)) // 2


def get_kmin(x: int, y: int) -> int:
    """
    计算给定 chord 向量 (x,y) 的最小 k（k_min），使得：
      对所有落在 chord 线段 AB 的同侧、且“最靠近 AB 的格点”，圆不包含它们在内部。
    从而保证透镜内部不含格点。

    几何背景（对应题解的“严格版本”）：
      - 固定两交点 A=(0,0), B=(x,y)，圆心在 AB 的垂直平分线上
      - 当我们用整数 k 枚举垂直平分线上的格点圆心 M_k 时，
        k 越大圆心离 AB 越远，对侧圆弓形越薄，越不容易包进格点
      - 因此存在最小 k_min，使得对所有 k >= k_min 都“安全”

    这里采用的严格计算方法：
      对每一列 u = 0..x-1，取该列中“刚好在直线 AB 上方的最邻近格点”
        P_u = (u, floor(u*y/x) + 1)
      对每个 P_u，写出点幂/代入展开得到不等式：
        F_k(P_u) = C(P_u) + 2k * Delta(P_u) >= 0
      其中
        Delta(P) = (AB x AP) = v*x - u*y  > 0  (保证 P 在 AB 左侧/上方那一侧)
        C(P)     = u^2 + v^2 - u(x+y) + v(x-y)

      因此对每个 P_u，
        k >= ceil( -C(P_u) / (2*Delta(P_u)) )
      取所有列 u 的最大值，即得到 k_min。

    为什么只用这些 P_u 就够：
      它们构成“紧贴 AB 的那层格点边界”，任何更远的格点需要更厚的圆弓形才包得住，
      所以若这层边界都在圆外/圆上，则更远的格点也不会在圆内。
    """
    km = 0
    xp = x + y  # 预先缓存 x+y
    xm = x - y  # 预先缓存 x-y

    # 枚举 u=0..x-1，构造每一列刚好在 AB 上方的最近格点 P_u
    for u in range(x):
        v = (u * y) // x + 1  # floor(u*y/x)+1，确保点在直线 AB 上方（同侧）
        # C = u^2+v^2 - u(x+y) + v(x-y)
        # 不等式：C + 2k*Delta >= 0  => k >= -C/(2*Delta)
        C = u * u + v * v - u * xp + v * xm
        delta = v * x - u * y  # Delta = v*x - u*y，按构造必为正
        k = ceil_div(-C, 2 * delta)
        if k > km:
            km = k

    return km


def get_kmax(b: int, N2: int) -> int:
    """
    计算给定 b=x^2+y^2 时，满足 r^2 <= N^2 的最大 k（k_max）。

    r^2 = b/2 * (2k^2+2k+1) <= N^2
    等价于：
      2k^2+2k+1 <= floor(2N^2 / b) = s

    解二次不等式可得：
      k <= (sqrt(2s-1)-1)/2

    为避免浮点误差，这里用整数 sqrt：
      t = isqrt(2s-1)
      k = (t-1)//2
    再用 while 做一两步安全修正，确保 rr_from_b_k(b,k) <= N2 且 k+1 不再满足。
    """
    s = (2 * N2) // b
    if s <= 0:
        return -1

    # 2k^2+2k+1 <= s  =>  k <= (sqrt(2s-1)-1)/2
    t = int_sqrt(2 * s - 1)
    k = (t - 1) // 2

    # 安全修正（理论上不多次循环，通常 0~1 次）
    while k >= 0 and rr_from_b_k(b, k) > N2:
        k -= 1
    while rr_from_b_k(b, k + 1) <= N2:
        k += 1
    return k


def gen_chords(N: int):
    """
    枚举所有可能产生至少一个可行半径的 chord 向量 (x,y)。

    结论（来自几何/格点条件）：
      1) 交点差向量 (x,y) 必须为奇奇（x,y 均为奇数），否则垂直平分线上没有格点圆心
      2) gcd(x,y)=1，否则 AB 上有内部格点，透镜必含格点

    利用对称性：
      chord (x,y) 与 (y,x) 其实是旋转/对称等价，因此只取 1 <= x < y，
      另外特殊处理 (1,1)（它不满足 x<y）。

    终止条件：
      对每个 y（奇数递增），先检查 x=1 是否还存在可行 k 区间；
      若 (1,y) 已无可行半径，则更大的 y（以及同 y 的更大 x）也不会有（本题经验上单调），可 break。
    """
    N2 = N * N
    chords = []

    # 单独加入 (1,1)
    x = y = 1
    b = 2
    km = get_kmin(x, y)
    kx = get_kmax(b, N2)
    if kx >= km:
        chords.append((x, y, b, km, kx))

    # 枚举 y 为奇数：3,5,7,...
    y = 3
    while True:
        # 先用 (1,y) 做“是否可以停止”的判定（经验上单调足够）
        btest = 1 + y * y
        km_test = get_kmin(1, y)
        kx_test = get_kmax(btest, N2)
        if kx_test < km_test:
            break

        # 对固定 y，枚举所有奇数 x < y
        for x in range(1, y, 2):
            if gcd(x, y) != 1:
                continue

            b = x * x + y * y
            kx = get_kmax(b, N2)
            if kx < 0:
                continue

            km = get_kmin(x, y)
            if kx >= km:
                # 记录该 chord 对应的 b 以及可行 k 区间 [km,kx]
                chords.append((x, y, b, km, kx))

        y += 2

    return chords


def add_all_subsets_ge2(ids, mp):
    """
    给定一个 chord-id 列表 ids（已去重并排序），
    对其所有大小 >=2 的子集 T，令 mp[T] += 1。

    这里的语义用于包含-排除（PIE）：
      - 某个半径平方 rr 可能出现在多个 chord 的可行集合 S_{x,y} 中
      - 对这个 rr，我们得到它出现过的 chord-id 集合 I(rr)
      - 那么对于 I(rr) 的任意子集 T，rr 都属于所有 chord in T 的交集
        => 交集计数 c_T 需要 +1

    最终 mp[T] 存的就是 c_T = | ∩_{i in T} S_i |（仅统计“重复 rr”带来的部分）。
    """
    m = len(ids)
    if m < 2:
        return

    cur = []

    def dfs(start):
        for i in range(start, m):
            cur.append(ids[i])
            if len(cur) >= 2:
                t = tuple(cur)
                mp[t] += 1
            dfs(i + 1)
            cur.pop()

    dfs(0)


def L(N: int) -> int:
    """
    主函数：计算 L(N)，即所有不同 lenticular pair (r1,r2)（r1<=r2<=N）的个数。

    总体思路：
      - 对每个 chord (x,y)，它对应一段可行 k 区间 [kmin,kmax]，
        生成一组半径平方 S_{x,y}
      - 在该 chord 内，任取两条半径（允许相同）都能组成一个“半径对”
        => 贡献 |S|*(|S|+1)/2
      - 但不同 chord 的集合可能有交集（相同 rr 出现多次），导致 (r1,r2) 重复计数
      - 通过“两遍扫描 + PIE 修正”精确去重

    两遍扫描：
      Pass1:
        - 计算 baseline = sum over chord of f*(f+1)/2
        - 同时找出哪些 rr（半径平方）在不同 chord 中出现过不止一次 => multiples 集合
      Pass2:
        - 仅对 rr in multiples，记录它出现在哪些 chord（形成 I(rr)）
        - 对每个 rr 的 chord 列表 I(rr)，对其所有子集 T (|T|>=2) 做 c_T++
      PIE：
        - 对每个 T，有 c_T 个共同 rr，它们内部产生的半径对数为 C(c_T+1,2)
        - 并集计数的包含-排除修正项为：
            + sum_{|T| odd}  C(c_T+1,2)
            - sum_{|T| even} C(c_T+1,2)
        - 把该修正项加到 baseline 上，即得到去重后的答案
    """
    chords = gen_chords(N)

    # Pass 1: baseline 累加 + 找到重复出现的 rr
    unique = set()     # 所有出现过的 rr
    multiples = set()  # 出现次数 >=2 的 rr
    ans = 0

    for x, y, b, km, kx in chords:
        f = kx - km + 1
        ans += f * (f + 1) // 2  # chord 内的半径对数量（允许相同半径）

        # 逐个生成 rr(k)，并检测是否重复
        rr = rr_from_b_k(b, km)
        for k in range(km, kx + 1):
            if rr in unique:
                multiples.add(rr)
            else:
                unique.add(rr)

            # 递推：rr(k+1) = rr(k) + 2*b*(k+1)
            # 因为：
            #   2k^2+2k+1 的差分为 4(k+1)
            #   rr = b/2 * (...)，所以差分为 b/2 * 4(k+1) = 2*b*(k+1)
            rr += 2 * b * (k + 1)

    # Pass 2: 对每个重复 rr，收集它出现的 chord-id 列表
    in_map = defaultdict(list)
    for x, y, b, km, kx in chords:
        cid = (x << 16) | y  # 一个可哈希/可比较的 chord 编码
        rr = rr_from_b_k(b, km)
        for k in range(km, kx + 1):
            if rr in multiples:
                in_map[rr].append(cid)
            rr += 2 * b * (k + 1)

    # 对每个 rr 的 chord 列表 I(rr)，对所有子集 T (|T|>=2) 令 c_T++
    count = defaultdict(int)
    for lst in in_map.values():
        ids = sorted(set(lst))
        add_all_subsets_ge2(ids, count)

    # Inclusion-Exclusion 修正：
    # 对每个子集 T，有 c = c_T 个共同 rr，则它们内部的半径对数为 c*(c+1)/2。
    # 若 |T| 为奇数 => 加；若 |T| 为偶数 => 减。
    for subset, c in count.items():
        add = c * (c + 1) // 2
        if len(subset) % 2 == 0:
            ans -= add
        else:
            ans += add

    return ans


if __name__ == "__main__":
    # 题目给出的校验值
    print(L(10))        # 30
    print(L(100))       # 3442

    # 目标
    print(L(100000))    # 4884650818