Project Euler 474

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

Project Euler 474

题目

Last digits of divisors

For a positive integer \(n\) and digits \(d\), we define \(F(n, d)\) as the number of the divisors of \(n\) whose last digits equal \(d\).

For example, \(F(84, 4) = 3\). Among the divisors of \(84 (1, 2, 3, 4, 6, 7, 12, 14, 21, 28, 42, 84)\), three of them \((4, 14, 84)\) have the last digit \(4\).

We can also verify that \(F(12!, 12) = 11\) and \(F(50!, 123) = 17888\).

Find \(F(10^6!, 65432)\) modulo \((10^{16} + 61)\).

解决方案

令\(N=10^6,d=65432,a=v_2(d),b=v_5(d),k=5\)，即\(k\)是\(d\)的十进制位数。并记\(g=2^a5^b.\)

首先证明引理1：设 \(p\) 为素数，\(k\ge 1\)，且 \(t=v_p(d)<k\)。若\(x\equiv d\pmod{p^k},\)则必有\(v_p(x)=t.\)

证明： 写 \(d=p^t u\)，其中 \(p\nmid u\)。由 \(x\equiv d\pmod{p^k}\) 得 \(p^k\mid (x-d)\)，特别地 \(p^{t+1}\mid(x-d)\)，于是\(x=p^t u+p^{t+1}q=p^t(u+pq).\)

括号内对 \(p\) 不为 \(0\)（因为 \(u\not\equiv 0\pmod p\)），故 \(v_p(x)=t\)。证毕。

在本题中有\(a<k, b<k.\) 因此对任意满足 \(x\equiv d\pmod{10^k}\) 的 \(x\)，由于同时有 \(x\equiv d\pmod{2^k}\) 与 \(x\equiv d\pmod{5^k}\)，由上面的引理得到\(v_2(x)=a, v_5(x)=b.\)

这说明：所有符合末位条件的因子，都必须包含完全相同的 \(2\) 与 \(5\) 因子部分 \(g\)。

令\(M=\dfrac{10^k}{g}, D=\dfrac{d}{g}.\)把因子写成\(x=g\cdot u.\)

接下来证明引理2：在 \(a<k,b<k\) 的前提下，有

\[ x\equiv d\pmod{10^k} \Longleftrightarrow \begin{cases} x=g\cdot u,\\ u\equiv D\pmod M,\\ \gcd(u,10)=1. \end{cases} \]

证明：

  • “\(\Rightarrow\)”：由引理1可知 \(v_2(x)=a,v_5(x)=b\)，故 \(x\) 可写为 \(g\cdot u\) 且 \(v_2(u)=v_5(u)=0\)，即 \(\gcd(u,10)=1\)。又 \(x\equiv d\pmod{10^k}\) 且 \(g\mid 10^k\)、\(g\mid d\)，两边同除以 \(g\) 得 \(u\equiv D\pmod M\)。
  • “\(\Leftarrow\)”：若 \(x=g u\) 且 \(u\equiv D\pmod M\)，则 \(x\equiv gD=d\pmod{gM}=10^k\)。证毕。

因此我们只需统计：从 \(N!\) 中除去 \(2,5\) 的部分后，能组成多少个满足\(u\equiv D\pmod M, \gcd(u,M)=1\)的因子。

对每个素数 \(p\le N\)，根据勒让德公式，令\(\displaystyle{e_p=v_p(N!)=\sum_{t\ge 1}\left\lfloor\frac{N}{p^t}\right\rfloor.}\) 那么任意因子可写成\(\displaystyle{x=\prod_{p\le N} p^{\alpha_p}, 0\le \alpha_p\le e_p.}\)

由引理2约化，\(u\) 只允许使用 \(p\ne 2,5\) 的素数：

\[ u=\prod_{\substack{p\le N\\p\ne 2,5}} p^{\alpha_p}, 0\le \alpha_p\le e_p, \]

并要求 \(u\bmod M\) 等于指定的单位余数 \(D\)。

记单位集合

\[ U=(\mathbb Z/M\mathbb Z)^\times={r\in{1,\dots,M-1}:\gcd(r,M)=1}. \]

