Project Euler 714

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

Project Euler 714

题目

Duodigits

We call a natural number a duodigit if its decimal representation uses no more than two different digits.

For example, \(12\), \(110\) and \(33333\) are duodigits, while \(102\) is not.

It can be shown that every natural number has duodigit multiples. Let \(d(n)\) be the smallest (positive) multiple of the number \(n\) that happens to be a duodigit. For example, \(d(12)=12\), \(d(102)=1122\), \(d(103)=515\), \(d(290)=11011010\) and \(d(317)=211122\).

Let \(\displaystyle D(k)=\sum_{n=1}^k d(n)\). You are given \(D(110)=11\,047\), \(D(150)=53\,312\) and \(D(500)=29\,570\,988\).

Find \(D(50\,000)\). Give your answer in scientific notation rounded to \(13\) significant digits (\(12\) after the decimal point). If, for example, we had asked for \(D(500)\) instead, the answer format would have been \(2.957098800000e7\).

解决方案

考虑长度为 \(L\) 的 某个 duodigit \(w\)，使用的两种数字是 \(a<b\)。令：\(R_L = \underbrace{111\cdots 1}_{L\text{ times}} = \dfrac{10^L-1}{9}\)

再构造一个长度同为 \(L\) 的 0/1 十位数字 \(x\)：当 \(w\) 的某位取 \(b\) 时，\(x\) 的该位取 \(1\)；取 \(a\) 时取 \(0\)。

那么逐位写出：若该位是 \(a\)，这个数位贡献 \(a\)；若该位是 \(b\)：这个数位贡献 \(a+(b-a)=b\)，所以\(w = a\cdot R_L + (b-a)\cdot x\)恒成立。

记 \(e=b-a\in\{1,2,\dots,9\}\)，得到\(w = aR_L + e x\)

因此，只要能找到某个长度 \(L\)、某个 \(a<b\)、某个 0/1 数字 \(x\)，使得 \(aR_L+ex\) 是 \(n\) 的倍数，就得到一个 duodigit 倍数。

要求：\(aR_L + e x \equiv 0 \pmod n\)，设：\(r_a \equiv aR_L \pmod n\)，则需要：\(e x \equiv -r_a \pmod n\)。计算这个式子的困难点在于：\(x\) 是十进制“0/1 数字”，不是二进制，所以不能简单枚举。

我们可以使用动态扩展余数表法解决这个问题。任何长度为 \(L\) 的 0/1 十进制数都可以写成：\(x = d_{L-1}\cdot 10^{L-1} + \cdots + d_1\cdot 10 + d_0, d_i\in\{0,1\}\)。为了方便，我们定义 0/1 十进制数的集合为\(X\)。

