Project Euler 372

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

Project Euler 372

题目

Pencils of rays

Let \(R(M, N)\) be the number of lattice points \((x, y)\) which satisfy \(M<x\le N, M<y\le N\) and \(\left\lfloor \dfrac{y^2}{x^2} \right\rfloor\) is odd.

We can verify that \(R(0, 100) = 3019\) and \(R(100, 10000) = 29750422\).

Find \(R(2\cdot 10^6, 10^9)\).

Note : \(\lfloor x \rfloor\) represents the floor function.

解决方案

令\(L=M+1,\)那么 \(x,y\) 的范围等价于\(L\le x\le N, L\le y\le N.\)

对固定整数 \(k\ge 0\)，条件\(\left\lfloor \dfrac{y^2}{x^2}\right\rfloor = k\Longleftrightarrow k\le \dfrac{y^2}{x^2} < k+1 \Longleftrightarrow \sqrt{k}x \le y < \sqrt{k+1}x.\)

由于 \(y\) 是整数，区间整数解个数为

\[ \#\{y\in\mathbb Z:\sqrt{k}x \le y < \sqrt{k+1}x\}=\max(0,\min(N,\lceil \sqrt{k+1}x\rceil-1)-\max(L,\lceil \sqrt{k}x\rceil)+1). \]

题目只要奇数层，因此

\[ R(M,N)=\sum_{\substack{k\ge 1\\k\equiv 1\pmod 2}}\sum_{x=L}^{N}\#\{y:\sqrt{k}x \le y < \sqrt{k+1}x,L\le y\le N\}. \]

对于给定 \(k\)，要有解必须存在 \(y\le N\) 且 \(y\ge \sqrt{k}x\)，即\(\displaystyle{\sqrt{k}x\le N\Longleftrightarrow x\le \left\lfloor \frac{N}{\sqrt{k}}\right\rfloor.}\) 记\(X_0=\left\lfloor \dfrac{N}{\sqrt{k}}\right\rfloor,X_1=\left\lfloor \dfrac{N}{\sqrt{k+1}}\right\rfloor.\) 因为 \(\sqrt{k+1}>\sqrt{k}\)，所以 \(X_1\le X_0\)。

于是：

  • 当 \(x\le X_1\) 时，\(\sqrt{k+1},x\le N\)，上界不触顶；
  • 当 \(X_1 < x \le X_0\) 时，\(\sqrt{k+1},x>N\)，上界被 \(N\) 截断；
  • 当 \(x> X_0\) 时，无解。

因此只需处理 \(x\in[L,X_0]\)。

对整数 \(x\)，令\(A_k(x)=\left\lceil x\sqrt{k}\right\rceil,B_k(x)=\left\lceil x\sqrt{k+1}\right\rceil-1.\)

如果区间 \(L\le x\le X_1\)，此时 \(B_k(x)\le N\)，不需要截断，上层可取到 \(B_k(x)\)，于是该 \(x\) 贡献为\((B_k(x)-A_k(x)+1)=\left(\lceil x\sqrt{k+1}\rceil-1\right)-\lceil x\sqrt{k}\rceil+1=\lceil x\sqrt{k+1}\rceil-\lceil x\sqrt{k}\rceil.\)

对于区间 \(X_1 < x\le X_0\)，此时上界被 \(N\) 截断，贡献为\(N-\lceil x\sqrt{k}\rceil+1=(N+1)-\lceil x\sqrt{k}\rceil.\)

于是对固定奇数 \(k\)，总贡献为\(\displaystyle{S_k=\sum_{x=L}^{X_1}(\lceil x\sqrt{k+1}\rceil-\lceil x\sqrt{k}\rceil\bigr)+\sum_{x=X_1+1}^{X_0}\bigl((N+1)-\lceil x\sqrt{k}\rceil).}\)整理得

\[ S_k =\sum_{x=L}^{X_1}\lceil x\sqrt{k+1}\rceil +(X_0-X_1)(N+1) -\sum_{x=L}^{X_0}\lceil x\sqrt{k}\rceil. \]

所以 \[ R(M,N)=\sum_{\substack{1\le k\le k_{\max} \\ k\equiv 1\pmod 2}} S_k. \]

接下来计算\(k\)能得到的最大取值\(k_{\max}\)。最大可能的比值发生在最小 \(x=L\)、最大 \(y=N\)：\(\max\left\{\dfrac{y^2}{x^2}\right\}=\dfrac{N^2}{L^2}.\)因此\(k_{\max}=\left\lfloor \left(\dfrac{N}{L}\right)^2\right\rfloor\)

对非平方整数 \(d\)，\(x\sqrt d\) 从不为整数，所以\(\lceil x\sqrt d\rceil=\lfloor x\sqrt d\rfloor+1.\)于是\(\displaystyle{\sum_{x=L}^{T}\lceil x\sqrt d\rceil=\sum_{x=L}^{T}\lfloor x\sqrt d\rfloor + (T-L+1).}\)

而当 \(d\) 为完全平方数 \(d=r^2\) 时，\(\displaystyle{\sum_{x=L}^{T}\lceil x\sqrt d\rceil=\sum_{x=L}^{T}rx=r\cdot \frac{(L+T)(T-L+1)}{2}.}\)所以问题归结为高效计算\(\displaystyle{F(d,T)=\sum_{x=1}^{T}\lfloor x\sqrt d\rfloor.}\)

我们接下来基于 Beatty / Euclid思想 计算 \(F(d,T)\)。

令\(\alpha=\dfrac{a+b\sqrt d}{c}, c>0,\)考虑\(\displaystyle{S(T,\alpha)=\sum_{x=1}^{T}\lfloor x\alpha\rfloor.}\)

如果 \(\alpha\ge 1\)，设 \(q=\lfloor\alpha\rfloor\)，\(\alpha'=\alpha-q\in[0,1)\)，则 \(\lfloor x\alpha\rfloor=\lfloor x(\alpha'+q)\rfloor=\lfloor x\alpha'\rfloor+qx,\)因此

\[ S(T,\alpha)=S(T,\alpha')+q\cdot\frac{T(T+1)}{2}. \]

如果 \(0<\alpha<1\)，设\(m=\lfloor T\alpha\rfloor,\)经典恒等式给出 \[ S(T,\alpha)=Tm-\sum_{y=1}^{m}\left\lfloor \frac{y}{\alpha}\right\rfloor. \]

而\(\dfrac1\alpha=\dfrac{c}{a+b\sqrt d}=\dfrac{c(a-b\sqrt d)}{a^2-b^2d},\)仍是同类型的二次无理数，因此可以继续递归。递归深度是 \(O(\log T)\) 量级。

最终，我们进行如下计算：

  1. 设 \(L=M+1\)，计算\(k_{\max}=\left\lfloor\dfrac{N^2}{L^2}\right\rfloor.\)

  2. 对每个奇数 \(k=1,3,5,\dots,k_{\max}\)：

    • 计算\(X_0=\left\lfloor\dfrac{N}{\sqrt{k}}\right\rfloor,X_1=\left\lfloor\dfrac{N}{\sqrt{k+1}}\right\rfloor.\) 若 \(X_0<L\) 跳过。
    • 计算两次前缀的 \(\sum\lceil x\sqrt d\rceil\)（\(d=k\) 与 \(d=k+1\)），代入\(\displaystyle{S_k=\sum_{x=L}^{X_1}\lceil x\sqrt{k+1}\rceil+(X_0-X_1)(N+1)-\sum_{x=L}^{X_0}\lceil x\sqrt{k}\rceil.}\)

  3. 累加所有 \(S_k\) 得 \(R(M,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
from tools import gcd, int_sqrt
M = 2_000_000
N = 1_000_000_000

def floor_quad(a: int, b: int, c: int, d: int) -> int:
    """
    q = floor((a + b*sqrt(d))/c), with integer arithmetic.
    Works for any integers a,b and nonzero integer c (we normalize c>0).
    """
    if c < 0:
        a, b, c = -a, -b, -c

    if b >= 0:
        t = int_sqrt(b * b * d)  # floor(b*sqrt(d))
        return (a + t) // c

    bb = -b
    t = int_sqrt(bb * bb * d)   # floor(bb*sqrt(d))
    q = (a - t) // c         # may be 1 too large
    # Fix if needed (at most a few steps; in theory only 0/1 step)
    while True:
        r = a - q * c
        # need q*c <= a - bb*sqrt(d)  <=>  bb*sqrt(d) <= r
        if r > 0 and r * r >= bb * bb * d:
            return q
        q -= 1


def sum_floor_quad(n: int, a: int, b: int, c: int, d: int) -> int:
    """
    Sum_{i=1..n} floor(i*(a + b*sqrt(d))/c) using a Euclid/Beatty-style recursion,
    implemented iteratively with a stack.
    """
    if n <= 0:
        return 0

    res = 0
    stack = [(n, a, b, c, 1)]  # (n,a,b,c, sign)

    while stack:
        n, a, b, c, sgn = stack.pop()
        while n > 0:
            if c < 0:
                a, b, c = -a, -b, -c

            g = gcd(gcd(abs(a), abs(b)), c)
            if g > 1:
                a //= g
                b //= g
                c //= g

            q = floor_quad(a, b, c, d)
            if q != 0:
                res += sgn * q * n * (n + 1) // 2
                a -= q * c
                continue

            # now 0 < alpha < 1
            m = floor_quad(a * n, b * n, c, d)  # floor(n*alpha)
            res += sgn * n * m

            # 1/alpha = c/(a + b*sqrt(d)) = (c*a - c*b*sqrt(d)) / (a^2 - b^2*d)
            a2 = c * a
            b2 = -c * b
            c2 = a * a - b * b * d
            stack.append((m, a2, b2, c2, -sgn))
            break

    return res


def sum_floor_sqrt(d: int, n: int) -> int:
    """Sum_{i=1..n} floor(i*sqrt(d))."""
    if n <= 0:
        return 0
    r = int_sqrt(d)
    if r * r == d:
        return r * n * (n + 1) // 2
    return sum_floor_quad(n, 0, 1, 1, d)


def floor_div_sqrt(N: int, d: int) -> int:
    """floor(N/sqrt(d)) with integer arithmetic."""
    r = int_sqrt(d)
    if r * r == d:
        return N // r
    n2 = N * N
    x = int_sqrt(n2 // d)
    while (x + 1) * (x + 1) * d <= n2:
        x += 1
    while x * x * d > n2:
        x -= 1
    return x


def prefix_ceil_sqrt(d: int, t: int, L: int) -> int:
    """Sum_{x=L..t} ceil(x*sqrt(d))."""
    if t < L:
        return 0
    r = int_sqrt(d)
    if r * r == d:
        cnt = t - L + 1
        return r * (L + t) * cnt // 2
    cnt = t - L + 1
    # non-square: ceil(x*sqrt(d)) = floor(x*sqrt(d)) + 1
    return (sum_floor_sqrt(d, t) - sum_floor_sqrt(d, L - 1)) + cnt


def R(M: int, N: int) -> int:
    """
    Count lattice points (x,y) with M < x <= N, M < y <= N, and floor(y^2/x^2) odd.
    """
    L = M + 1
    # ✅ correct: kmax = floor((N/L)^2) = floor(N^2/L^2)
    kmax = (N * N) // (L * L)

    ans = 0
    for k in range(1, kmax + 1, 2):
        X0 = floor_div_sqrt(N, k)
        if X0 < L:
            continue
        X1 = floor_div_sqrt(N, k + 1)

        if X1 >= L:
            ans += (
                prefix_ceil_sqrt(k + 1, X1, L)
                + (X0 - X1) * (N + 1)
                - prefix_ceil_sqrt(k, X0, L)
            )
        else:
            ans += (X0 - L + 1) * (N + 1) - prefix_ceil_sqrt(k, X0, L)

    return ans


if __name__ == "__main__":
    print(R(M,N))