Project Euler 838

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

Project Euler 838

题目

Not Coprime

Let \(f(N)\) be the smallest positive integer that is not coprime to any positive integer \(n \le N\) whose least significant digit is \(3\).

For example \(f(40)\) equals to \(897 = 3 \cdot 13 \cdot 23\) since it is not coprime to any of \(3,13,23,33\). By taking the natural logarithm (log to base \(e\)) we obtain \(\ln f(40) = \ln 897 \approx 6.799056\) when rounded to six digits after the decimal point.

You are also given \(\ln f(2800) \approx 715.019337\).

Find \(f(10^6)\). Enter its natural logarithm rounded to six digits after the decimal point.

解决方案

令\(M=10^6,A_N=\{n\in\mathbb Z_{>0}\mid n\le N,n\equiv 3\pmod{10}\}.\)

若 \(m\) 满足 \(\gcd(m,n)>1\)，则对每个 \(n\) 必存在某个素数 \(p\) 使得 \(p\mid m\) 且 \(p\mid n\)。

因此只需要关心 \(m\) 的不同素因子集合：若 \(m\) 含有 \(p^e\)，把它替换为 \(p\) 不会改变与任何 \(n\) 的最大公因子是否大于 \(1\)。所以最优解可取为平方自由形式\(\displaystyle{m=\prod_{p\in S}p}\)，其中 \(S\) 为某个素数集合。将目标“最小化 \(m\)”改写为“最小化 \(\ln m\)”，因为 \(\ln\) 单调：\(\displaystyle{\ln m=\sum_{p\in S}\ln p.}\)

对每个 \(n\in A_N\)，设其素因子集合为\(P(n)=\{p\in\Omega: p\mid n\}.\)其中\(\Omega\)是素数集合。条件 \(\gcd(m,n)>1\) 等价于\(P(n)\cap S\ne\varnothing.\)于是问题变成一个带权集合覆盖：选取素数集合 \(S\)，使得对所有 \(n\in A_N\) 都被“命中”，并最小化 \(\displaystyle{\sum_{p\in S}\ln p}\)。

接下来有两类素数必选。

  1. 所有 \(p\equiv 3\pmod{10}\) 且 \(p\le N\) 必须选。 若 \(p\) 是素数且 \(p\equiv 3\pmod{10}\)，则 \(p\in A_N\)（当 \(p\le N\)）。要满足 \(\gcd(m,p)>1\)，只能是 \(p\mid m\)，因此\(p\in S.\) 所以所有这类素数必须选入。
  2. 若 \(p\equiv 7\pmod{10}\) 且 \(p^3\le N\)，则 \(p\) 也必须选。 因为\(p\equiv 7\pmod{10}\Rightarrow p^3\equiv 7^3\equiv 3\pmod{10},\)所以当 \(p^3\le N\) 时，\(p^3\in A_N\)，且 \(P(p^3)={p}\)。为了命中 \(p^3\)，必须有 \(p\in S\)。

因此强制选入的部分是：所有 \(p\le N\) 且 \(p\equiv 3\pmod{10}\)；所有 \(p^3\le N\) 且 \(p\equiv 7\pmod{10}\)。

把这些素数的对数和记为 \(L_1\)，并将所有含它们因子的 \(n\in A_N\) 视为已覆盖，剩下只需处理“未覆盖”的 \(n\)。

接下来只看仍未被覆盖的 \(n\in A_N\)。此时 \(n\) 不能含有上面强制选入的素因子。注意 \(n\equiv 3\pmod{10}\) 一定是奇数且不被 \(5\) 整除，因此其所有素因子都在 \({1,3,7,9}\pmod{10}\) 中。

接下来证明，含有 \(p\equiv 1\pmod{10}\) 的数是冗余的。若 \(n\in A_N\) 且存在素数 \(p\equiv 1\pmod{10}\) 使 \(p\mid n\)，则令 \(n_1=n/p\)，有\(n_1\le n\le N, n_1\equiv 3\pmod{10}.\)并且\(P(n_1)\subseteq P(n).\)因此覆盖了 \(n_1\) 必然覆盖 \(n\)，因此这类 \(n\) 的约束可忽略。于是可只保留不含 \(1\pmod{10}\) 素因子的约束。

剔除 \(1\) 类后，且 \(3\) 类已强制全选，未覆盖的 \(n\) 的素因子只能来自 \(\{7,9\}\pmod{10}\)。

