Project Euler 383

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

Project Euler 383

题目

Divisibility comparison between factorials

Let \(f_5(n)\) be the largest integer \(x\) for which \(5^x\) divides \(n\).

For example, \(f_5(625000) = 7\).

Let \(T_5(n)\) be the number of integers \(i\) which satisfy \(f_5((2\cdot i-1)!) < 2\cdot f_5(i!)\) and \(1 \le i \le n\).

It can be verified that \(T_5(10^3) = 68\) and \(T_5(10^9) = 2408210\).

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

解决方案

令\(N=10^{18}\)，\(v_5(n)\) 表示 \(5\) 在 \(n\) 中的指数，即最大的 \(x\) 使得 \(5^x\mid n\)。并记\(\displaystyle{v_5(n!)=\sum_{k\ge 1}\left\lfloor \frac{n}{5^k}\right\rfloor.}\)

先把 \((2i-1)!\) 换成 \((2i)!\)：\(v_5((2i-1)!)=v_5((2i)!)-v_5(2i)=v_5((2i)!)-v_5(i),\)因为 \(2\) 与 \(5\) 互素，所以 \(v_5(2i)=v_5(i)\)。

于是原不等式等价于\(v_5((2i)!)-v_5(i)<2v_5(i!)\Longleftrightarrow v_5((2i)!)-2v_5(i!)<v_5(i).\)

注意左边正是二项式系数的 \(5\)-进指数：\(v_5\dbinom{2i}{i}=v_5((2i)!)-2v_5(i!).\)因此条件等价于\(v_5\dbinom{2i}{i}<v_5(i).\)

这里需要用到一个在题目323中也出现过的结论：对素数 \(p\)，二项式系数的 \(p\)-进指数等于对应的进位次数（Kummer 定理的一种表述）。

具体到 \(p=5\)，令 \(c(i)\) 表示把 \(i\) 与 \(i\) 用 五进制相加时产生的进位次数，则有\(v_5\dbinom{2i}{i}=c(i).\)

直观上，在五进制做加法时，每产生一次进位，数位和就会比“未进位时的数位和”少 \(5-1=4\)；而 \(v_5\dbinom{2i}{i}\) 正是在度量这种“被进位消耗掉的 \(5\) 的次数”，两者因此完全一致。于是原判定条件立刻化为\(c(i) < v_5(i).\)

其中 \(v_5(i)\) 也等价于 \(i\) 的五进制表示末尾连续 \(0\) 的个数（尾随零个数）。所以若 \(5\nmid i\)，则 \(v_5(i)=0\)，不可能满足 \(c(i)<0\)，因此只有 \(5\) 的倍数才会被计入。因此我们要数的就是：五进制下把 \(i\) 翻倍的进位次数，严格小于 \(i\) 的尾随零个数。

那么设 \(t=v_5(i)\ge 1\)，写成\(i=5^t m, 5\nmid m.\)因为末尾附加 \(t\) 个 \(0\) 不会引入任何新的进位（低位全是 \(0\)，从最低位开始翻倍不会产生进位，且不会把进位“传”进来），因此\(c(i)=c(m).\)条件 \(c(i)<v_5(i)\) 等价于\(c(m) < t.\)