由 \(\gcd(u,M)=1\) 可知最终余数一定落在 \(U\)，因此状态只需开在 \(U\) 上。

设当前从小到大已经处理完的素数集合为 \(\mathcal P\)，定义状态

\[ f_{\mathcal P}(r) =\#\left\{(\alpha_p)_{p\in\mathcal P}:\ 0\le \alpha_p\le e_p,\ \prod_{p\in\mathcal P}p^{\alpha_p}\equiv r\pmod M\right\},r\in U. \]

表示处理完\(\mathcal P\)后，各个不同模数\(r\)下因子的个数。

那么，初始条件为

\[ f_{\varnothing}(r)= \left\{\begin{aligned} &1 && \text{if}\quad r\equiv 1\pmod M,\\ &0 && \text{if}\quad r\not\equiv 1\pmod M, \end{aligned}\right. \quad (r\in U). \]

则对任意 \(r\in U\) 有：

\[ f_{\mathcal{P}\cup \{p\}}(r)= \displaystyle\sum_{t=0}^{e_p} f_{\mathcal P}(r\cdot p^{-t})\ \]

其中 \(p^{-t}\) 表示模 \(M\) 下的逆元幂；由于 \(\gcd(p,M)=1\)，逆元总是存在的。不过，难点在于 \(e_p\) 很大，不能直接做 \(O(|U|\cdot e_p)\) 的卷积。

由此，我们定义置换\(T_p(r)\equiv r\cdot p\pmod M (r\in U),\)并取它的某个循环\(\mathcal C=(c_0,c_1,\dots,c_{m-1}), c_{j+1}\equiv T_p(c_j)\equiv c_j\cdot p\pmod M,\)下标按模 \(m\) 计算（即 \(c_{j+m}=c_j\)）。

将 \(f_{\mathcal P}\) 限制到该循环上，记\(g(j)=f_{\mathcal P}(c_j), g'(j)=f_{\mathcal P'}(c_j).\)由 \(c_j\cdot p^{-t}=c_{j-t}\)，转移在循环上化为\(\displaystyle{g'(j)=\sum_{t=0}^{C-1} g(j-t).}\)

于是，我们先定义窗口和 \(w(j)\)，再写 \(g'(j)\)。令\(C=qm+r, 0\le r<m,\)并定义\(\displaystyle{S=\sum_{i=0}^{m-1} g(i).}\)

定义长度为 \(r\) 的窗口和函数\(\displaystyle{w(j)=\sum_{i=0}^{r-1} g(j-i).}\)那么 \(w(j)\) 可用以下分段递推计算：

\[ w(j)= \left\{\begin{aligned} &\sum_{i=0}^{r-1} g(-i) && \text{if}\quad j=0,\\ &w(j-1)+g(j)-g(j-r) && \text{if}\quad 1\le j\le m-1. \end{aligned}\right. \]

说明：\(g(-i)\)、\(g(j-r)\) 的下标都按模 \(m\) 解释，例如 \(g(-i)=g(m-i)\)。

因为长度为 \(C\) 的环形窗口包含 \(q\) 次整圈加 \(r\) 个额外项，所以\(g'(j)= qS+w(j).\)

最终，对循环 \(\mathcal C\) 上的每个 \(j\)，最后回代为\(f_{\mathcal{P}\cup\{p\}}(c_j)=g'(j).\)

对所有循环做同样计算，即完成一次素数 \(p\) 的整体更新。

代码

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
#include <bits/stdc++.h>
#include "tools.h"
using namespace std;

typedef long long ll;
const int N = 1000000;
const ll d0 = 65432;
const ll MOD = 10000000000000061ll; // 10^16 + 61

ll mulmod(ll a, ll b) {
    __int128_t t = ( __int128_t )a * b;
    t %= MOD;
    return (ll)t;
}


ll legendre_exponent(int N, int p)
{
    ll e = 0;
    while (N)
    {
        N /= p;
        e += N;
    }
    return e;
}

int main()
{

    // k = digits(d), L = 10^k
    ll tmp = d0;
    ll L = 1;
    int k = 0;
    while (tmp > 0)
    {
        tmp /= 10;
        L *= 10;
        ++k;
    }

    // vp2, vp5 of d
    ll t2 = d0, t5 = d0;
    int vp2 = 0, vp5 = 0;
    while ((t2 & 1LL) == 0)
    {
        t2 >>= 1;
        ++vp2;
    }
    while (t5 % 5 == 0)
    {
        t5 /= 5;
        ++vp5;
    }

    // // 本题 d=65432, k=5, vp2=3<5, vp5=0<5
    // // 这种情况下 v2,v5 被同余条件“锁死”，可以直接约掉 g。
    // if (vp2 >= k || vp5 >= k)
    // {
    //     // Euler 474 这组数据不会走到这里；若要泛化需另外处理“指数可>=k”的情况。
    //     cerr << "Unhandled case: vp2>=k or vp5>=k\n";
    //     return 0;
    // }

    ll g = 1;
    for (int i = 0; i < vp2; ++i)
        g *= 2;
    for (int i = 0; i < vp5; ++i)
        g *= 5;

    int M = (int)(L / g);  // 12500
    int D = (int)(d0 / g); // 8179
    // D 与 M 互素（单位）
    if (std::gcd(D, M) != 1)
    {
        cout << 0 << "\n";
        return 0;
    }

    // 枚举单位元素（不被2或5整除），建立 residue -> index 映射
    vector<int> units;
    vector<int> idx(M, -1);
    units.reserve(M);
    for (int r = 1; r < M; ++r)
    {
        if ((r & 1) && (r % 5 != 0))
        {
            idx[r] = (int)units.size();
            units.push_back(r);
        }
    }
    const int PHI = (int)units.size(); // 5000

    // dp[pos]：当前已处理素数的情况下，余数=units[pos] 的方案数
    vector<ll> dp(PHI, 0), ndp(PHI);
    dp[idx[1]] = 1;

    // 预分配，避免频繁分配
    vector<int> nxt(PHI);
    vector<int> mark(PHI, 0);
    int token = 1;
    vector<int> cyc(PHI);

    // 筛素数并做DP
    auto primes = Factorizer(N).primes;
    for (int p : primes)
    {
        if (p == 2 || p == 5)
            continue;

        ll e = legendre_exponent(N, p);
        ll C = e + 1; // 0..e 一共 C 种

        // 置换：pos -> pos' where residue' = residue * p mod M
        for (int i = 0; i < PHI; ++i)
        {
            int r = (int)((ll)units[i] * p % M);
            nxt[i] = idx[r]; // 一定存在（r 仍是单位）
        }

        ++token;
        // 分解成若干个循环，每个循环用滑窗做更新
        for (int i = 0; i < PHI; ++i)
        {
            if (mark[i] == token)
                continue;

            int len = 0;
            int x = i;
            do
            {
                mark[x] = token;
                cyc[len++] = x;
                x = nxt[x];
            } while (x != i);

            ll S = 0;
            for (int t = 0; t < len; ++t)
                S = (S + dp[cyc[t]]) % MOD;

            ll q = C / len;
            int r = (int)(C % len);
            ll base = mulmod(q, S);

            if (r == 0)
            {
                for (int t = 0; t < len; ++t)
                    ndp[cyc[t]] = base;
            }
            else
            {
                // window(t)=sum_{i=0..r-1} dp[cyc[(t-i) mod len]]
                ll win = 0;
                for (int j = 0; j < r; ++j)
                {
                    int id = (len - j) % len; // 0, len-1, len-2, ...
                    win = (win + dp[cyc[id]]) % MOD;
                }
                ndp[cyc[0]] = (base + win) % MOD;

                for (int t = 1; t < len; ++t)
                {
                    win = (win + dp[cyc[t]]) % MOD;
                    int old = t - r;
                    if (old < 0)
                        old += len;
                    ll sub = dp[cyc[old]];
                    if (win < sub)
                        win += MOD;
                    win -= sub;
                    ndp[cyc[t]] = (base + win) % MOD;
                }
            }
        }

        dp.swap(ndp);
    }

    cout << dp[idx[D]] % MOD << "\n";
    return 0;
}