把长度从 \(L\) 扩展到 \(L+1\)，只需要考虑新增最高位为 1 的情况：\(x' = x + 10^L\)。因为新增最高位为 0 等价于前面补 0，数值不变，不必生成。

我们维护：\(\text{rm}[r] = \min\{x \mid x\equiv r \pmod{n},x\in X，初始化可以从短的 0/1 数开始（长度\)$）：\(0,1,10,11\)。每次长度增加，把所有已有 \(x\) 变为 \(x+10^L\)，得到一批新的余数。

因为 \(e=b-a\le 9\) 极小，可以不做扩展欧几里得/逆元。因此，我们希望找到一个 整数代表 \(q\in[0,n)\)，满足：

\[ q \equiv x \pmod n, e q \equiv -r_a \pmod n \]

等价于存在某个 \(j\) 使：\(e q = jn - r_a\)。因为 \(1\le j\le e\) 时：\(0 < jn - r_a \le en\)，从而\(q = \dfrac{jn-r_a}{e}\)。

只要 \(jn-r_a\) 能被 \(e\) 整除，\(q\) 就落在 \([0,n]\) 内（边界 \(q=n\) 视作余数 0）。

于是对每个 \(a<b\)，只需枚举：\(e=b-a,j=1,2,\dots,e\)。若：\((jn-r_a)\bmod e = 0\)就得到候选余数 \(q\)，若 \(\text{rm}[q]\) 存在，则构造出 duodigit：

\[ w = aR_L + e\cdot \text{rm}[q] \]

可见，字典里 \(\text{rm}[0]=0\) 永远存在，但当 \(a=0\) 且需要 \(x\equiv 0\pmod n\) 时，不能用 \(x=0\)（会得到 0，不是正整数）。

因此需要额外维护x0pos：目前已发现的、满足 \(x\equiv 0\pmod n\) 的 最小正 0/1 数。一旦在扩展过程中首次得到这样的候选值，之后继续增加位数只会使 0/1 数的数值更大，因此只需保留最小的那个即可。

代码

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
N = 50000
def d(n: int) -> int:
    """
    Return d(n): the smallest positive multiple of n that is a duodigit.
    """

    # 对 n < 102，n 本身一定是 duodigit：
    # 1..99 都最多两种数字；100,101 也满足。
    if n < 102:
        return n

    # rm[r] = minimal 0/1-duodigit x with x % n == r (leading zeros allowed)
    rm = {0: 0, 1 % n: 1, 10 % n: 10, 11 % n: 11}

    # rep = R_L = 111...1 (current length L), bas = 10^(L-1)
    rep = 11
    bas = 10
    rep_mod = rep % n
    bas_mod = bas % n

    # smallest positive 0/1-multiple of n found so far
    x0pos = None

    while True:
        best = None

        # Search all digit pairs (a,b) for current length
        for a in range(10):
            ra = (a * rep_mod) % n

            # case: only digit a is used -> a * repunit
            if a and ra == 0:
                cand = a * rep
                if best is None or cand < best:
                    best = cand

            for b in range(a + 1, 10):
                e = b - a  # e in 1..9

                # find q = (j*n - ra)/e for some j in [1..e]
                for j in range(1, e + 1):
                    num = j * n - ra
                    q, r = divmod(num, e)
                    if r:
                        continue
                    if q == n:  # treat boundary as remainder 0
                        q = 0

                    # For a==0 and q==0, must use positive x ≡ 0 mod n
                    if a == 0 and q == 0:
                        x = x0pos
                    else:
                        x = rm.get(q)

                    if x is None:
                        continue

                    cand = a * rep + e * x
                    if cand == 0:
                        continue
                    if best is None or cand < best:
                        best = cand

        if best is not None:
            return best

        # Increase length L -> L+1
        bas *= 10
        bas_mod = (bas_mod * 10) % n
        rep = rep * 10 + 1
        rep_mod = (rep_mod * 10 + 1) % n

        # Update rm by adding the new most significant digit 1 (+ bas)
        snapshot = list(rm.items())
        new = {}
        for r, x in snapshot:
            r2 = r + bas_mod
            if r2 >= n:
                r2 -= n
            v2 = x + bas

            # record smallest positive 0/1 multiple
            if r2 == 0 and v2 != 0:
                if x0pos is None or v2 < x0pos:
                    x0pos = v2

            old = rm.get(r2)
            if (old is None or v2 < old) and (r2 not in new or v2 < new[r2]):
                new[r2] = v2

        rm.update(new)


def D(N: int) -> int:
    s = 0
    for n in range(1, N + 1):
        s += d(n)
    return s


def sci13(x: int) -> str:
    """
    Scientific notation with 13 significant digits (12 after decimal), rounded.
    Output example: 2.957098800000e7
    """
    if x == 0:
        return "0.000000000000e0"

    s = str(x)
    exp = len(s) - 1

    # take first 13 digits and round by next digit
    if len(s) <= 13:
        main = int(s) * (10 ** (13 - len(s)))
    else:
        main = int(s[:13])
        if ord(s[13]) >= ord("5"):
            main += 1

    # handle carry like 9.999... -> 10.000...
    if main == 10**13:
        main //= 10
        exp += 1

    ms = str(main).zfill(13)
    return f"{ms[0]}.{ms[1:]}e{exp}"


ans = D(N)
print(sci13(ans))