同时范围 \(i\le N\)等价于\(m \le \left\lfloor \dfrac{N}{5^t}\right\rfloor.\)因此我们可以按 \(t\) 求和：令 \(M_t=\left\lfloor \dfrac{N}{5^t}\right\rfloor\)，则\(\displaystyle{T_5(N)=\sum_{t\ge 1}\#\left\{1\le m\le M_t:5\nmid m,c(m)\le t-1\right\}.}\)

而 \(t\) 最大只到 \(T=\left\lfloor \log_5 N\right\rfloor\)。因此\(t\)的值可以枚举计算。

定义计数函数\(A(M,s)=\#\{0\le x\le M:c(x)\le s\}.\)那么我们要的“去掉 \(5\) 的倍数”的计数可以用一个非常省事的差分完成：

  • \(5\mid x \iff x=5y\)；
  • 且 \(c(5y)=c(y)\)（尾部补一个 \(0\) 不产生进位）。

所以\(\#\{1\le x\le M:5\nmid x,c(x)\le s\}= A(M,s)-A(\lfloor M/5\rfloor,s)\)

因为 \(x=0\) 会在两边同时出现并抵消。于是第 \(t\) 层贡献就是\(A(M_t,t-1)-A\left(\left\lfloor \dfrac{M_t}{5}\right\rfloor,t-1\right).\)剩下的任务就是快速算 \(A(M,s)\)。

因此，我们在 五进制下做“翻倍”，从低位到高位处理。设某一位的输入进位为 \(c\in{0,1}\)，该位数字为 \(d\in\{0,1,2,3,4\}\)，则这一位计算后向更高位输出的进位为\(c'=\left\lfloor\dfrac{2d+c}{5}\right\rfloor\in\{0,1\}.\)

并且“这一位是否产生进位”恰好等于 \(c'\)（因为只有 \(c'=1\) 时才发生进位）。

令 \(\ell\) 表示后缀长度（位数，允许前导 \(0\)）。定义三维 DP：\(f(\ell,c,k)\)表示所有长度为 \(\ell\) 的五进制后缀 \(x\) 满足如下条件的个数: 把 \(x\) 翻倍时，\(c\in\{0,1\}\) 表示处理完这 \(\ell\) 位之后，向更高一位输出的进位；\(k\) 表示在这 \(\ell\) 位内部一共发生了多少次进位。

长度为 \(0\) 的空后缀没有任何进位，因此有\(f(0,0,0)=1, f(0,1,k)=0, f(0,0,k)=0(k\ne 0).\)

从长度 \(\ell-1\) 扩展到 \(\ell\)：在最低位新添一个数字 \(d\in\{0,1,2,3,4\}\)。设新添最低位的输入进位为 \(c_{\text{in}}\)，输出进位为\(\displaystyle{c_{\text{out}}=\left\lfloor \frac{2d+c_{\text{in}}}{5}\right\rfloor.}\)

由于“这一位产生的进位次数”就是 \(c_{\text{out}}\)，因此总进位次数也随之增加 \(c_{\text{out}}\)。于是有转移：

\[ f(\ell,c_{\text{out}},k+c_{\text{out}})\leftarrow f(\ell-1,c_{\text{in}},k),\forall d\in\{0,1,2,3,4\},c_{\text{in}}\in\{0,1\}. \]

为了后续能在 \(O(1)\) 得到“进位次数不超过 \(K\) 的数量”，我们对每个 \((\ell,c)\) 做前缀和：\(\displaystyle{F(\ell,c,K)=\sum_{k\le K} f(\ell,c,k).}\)

我们需要根据DP结果计算\(A(M,s)\)。把 \(M\) 写成五进制数位\(M=(d_0d_1\dots d_{L-1})_5,\)其中 \(d_0\) 为最高位，\(d_{L-1}\) 为最低位。

按“首个变小位”的典型写法：枚举第一个位置 \(i\)，在该位选一个 \(e<d_i\)，那么：

  • 左侧高位前缀 \((d_0\dots d_{i-1}e)\) 被固定；
  • 右侧低位一共还有 \(\ell=L-i-1\) 位，且因为已经变小了，这 \(\ell\) 位可以任意取 \(0\sim 4\)，也就是一个“自由后缀”。

但进位是从低位向高位传播的，所以自由后缀不仅贡献自己的进位次数，还会向左侧前缀“吐出”一个进位 \(c\in\{0,1\}\)，作为前缀最低位的输入进位。

这正好与针对\(f\)的后缀 DP 对接：

  • 用 \(F(\ell,c,K)\) 统计“长度为 \(\ell\) 的自由后缀，最终向左输出进位为 \(c\)，且内部进位次数 \(\le K\)”的数量；
  • 对每个可能的 \(c\)，我们在固定前缀上直接模拟一遍翻倍的进位次数。

固定了某个位置 \(i\) 和该位取值 \(e<d_i\) 后：

  • 前缀为 \(P=(d_0\dots d_{i-1}e)\)，长度为 \(i+1\)；
  • 后缀长度为 \(\ell=L-i-1\)；
  • 枚举后缀向左输出的进位 \(c\in\{0,1\}\)。

令 \(c_P(P,c)\) 表示“把前缀 \(P\) 翻倍时，若其最低位输入进位为 \(c\)，前缀内部产生的进位次数”。那么为了让整体进位次数 \(\le s\)，后缀允许的进位次数上界为\(K = s - c_P(P,c).\)只要 \(K\ge 0\)，该组 \((i,e,c)\) 的贡献就是\(F(\ell,c,K).\)

把所有情况相加，并且最后还要单独判断 \(x=M\) 本身是否满足 \(c(M)\le s\)，即可得到

\[ A(M,s)=\left(\sum_{i=0}^{L-1}\sum_{e=0}^{d_i-1}\sum_{c\in\{0,1\}}\mathbf 1\{K\ge 0\}\cdot F(L-i-1,c,K)\right)+\mathbf 1\{c(M)\le s\}. \]

代码

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

N = 10**18

def base5_digits(n: int) -> List[int]:
    """Return base-5 digits of n in MSB -> LSB order."""
    if n == 0:
        return [0]
    ds = []
    while n:
        n, r = divmod(n, 5)
        ds.append(r)
    return ds[::-1]

def precompute_suffix_cum(Lmax: int):
    """
    suf_cum[len][carry_out][k] = sum_{j<=k} count of len-digit base-5 suffixes
    (allow leading zeros), whose doubling from LSB with initial carry=0
    results in:
      - carry_out to the next higher digit = carry_out (0/1)
      - total number of carries inside these len digits = j
    """
    suf = []
    suf_cum = []

    # len = 0: empty suffix -> carry_out=0, carries=0
    suf.append([[1], [0]])
    suf_cum.append([[1], [0]])

    for length in range(1, Lmax + 1):
        prev0, prev1 = suf[length - 1]  # arrays of size length
        dp0 = [0] * (length + 1)
        dp1 = [0] * (length + 1)

        for carry_in, prev in ((0, prev0), (1, prev1)):
            for k, val in enumerate(prev):
                if val == 0:
                    continue
                for d in range(5):
                    s = 2 * d + carry_in
                    carry_out = 1 if s >= 5 else 0
                    k2 = k + carry_out
                    if carry_out == 0:
                        dp0[k2] += val
                    else:
                        dp1[k2] += val

        suf.append([dp0, dp1])

        cum0 = [0] * (length + 1)
        cum1 = [0] * (length + 1)
        run = 0
        for k in range(length + 1):
            run += dp0[k]
            cum0[k] = run
        run = 0
        for k in range(length + 1):
            run += dp1[k]
            cum1[k] = run
        suf_cum.append([cum0, cum1])

    return suf_cum

def carry_count_digits(digits_msb: List[int], carry_in_lsb: int = 0) -> int:
    """
    Count number of carries when doubling the number represented by digits_msb (MSB->LSB),
    assuming the carry into the least significant digit equals carry_in_lsb (0/1).
    """
    carry = carry_in_lsb
    cnt = 0
    for d in reversed(digits_msb):  # LSB -> MSB
        s = 2 * d + carry
        carry = 1 if s >= 5 else 0
        cnt += carry
    return cnt

def A_leq(M: int, s: int, suf_cum) -> int:
    """
    A(M,s) = #{ 0 <= x <= M : carry_count(2*x in base 5) <= s }.
    """
    if M < 0:
        return 0
    digits = base5_digits(M)
    L = len(digits)

    ans = 0

    for idx in range(L):
        rem = L - idx - 1
        bd = digits[idx]
        prefix_base = digits[:idx]

        # choose digit smaller than bound at this position
        for d in range(bd):
            prefix_digits = prefix_base + [d]
            # carry_in to this prefix (its LSB) is the carry_out from the free suffix
            for cin in (0, 1):
                cp = carry_count_digits(prefix_digits, cin)
                allowed = s - cp
                if allowed < 0:
                    continue
                if allowed > rem:
                    allowed = rem
                ans += suf_cum[rem][cin][allowed]

    # include M itself if it satisfies
    if carry_count_digits(digits, 0) <= s:
        ans += 1

    return ans

def T5(N: int) -> int:
    """
    Compute T_5(N) = #{ 1<=i<=N : v5((2i-1)!) < 2*v5(i!) }.
    Using equivalence: carry_count_base5(i+i) < v5(i).
    """
    Lmax = len(base5_digits(N))
    suf_cum = precompute_suffix_cum(Lmax)

    res = 0
    p5 = 5
    t = 1
    while p5 <= N:
        M = N // p5
        s = t - 1
        res += A_leq(M, s, suf_cum) - A_leq(M // 5, s, suf_cum)
        t += 1
        p5 *= 5
    return res


print(T5(N))