要使乘积末位为 \(3\)，最小形式是\(7\cdot 9\equiv 3\pmod{10}.\)

若 \(n\) 含有多于两个素因子（例如还多乘了别的 \(7\) 或 \(9\)），则其某个因子子集（例如取其中一个 \(7\) 与一个 \(9\)）的乘积仍为 \(3\pmod{10}\)，且不超过 \(n\le N\)，从而给出一个更小的约束集合，导致原约束被“子集支配”而可忽略。

再结合 \(p\equiv 7\pmod{10}\) 且 \(p^3\le N\) 已被强制选入（因此任何含这类 \(p\) 的数都已覆盖），最终留下的本质约束只来自形如

\[ n=p_7\cdot p_9\le N, p_7\equiv 7\pmod{10},\quad p_9\equiv 9\pmod{10}, \]

并且 \(p_7^3>N\)（否则已被强制覆盖）。

为解决这个问题，我们构造二分图。令

  • 左侧点集 \(U\)：所有满足 \(p\le N,p\equiv 7\pmod{10},p^3>N\) 的素数 \(p\)；
  • 右侧点集 \(V\)：所有满足 \(q\le N,q\equiv 9\pmod{10}\) 的素数 \(q\)；
  • 若 \(p\in U,q\in V\) 且 \(pq\le N\)，连一条边。

每条边对应一个需要命中的数 \(pq\)，要求选点集合 \(C\subseteq U\cup V\) 使得每条边至少有一个端点在 \(C\) 中；代价为\(\displaystyle{\sum_{p\in C\cap U}\ln p+\sum_{q\in C\cap V}\ln q.}\) 这就是一个带权二分图最小点覆盖，且权重是 \(\ln x\)。

关键是该图具有单调结构：把 \(U,V\) 按素数从小到大排序。对固定 \(p\)，其邻居是所有满足\(q\le \left\lfloor\dfrac{N}{p}\right\rfloor\)的 \(q\)，因此邻居集合是 \(V\) 的一个前缀；且 \(p\) 越大，前缀越短。由于 \(\ln\) 随素数增大而增大，可以把任意解压缩成前缀形式：存在最优解满足

  • 选取 \(U\) 的前 \(i\) 个素数；
  • 选取 \(V\) 的前 \(j\) 个素数。

若选了 \(U\) 的前 \(i\) 个，那么未选的最小左点是 \(p=U[i]\)（若 \(i=|U|\) 则无未选点）。为了覆盖所有与未选左点相连的边，必须选取所有\(q\le \left\lfloor\dfrac{N}{U[i]}\right\rfloor\)的右点，也就是某个最小的前缀长度 \(j(i)\)。于是总代价为\(\displaystyle{c(i)=\sum_{t< i}\ln U[t]+\sum_{t< j(i)}\ln V[t],}\)对所有 \(i\) 取最小即可。\(j(i)\) 随 \(i\) 单调不增，可以双指针线性扫描完成。

最终有

\[ \ln f(N)=L_1+\min_i\{c(i)\}. \]

代码

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
from math import log

from tools import Factorizer

N = 10 ** 6
def ln_f(N: int) -> float:
    ps = Factorizer(N).primes

    p3 = []
    p7_small = []
    p7_large = []
    p9 = []

    for p in ps:
        r = p % 10
        if r == 3:
            p3.append(p)
        elif r == 7:
            if p * p * p <= N:
                p7_small.append(p)
            else:
                p7_large.append(p)
        elif r == 9:
            p9.append(p)

    base = 0.0
    for p in p3:
        base += log(p)
    for p in p7_small:
        base += log(p)

    A = p7_large
    B = p9

    prefA = [0.0]
    s = 0.0
    for p in A:
        s += log(p)
        prefA.append(s)

    prefB = [0.0]
    s = 0.0
    for p in B:
        s += log(p)
        prefB.append(s)

    best = float('inf')
    j = len(B)
    for i in range(len(A) + 1):
        if i == len(A):
            j = 0
        else:
            thr = N // A[i]
            while j > 0 and B[j - 1] > thr:
                j -= 1
        cost = prefA[i] + prefB[j]
        if cost < best:
            best = cost

    return base + best

ans = ln_f(N)
print(f"{ans:.6f}")