Project Euler 466

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

Project Euler 466

题目

Distinct terms in a multiplication table

Let \(P(m,n)\) be the number of distinct terms in an \(m\times n\) multiplication table.

For example, a \(3\times4\) multiplication table looks like this:

\(\times\) \(1\) \(2\) \(3\) \(4\)
\(1\) \(1\) \(2\) \(3\) \(4\)
\(2\) \(2\) \(4\) \(6\) \(8\)
\(3\) \(3\) \(6\) \(9\) \(12\)

There are 8 distinct terms \(\{1,2,3,4,6,8,9,12\}\), therefore \(P(3,4) = 8\).

You are given that:

\(P(64,64) = 1263,\)

\(P(12,345) = 1998,\) and

\(P(32,10^{15}) = 13826382602124302\).

Find \(P(64,10^{16})\).

解决方案

把乘法表按第一维 \(a\) 分成 \(m\) 行，第 \(a\) 行的所有乘积集合为\(R_a=\{a\cdot b:1\le b\le n\}\)，于是\(\displaystyle{P(m,n)=\left|\bigcup_{a=1}^{m}R_a\right|}\)，直接求并集困难，但注意到 \(m=64\) 很小，可以把每个乘积分配给一个唯一的行来计数。

固定 \(x<y\)，考虑一个数 \(x\cdot i\) 是否也能写成 \(y\cdot j\)。令\(g=\gcd(x,y), x=gx', y=gy', \gcd(x',y')=1\)。若\(x i=y j\)代入得到\(gx'i=gy'j\Longrightarrow x'i=y'j\)。因为 \(\gcd(x',y')=1\)，所以 \(y'\mid i\)，即\(y\cdot j=x\cdot i\Longleftrightarrow\dfrac{y}{\gcd(x,y)}\mid i\)

因此，对固定的 \(x\)，只要存在某个 \(y\in\{x+1,\dots,m\}\) 使得\(d_{x,y}=\dfrac{y}{\gcd(x,y)}\)整除 \(i\)，那么 \(x i\) 就会在第 \(y\) 行出现，从而不是“独占”的乘积。

对固定的 \(x\)，定义禁止整除集合\(D_x=\left\{\dfrac{y}{\gcd(x,y)}:x<y\le m\right\}\)，那么第 \(x\) 行中那些不会在更大行号重复的项对应的 \(i\) 是：\(1\le i\le n, \forall d\in D_x,d\nmid i\)。

记这类 \(i\) 的个数为 \(U_x(n)\)，则

\[P(m,n)=\sum_{x=1}^{m}U_x(n)\]

原因是：任意乘积 \(t\) 在 \({1,\dots,m}\) 行中出现的所有行号集合有最大值 \(x_{\max}\)，它会且只会被 \(x_{\max}\) 这一行计入一次，从而恰好统计出不同乘积总数。

仙子我们需要计算\(U_x(n)=\#\{1\le i\le n:\forall d\in D_x,d\nmid i\}\)，等价于\(U_x(n)=n-\#\{1\le i\le n:\exists d\in D_x,d\mid i\}\)，即典型的“被若干个数整除的并集计数”，可用容斥思想，但 \(D_x\) 中元素并不互素，需要更一般的做法。

定义一个更一般的函数：\(C(t,\mathcal S)=\#\{1\le i\le n:t\mid i,\forall s\in\mathcal S,s\nmid i\}\)，那么我们最终需要的是\(U_x(n)=C(1,D_x)\)。取 \(\mathcal S\) 的第一个元素为 \(k\)，把满足 \(t\mid i\) 的数分成两类：

  • 不被 \(k\) 整除的：贡献 \(C(t,\mathcal S\setminus\{k\})\)
  • 被 \(k\) 整除的：等价于被 \(\mathrm{lcm}(t,k)\) 整除，贡献 \(C(\mathrm{lcm}(t,k),\mathcal S\setminus\{k\})\)

因此得到递推：

  • 若 \(t\equiv 0\pmod k\)，则所有 \(t\) 的倍数都自动是 \(k\) 的倍数，会被排除，所以\(C(t,\mathcal S)=0\)
  • 否则\(C(t,\mathcal S)=C(t,\mathcal S\setminus\{k\})-C(\mathrm{lcm}(t,k),\mathcal S\setminus\{k\})\)

可见，递归边界为\(C(t,\varnothing)=\left\lfloor \dfrac{n}{t}\right\rfloor\)。这样就完成了“对任意整除集合”的计数。

若 \(a,b\in D_x\) 且 \(a\mid b\)，那么\(\{i:b\mid i\}\subseteq \{i:a\mid i\}\)。因此在“禁止被任意一个整除”的问题里，\(b\) 永远不会增加额外限制，可以删除。最终保留一个在整除意义下的“极小集合”，递归树会显著变小，速度非常快。

代码

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
from functools import lru_cache
from typing import List

from tools import gcd
M = 64
N = 10 ** 16

def simplify_divisors(divs: List[int]) -> List[int]:
    divs = sorted(set(divs), reverse=True)
    res = []
    for d in divs:
        bad = False
        for r in res:
            if d % r == 0:
                bad = True
                break
        if bad:
            continue
        i = 0
        while i < len(res):
            if res[i] % d == 0:
                res.pop(i)
            else:
                i += 1
        res.append(d)
    res.sort(reverse=True)
    return res


def count_not_divisible(n: int, divs: List[int]) -> int:
    divs = tuple(divs)

    @lru_cache(maxsize=None)
    def rec(idx: int, cur: int) -> int:
        if cur > n:
            return 0
        if idx == len(divs):
            return n // cur
        k = divs[idx]
        if cur % k == 0:
            return 0
        nxt = cur // gcd(cur, k) * k
        return rec(idx + 1, cur) - rec(idx + 1, nxt)

    return rec(0, 1)


def P(m: int, n: int) -> int:
    total = 0
    for x in range(m, 0, -1):
        divs = [y // gcd(x, y) for y in range(x + 1, m + 1)]
        divs = simplify_divisors(divs)
        total += count_not_divisible(n, divs)
    return total

print(P(M,N))