Project Euler 604

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

Project Euler 604

题目

Convex path in square

Let \(F(N)\) be the maximum number of lattice points in an axis-aligned \(N\times N\) square that the graph of a single strictly convex increasing function can pass through.

You are given that \(F(1) = 2, F(3) = 3, F(9) = 6, F(11) = 7, F(100) = 30\) and \(F(50000) = 1898\).

Below is the graph of a function reaching the maximum \(3\) for \(N=3\):

Find \(F(10^{18})\).

解决方案

在轴对齐的 \(N\times N\) 正方形里，取所有曲线经过的整点按 \(x\) 坐标递增排列：\((x_1,y_1),(x_2,y_2),\dots,(x_m,y_m), x_1<\cdots<x_m.\)因为函数严格递增，所以\(y_1<y_2<\cdots<y_m.\)

令相邻两点差分向量为\((\Delta x_i,\Delta y_i)=(x_{i+1}-x_i,y_{i+1}-y_i), i=1,\dots,m-1.\)则\(\Delta x_i>0, \Delta y_i>0.\)

严格凸的等价离散条件是：相邻割线斜率严格递增\(\dfrac{\Delta y_1}{\Delta x_1}<\dfrac{\Delta y_2}{\Delta x_2}<\cdots<\dfrac{\Delta y_{m-1}}{\Delta x_{m-1}}.\)因此，整条曲线在这些整点之间的骨架就是一条斜率严格递增的格点折线（凸折线）。

同时，因为所有点都在 \(N\times N\) 正方形中，整体横向、纵向位移都不超过 \(N\)，于是\(\displaystyle\sum_{i=1}^{m-1}\Delta x_i\le N,\sum_{i=1}^{m-1}\Delta y_i\le N.\)

如果某一段 \((\Delta x,\Delta y)\) 不是互素的，令 \(g=\gcd(\Delta x,\Delta y)>1\)，则该线段内部存在整点（把向量等分）。这样就会出现三点共线，从而相邻两段割线斜率相等，违背严格凸。所以在最优解中可以假设每一段都满足\(\gcd(\Delta x_i,\Delta y_i)=1.\)

到这里，问题已经变成纯组合优化：选取尽可能多的正整数对 \((a_i,b_i)\)，满足\(\gcd(a_i,b_i)=1\)；斜率 \(b_i/a_i\) 两两不同；\(\sum a_i\le N,\sum b_i\le N\)。设最多能选到 \(k\) 对，则\(F(N)=k+1.\)最后的 \(+1\) 是因为 \(k\) 段折线连接 \(k+1\) 个整点。

最终，上面的条件里唯一看似几何的是凸性，但它已经被斜率互异吸收掉了：只要选到一组互素向量且斜率互不相同，就把它们按斜率从小到大排序；从 \((0,0)\) 开始累加这些向量得到折线顶点；由于斜率严格递增，这条折线就是严格凸的。因此，只需要最大化可选向量数，排序会自动给出严格凸折线。

接下来，我们尝试按 \(s=a+b\) 分层。对每个向量 \((a,b)\)，引入层号\(s=a+b.\)固定 \(s\)，所有正整数解为 \((a,s-a)\)，其中 \(1\le a\le s-1\)。互素条件为\(\gcd(a,s-a)=\gcd(a,s)=1.\)因此该层可用向量数量为\(|V_s|=\varphi(s).\)

更关键的是，如果把这一层全部取完，那么横向消耗与纵向消耗完全相同。原因是互素的 \(a\) 在映射 \(a\mapsto s-a\) 下成对出现，于是\(\displaystyle\sum_{(a,b)\in V_s} a=\sum_{(a,b)\in V_s} b=\frac{s\varphi(s)}{2}.\)所以：取满一层 \(s\) 会带来段数增加 \(\varphi(s)\)；两轴预算各减少 \(\dfrac{s\varphi(s)}{2}\)。

可见，每条段都要花费 \((a,b)\) 的预算，而 \(a+b=s\) 越小，单位段数的资源消耗越低，因此应优先选小 \(s\)。并且整层取满还有一个重要性质：它让\(\sum a=\sum b\)始终保持完全平衡，不会出现某一轴先耗尽导致另一轴剩余被浪费。因此最优解具有这样的形态：对 \(s=2,3,\dots,S\) 全部取满；对 \(s=S+1\) 取一部分；对更大的 \(s\) 不再取。

记\(\displaystyle L(S)=\sum_{s=2}^{S}\frac{s\varphi(s)}{2},C(S)=\sum_{s=2}^{S}\varphi(s).\)取满到 \(S\) 时已经获得 \(C(S)\) 条段，并在每个方向消耗了 \(L(S)\)。取最大的 \(S\) 使得 \(L(S)\le N\)。剩余预算为\(w=N-L(S),\)下一层为 \(s_0=S+1\)。

在层 \(s_0\) 中，若 \((a,s_0-a)\) 可用，则 \((s_0-a,a)\) 也必可用（互素对称）。把这两个向量一起拿，会在两轴各消耗 \(s_0\)，并得到 2 条段。

因此一定可以拿到\(2\left\lfloor \dfrac{w}{s_0}\right\rfloor\)条段。令\(w=q s_0+r, 0\le r<s_0,\)则这部分贡献 \(2q\) 条段，并留下每轴 \(r\) 的预算。现在问：能否再多拿 1 条 \((a,s_0-a)\)？它需要\(a\le r, s_0-a\le r\)即\(s_0-r\le a\le r.\) 只要这个区间里存在某个 \(a\) 与 \(s_0\) 互素，就可以再加 1 条段。于是最终段数为

\[k=C(S)+2q+\mathbf{1}\{\exists a\in[s_0-r,r]\text{s.t.}\gcd(a,s_0)=1\}.\]

答案为\(F(N)=k+1.\)

可见，\(s\varphi(s)\)前缀和的渐进增长为：\(\displaystyle\sum_{s\le x} s\varphi(s)\sim \frac{2}{\pi^2}x^3,\)可得\(\displaystyle L(S)=\frac12\sum_{s\le S}s\varphi(s)\sim \frac{1}{\pi^2}S^3.\)令 \(L(S)\approx N\)，得到上界为\(S\approx (N\pi^2)^{1/3}.\)

代码

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
from tools import int_cubert, gcd
from math import pi

N = 10 ** 18
def phi_sieve(n):
    phi = list(range(n + 1))
    primes = []
    is_comp = bytearray(n + 1)
    phi[0] = 0
    if n >= 1:
        phi[1] = 1
    for i in range(2, n + 1):
        if not is_comp[i]:
            primes.append(i)
            phi[i] = i - 1
        for p in primes:
            v = i * p
            if v > n:
                break
            is_comp[v] = 1
            if i % p == 0:
                phi[v] = phi[i] * p
                break
            phi[v] = phi[i] * (p - 1)
    return phi

def F(N):
    limit = int_cubert(int(N * pi ** 2)) + 20
    phi = phi_sieve(limit)
    used = 0
    seg = 0
    s = 2
    while s <= limit:
        layer = s * phi[s] // 2
        if used + layer <= N:
            used += layer
            seg += phi[s]
            s += 1
        else:
            break
    res = N - used
    q, rem = divmod(res, s)
    seg += 2 * q
    if rem * 2 < s:
        return seg + 1
    l = s - rem
    r = rem
    mid = s // 2
    ok = 0
    for d in range(0, r - l + 1):
        a = mid - d
        if a >= l and gcd(a, s) == 1:
            ok = 1
            break
        a = mid + d
        if a <= r and gcd(a, s) == 1:
            ok = 1
            break
    return seg + 1 + ok

print(F(N))