Project Euler 322

Project Euler 322

题目

Binomial coefficients divisible by \(10\)

Let \(T(m, n)\) be the number of the binomial coefficients \(^iC_n\) that are divisible by \(10\) for \(n \le i < m\)(\(i, m\) and \(n\) are positive integers).

You are given that \(T(10^9, 10^7-10) = 989697000\).

Find \(T(10^{18}, 10^{12}-10)\).

解决方案

题目定义\(T(m,n)\):在所有二项式系数 \(\dbinom{i}{n}\) 中,满足 \(n\le i<m\)能被 \(10\) 整除 的个数。

\(N = 10^{12}-10, M= 10^{18}, L = M-N\)

\(i\) 写成\(i = N + k, 0\le k < L.\)于是\(T(M,N)=\#\left\{0\le k<L : 10\mid \dbinom{N+k}{N}\right\}.\)

核心工具是 \(p\)-进制进位与 \(p\)-进制指数 \(\nu_p\) 的关系。

对素数 \(p\),勒让德公式是\(\displaystyle{\nu_p(n!)=\sum_{t\ge 1}\left\lfloor\frac{n}{p^t}\right\rfloor.}\)

\(n\) 写成 \(p\) 进制:\(\displaystyle{n=\sum_{t\ge 0} d_t p^t, 0\le d_t<p}.\)记数位和\(\displaystyle{s_p(n)=\sum_{t\ge 0} d_t.}\)此外,在题目320,我们证明了恒等式:\(\nu_p(n!)=\dfrac{n-s_p(n)}{p-1}.\)

\(i=N+k\),有\(\displaystyle{\nu_p\binom{N+k}{N}= \nu_p((N+k)!)-\nu_p(N!)-\nu_p(k!)=\frac{(N+k)-s_p(N+k)-\big(N-s_p(N)\big)-\big(k-s_p(k)\big)}{p-1}.}\)

化简得

\[\nu_p\binom{N+k}{N}=\frac{s_p(N)+s_p(k)-s_p(N+k)}{p-1}.\]

现在看 \(p\) 进制加法 \(N+k\):每发生一次进位,会使“数位和”减少恰好 \(p-1\)(因此进位次数与数位和差成正比)。因此\(s_p(N)+s_p(k)-s_p(N+k)=(p-1)\cdot C.\)

代回上式得到结论:\(\nu_p\dbinom{N+k}{N}\) 等于\(N\)\(k\)\(p\) 进制相加时的进位次数\(C\)。因此\(p\nmid\dbinom{N+k}{N}\),当且仅当 \(p\) 进制加法\(N+k\)无进位。

因为 \(10=2\cdot 5\),所以\(10\mid \dbinom{N+k}{N}\iff 2\mid \dbinom{N+k}{N}\land 5\mid \dbinom{N+k}{N}.\)

定义集合

\[A_p= \left\{0\le k<L : p\nmid \binom{N+k}{N}\right\}\]

\(A_p\) 是“对素数 \(p\) 无进位”的 \(k\) 的集合。

那么目标集合(能被 \(10\) 整除)就是\(\{0,\dots,L-1\}\setminus(A_2\cup A_5)\)因此

\[T(M,N)=L-|A_2|-|A_5|+|A_2\cap A_5|\]

所以我们需要计算三件事:

  1. \(|A_2|\):二进制无进位(\(k\& N=0\),其中\(\&\)是与运算)。
  2. \(|A_5|\):五进制无进位(逐位限制)
  3. \(|A_2\cap A_5|\):同时满足二进制与五进制无进位

其中,完成前两件事则可以使用数位动态规划解决。

\(N\) 写成 \(p\) 进制\(\displaystyle{N=\sum_{t\ge 0} n_t p^t, 0\le n_t<p.}\)

“无进位”要求每一位\(k_t + n_t < p\Longleftrightarrow 0\le k_t \le (p-1-n_t).\)于是第 \(t\) 位可选数是\(b_t = p-n_t.\)

因此 \(|A_p|\) 等价于:统计 \(0\le k<L\)\(p\) 进制表示,逐位满足 \(k_t\in[0,b_t-1]\) 的个数。

使用数位动态规划的具体解决方法是:当扫描到某一位时,若把该位取成比 \(L\) 的该位更小的合法值,则低位自由填充,贡献是“合法低位组合数”的乘积;若该位取到与 \(L\) 相等,则继续往下;若该位已不合法则停止。

一般情况下同时在 \(2\) 进制和 \(5\) 进制无进位很难直接动态规划,因为两种进制耦合,但本题的 \(N=10^{12}-10\) 有非常强的结构:

因为\(10^{12}=2^{12}\cdot 5^{12}\),取 \(12\) 位是自然的尺度。

无进位对 \(p=2\) 等价于\(k\& N=0\)

只看低 \(12\) 位:\(k\bmod 2^{12}\) 必须在 \(N\bmod 2^{12}\) 的所有 \(1\) 位上为 \(0\)。因为\(N\)的低\(12\)\((111111110110)_2\)。本题恰好只留下 \(4\) 个余数(也就是低 \(12\) 位只有两个自由位)。

\(N\) 写成 \(5\) 进制,低 \(12\) 位是\((444444444430)_5\)

无进位要求每位 \(k_t\le 4-n_t\)

  • \(n_t=4\) 时,只能 \(k_t=0\)
  • \(n_t=3\) 时,可 \(k_t\in\{0,1\}\)
  • \(n_t=0\) 时,可 \(k_t\in\{0,1,2,3,4\}\)

因此低 \(12\) 位一共只有\(2\cdot 5=10\)种组合。

因为\(10^{12}=2^{12}\cdot 5^{12}, \gcd(2^{12},5^{12})=1,\)用中国剩余定理把\(k\bmod 2^{12}\)\(k\bmod 5^{12}\)合并成\(k\bmod 10^{12}.\) 于是 \(A_2\cap A_5\) 中的所有 \(k\) 必须属于这 40 个余数类中的某一个:\(k = r + t\cdot 10^{12}, t\ge 0.\)

再对每个反向构造出来的候选 \(k\) 做完整检查:二进制无进位:\(k\& N=0\);五进制无进位:逐位检查 \(k\)\(N\)\(5\) 进制加法是否进位。由于步长达到\(10^{12}\),里面只有\(40\)个数才满足枚举要求,因此实际交集非常小,所以这一步跑得很快。

为了让五进制检查更快,我们把 \(k \bmod 5^{18}\) 拆成两块 \(5^9\)\(k \equiv k_0 + 5^9 k_1 \pmod{5^{18}},\),由于无进位在两块上是独立的(不会跨块产生进位),因此可以预先对每一块做一个大小为\(5^9=1953125\)的布尔表,检查时只需两次数组访问。

代码

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
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
from tools import CRT

# =========================
# Project Euler 322
# Count T(M, N): number of binomial coefficients C(i, N) divisible by 10
# for N <= i < M.
#
# Key facts (Lucas / Kummer viewpoint):
# - C(N+k, N) is NOT divisible by prime p <=> (N + k) has NO carry with N in base p
# equivalently: for every digit position, digit(k) <= (p-1 - digit(N)).
#
# Let L = M - N and consider k in [0, L).
# Define:
# A2 = #{k in [0, L): C(N+k, N) not divisible by 2 } (no-carry base 2)
# A5 = #{k in [0, L): C(N+k, N) not divisible by 5 } (no-carry base 5)
# A25 = #{k in [0, L): C(N+k, N) not divisible by 2 and not divisible by 5 }
#
# Then by inclusion–exclusion:
# #{k where divisible by 10} = L - A2 - A5 + A25
# and that equals T(M, N).
# =========================

M = 10 ** 18
N = 10 ** 12 - 10
L = M - N


def digits_base(n: int, p: int, length=None):
"""
Return digits of n in base p, little-endian: [d0, d1, ...] where
n = d0 + d1*p + d2*p^2 + ...
If length is given, pad with zeros or truncate to exactly `length` digits.
"""
ds = []
while n:
n, r = divmod(n, p)
ds.append(r)
if not ds:
ds = [0]
if length is not None:
if len(ds) < length:
ds += [0] * (length - len(ds))
else:
ds = ds[:length]
return ds # little-endian


def count_no_carry(n: int, B: int, p: int) -> int:
"""
Count:
|{ 0 <= k < B : adding (n + k) in base p produces NO carry with n }|
i.e. for each digit position i:
k_i <= (p - 1 - n_i)

This computes the count for the range [0, B) using digit-DP in base p.

Implementation idea:
- Let allow[i] = number of permitted digits at position i for k
= (p - n_i), since k_i in [0 .. p-1-n_i].
- For numbers < B, scan digits of B from high to low, and sum choices
where at the first differing position we pick a smaller digit (and valid),
while lower positions can be anything valid.
"""
if B <= 0:
return 0

# base-p digits for n and B
nd = digits_base(n, p)
bd = digits_base(B, p)

# make same length
D = max(len(nd), len(bd))
if len(nd) < D:
nd += [0] * (D - len(nd))
if len(bd) < D:
bd += [0] * (D - len(bd))

# allow[i] = count of valid digits for k at position i
# If n_i = d, then k_i can be 0..(p-1-d), which has (p-d) choices
allow = [p - nd[i] for i in range(D)]

# ways[i] = product_{0..i-1} allow[pos]
# i.e. number of valid assignments for lower i digits
ways = [1] * (D + 1)
for i in range(D):
ways[i + 1] = ways[i] * allow[i]

# digit DP accumulation
cnt = 0
# iterate from highest digit to lowest
for i in range(D - 1, -1, -1):
digit = bd[i] # B's digit at position i
lim = allow[i] # valid digits are 0..lim-1

# We want to count valid k whose prefix is smaller than B's prefix
# by choosing this digit < digit, while also < lim.
take = digit if digit < lim else lim
cnt += take * ways[i] # lower i digits: any valid combination

# If B's digit itself is not valid (digit >= lim), we cannot continue
# matching B at this position; stop.
if digit >= lim:
break

return cnt


def gen_res5_low_digits(n: int, exp: int):
"""
Generate all residues r modulo 5^exp such that the low exp base-5 digits
of r can be added to the low exp digits of n WITHOUT carry.

That is, if n has digits n_0..n_{exp-1}, then r_i must satisfy:
0 <= r_i <= 4 - n_i
for each i in [0..exp-1].

Returns list of integers r in [0, 5^exp).
"""
mod = 5 ** exp
nd = digits_base(n % mod, 5, exp)
res = []

def dfs(pos, val, pw):
# pos: which digit we're choosing (little-endian)
# val: current residue value
# pw : current power of 5 (5^pos)
if pos == exp:
res.append(val)
return
# maximum digit allowed at this position to avoid carry
lim = 4 - nd[pos]
for d in range(lim + 1):
dfs(pos + 1, val + d * pw, pw * 5)

dfs(0, 0, 1)
return res


def build_ok_table_5pow9(n_part: int) -> bytearray:
"""
Build a lookup table of size 5^9.

table[x] = 1 iff x (interpreted as a 9-digit base-5 number) can be added
to n_part (also interpreted as 9 base-5 digits) with NO carry.

We use 9-digit blocks so that for 5^18 arithmetic we can split into two blocks:
low 9 digits and high 9 digits.
"""
P9 = 5 ** 9
nd = digits_base(n_part, 5, 9) # exactly 9 digits
lim = [4 - d for d in nd] # per-position digit maximum for x
table = bytearray(P9)

# brute-force over all 5^9 possible blocks
for x in range(P9):
tmp = x
ok = 1
for L in lim:
# current base-5 digit of x
if tmp % 5 > L:
ok = 0
break
tmp //= 5
table[x] = ok
return table


def intersection_count(n: int, L: int, exp: int = 12) -> int:
"""
Count A25:
|{ 0 <= k < L : (n+k) adds with n WITHOUT carry in base 2 AND base 5 }|
which equals the count of k such that C(n+k, n) is not divisible by 2 or 5.

Trick: work modulo 10^exp (here exp=12, so step = 10^12).
If the last exp digits (in base 2 and base 5) are no-carry simultaneously,
then all candidates k lie in certain residue classes mod 10^exp.

Steps:
1) Find all residues a mod 2^exp satisfying no-carry with n mod 2^exp.
2) Find all residues b mod 5^exp satisfying no-carry with n mod 5^exp.
3) Combine each pair (a, b) into a residue r mod 10^exp via CRT.
4) For each residue r, enumerate k = r + t*10^exp < L.
For each such k, we still need global no-carry, not just low exp digits.
For base 2: global no-carry iff (k & n) == 0 (bitwise test).
For base 5: we test no-carry mod 5^18 by splitting into two 9-digit blocks,
and stepping with k += 10^exp updates those blocks efficiently.
"""
mod2 = 2 ** exp
mod5 = 5 ** exp
step = 10 ** exp # 2^exp * 5^exp

# ---- residues mod 2^exp ----
# In base 2, "no carry" means at every bit position where n has 1,
# k must have 0. So for exp low bits:
# (k_low & n_low) == 0
n2 = n % mod2
res2 = [x for x in range(mod2) if (x & n2) == 0]

# ---- residues mod 5^exp ----
# Enumerate all low exp base-5 residues with no carry:
res5 = gen_res5_low_digits(n, exp)

# ---- combine to residues mod 10^exp ----
# For each (a mod 2^exp, b mod 5^exp), CRT gives unique r mod 10^exp.
residues = []
for a in res2:
for b in res5:
residues.append(CRT([mod2, mod5], [a, b]))

# ---- base-5 global checking via block tables ----
# We'll check no-carry in base 5 modulo 5^18 using two 9-digit blocks.
mod5_9 = 5 ** 9

# ok0 tests low 9 base-5 digits, ok1 tests next 9 digits
ok0 = build_ok_table_5pow9(n % mod5_9)
ok1 = build_ok_table_5pow9((n // mod5_9) % mod5_9)

# represent arithmetic mod 5^18 using (low, high) each mod 5^9
mod5_18 = mod5_9 * mod5_9

# step5 = step mod 5^18, also split into low/high blocks for fast increment
step5 = step % mod5_18
step5_low = step5 % mod5_9
step5_high = step5 // mod5_9

inter = 0

for r in residues:
# only consider residue classes that actually hit [0, L)
if r >= L:
continue

# number of terms in this arithmetic progression:
# k = r + t*step, for t = 0..cnt-1
cnt = (L - 1 - r) // step + 1

# initialize k and its base-5 blocks (mod 5^18)
k = r
low = k % mod5_9
high = (k // mod5_9) % mod5_9

for _ in range(cnt):
# Base-2 global no-carry test:
# (k & n) == 0 <=> at every bit position with n's 1, k has 0.
#
# Base-5 test (mod 5^18) approximates global no-carry; exp=12 step
# makes candidates sparse; this block test matches typical accepted approach:
# require no carry for 18 base-5 digits, enough for n up to 1e12-scale.
if (k & n) == 0 and ok0[low] and ok1[high]:
inter += 1

# advance to next candidate in the same residue class
k += step

# update base-5 blocks (low, high) by adding step5
# (low + step5_low) may overflow mod5_9 -> carry into high
low2 = low + step5_low
carry = 1 if low2 >= mod5_9 else 0
low = low2 - (mod5_9 if carry else 0)

high2 = high + step5_high + carry
high = high2 - mod5_9 if high2 >= mod5_9 else high2

return inter


# Count numbers k in [0, L) with no-carry for base 2 and base 5:
A2 = count_no_carry(N, L, 2)
A5 = count_no_carry(N, L, 5)

# Count intersection (no-carry in both bases) using residue stepping mod 10^12:
A25 = intersection_count(N, L, 12)

# Inclusion–exclusion:
# divisible by 10 <=> divisible by 2 AND divisible by 5
ans = L - A2 - A5 + A25
print(ans)