Project Euler 455

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

Project Euler 455

题目

Powers With Trailing Digits

Let \(f(n)\) be the largest positive integer \(x\) less than \(10^9\) such that the last \(9\) digits of \(n^x\) form the number \(x\) (including leading zeros), or zero if no such integer exists.

For example:

  • \(f(4) = 411728896 (4^{411728896} = \dots490\underline{411728896})\)
  • \(f(10) = 0\)
  • \(f(157) = 743757 (157^{743757} = \dots567\underline{000743757})\)
  • \(\sum f(n), 2 \le n \le 10^3 = 442530011399\)

Find \(\sum f(n), 2 \le n \le 10^6\).

解决方案

记 \(M=10^9=2^9\cdot 5^9\)，“末 9 位等于 \(x\)”等价于模 \(M\) 的同余：\(n^x \equiv x \pmod{M}.\)则问题变为：在区间 \(1\le x<M\) 内找满足\(n^x\equiv x\pmod M\)的最大解；若无解则返回 \(0\)。

若 \(10\mid n\)，则 \(2\mid n\) 且 \(5\mid n\)。由同余式可得 \(n^x\equiv x\pmod{10}\)。但 \(10\mid n\Rightarrow n^x\equiv 0\pmod{10}\)，于是必须有\(x\equiv 0\pmod{10}.\)

同理看模 \(10^2\)：当 \(x\ge 2\) 时 \(n^x\equiv 0\pmod{100}\)，因此必须 \(x\equiv 0\pmod{100}\)。继续递推得到必须\(x\equiv 0\pmod{M}.\)

但题设要求 \(0<x<M\)，不可能成立，因此\(10\mid n\Longrightarrow f(n)=0.\)后续只需处理 \(10\nmid n\) 的 \(n\)。

因为 \(\gcd(2^9,5^9)=1\)，有经典等价：

\[ n^x\equiv x\pmod{M} \iff \begin{cases} n^x\equiv x\pmod{2^9},\\ n^x\equiv x\pmod{5^9}. \end{cases} \]

因此可以分别求出两边的解，再用中国剩余定理（CRT）拼回模 \(M\) 的解。

不过在模 \(5^9\) 的分支中，不能简单把未知数理解为 \(x\bmod 5^9\)。当 \(\gcd(n,5)=1\) 时，由欧拉定理\(n^{x+\varphi(5^k)}\equiv n^x\pmod{5^k}, \varphi(5^k)=4\cdot 5^{k-1},\)所以 \(n^x\bmod 5^k\) 实际只依赖于 \(x\bmod (4\cdot 5^{k-1})\)。在从模 \(5^k\) 提升到模 \(5^{k+1}\) 时，我们会把解写成\(x'=x+4\cdot 5^k\alpha (\alpha=0,1,2,3,4),\)从而保证 \(n^{x'}\equiv n^x\pmod{5^{k+1}}\)，并在这五个候选中唯一选出满足新模数的解。因此提升过程自然要求我们把解看作定义在模 \(4\cdot 5^k\) 下，而非仅仅模 \(5^k\)。

令 \(p\) 为素数，\(k\ge 1\)。考虑\(n^x\equiv x\pmod{p^k}.\)接下来分\(2\)种情况讨论。

\(p\mid n\)

当 \(x\ge k\) 时 \(n^x\equiv 0\pmod{p^k}\)，因此方程要求 \(x\equiv 0\pmod{p^k}\)。

在模 \(p^k(p-1)\) 的意义下，满足 \(x\equiv 0\pmod{p^k}\) 的解恰好是\(x\equiv 0,p^k,2p^k,\dots,(p-2)p^k\pmod{p^k(p-1)},\)一共 \(p-1\) 个解。

\(p\nmid n\)

先看模 \(p\)。由于费马小定理 \(n^{p-1}\equiv 1\pmod p\)，所以 \(n^x\bmod p\) 只依赖于 \(x\bmod (p-1)\)。

固定一个余数类 \(i\in\{0,1,\dots,p-2\}\)，并附加条件\(x\equiv i\pmod{p-1}.\)

则 \(n^x\equiv n^i\pmod p\)，原方程模 \(p\) 就变成\(x\equiv n^i\pmod p.\)

于是我们同时拥有

\[ \begin{cases} x\equiv i\pmod{p-1},\\ x\equiv n^i\pmod p. \end{cases} \]

由于 \(\gcd(p,p-1)=1\)，CRT 保证该系统在模 \(p(p-1)\) 下唯一可解。当 \(i\) 从 \(0\) 到 \(p-2\) 变化时，得到恰好 \(p-1\) 个不同解：

\[ x\equiv x_i(n)\pmod{p(p-1)}. \]

其中上述方程\(x_i(n)\)表示由上面的方程组通过CRT恢复而来的解。

最终，假设已知某个 \(t\) 满足\(n^t\equiv t\pmod{p^k}, (p\nmid n).\)

我们要找到 \(x\) 满足：

  1. \(x\equiv t\pmod{p^k(p-1)}\)（保持指数周期一致性）
  2. \(n^x\equiv x\pmod{p^{k+1}}\)

写成\(x=t+p^k(p-1)\alpha, \alpha\in\{0,1,\dots,p-1\}.\)因为 \(\varphi(p^{k+1})=p^k(p-1)\)，且 \(p\nmid n\)，所以\(n^{x}\equiv n^{t}\pmod{p^{k+1}}.\)

因此只需满足\(n^t \equiv t+p^k(p-1)\alpha \pmod{p^{k+1}}.\)由于 \(n^t-t\) 被 \(p^k\) 整除，令\(n^t-t=p^k\cdot c, c\in\mathbb Z.\)把上式代入并除以 \(p^k\)，在模 \(p\) 下得到\(c\equiv (p-1)\alpha\pmod p.\)

因为 \(p-1\not\equiv 0\pmod p\)，它在模 \(p\) 下可逆，所以 \(\alpha\) 唯一确定，从而提升解唯一。

结论：对 \(\gcd(n,p)=1\)，每一个模 \(p^k\) 的解都会唯一提升到模 \(p^{k+1}\)，因此模 \(p^k(p-1)\) 下总共有且只有 \(p-1\) 个解。

回到本题，我们首先考虑\(p=2\)，于是 \(p-1=1\)，总解数为 \(1\)。直接在 \(x\in[0,2^5-1]\) 枚举即可找到唯一解（按 \(n\bmod 512\) 分类预处理）。

记该唯一解为\(x\equiv a(n)\pmod{2^5}.\)

现在考虑 \(p=5\)，所以理论上应该有 \(p-1=4\) 个解，而且它们自然落在\(x \bmod{4\cdot 5^9}\)的意义下。

  • 若 \(5\mid n\)，根据，四个解就是\(b\in\{0,5^9,2\cdot 5^9,3\cdot 5^9\}\bmod{4\cdot 5^9}.\)
  • 若 \(5\nmid n\)，先求底层模 \(5\) 的解并提升到模 \(5^9\)。

对于底层解（也就是模 20的解），此时 \(p=5\)，\(p-1=4\)。对每个 \(i\in\{0,1,2,3\}\)，规定 \(x\equiv i\pmod 4\)。则\(n^x\equiv n^i\pmod 5,\)需要满足 \(x\equiv n^i\pmod 5\)。于是系统为

\[ \begin{cases} x\equiv i\pmod 4,\\ x\equiv n^i\pmod 5. \end{cases} \]

在模 \(20\) 下唯一解，得到 \(4\) 个基解。

接下来逐级提升到 \(5^9\)。从 \(k=1\)（模 \(5\)）开始，按上面的公式依次提升到 \(5^2,5^3,\dots,5^9\)，每一步唯一。最终得到 4 个解：\(x\equiv b_j(n)\pmod{4\cdot 5^9}, j=1,2,3,4.\)

我们现在得到：

\[ \begin{cases} x\equiv a \pmod{512},\\ x\equiv b \pmod{4\cdot 5^9}, \end{cases} \]

其中\(b\in\{b_1,b_2,b_3,b_4\}.\)注意 \(\gcd(512,4\cdot 5^9)=4\)。因此该系统可解的必要充分条件是\(a\equiv b\pmod 4.\)一旦满足，就有唯一解模 \(\mathrm{lcm}(512,4\cdot 5^9)=2^9\cdot 5^9=M\)。

设 \(P=5^9\)。写 \(x=a+512t\)，代入第二个同余：\(a+512t\equiv b\pmod{4P}\Rightarrow 512t\equiv b-a\pmod{4P}.\)

因为 \(b-a\) 被 \(4\) 整除，除以 \(4\) 得\(128t\equiv \frac{b-a}{4}\pmod{P}.\)且 \(\gcd(128,P)=1\)，所以可逆。令\(t\equiv \dfrac{b-a}{4}\cdot 128^{-1}\pmod P,\)

即可得到\(x\equiv a+512t\pmod{M}.\)对四个 \(b\) 分别做一次 CRT，最多得到 4 个候选解，取其中最大的 \(x\in[0,M)\) 就是 \(f(n)\)。

代码

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

#include <bits/stdc++.h>
#include "tools.h"
using namespace std;

typedef long long ll;

const int N = 1000000;
const int POW2 = 1 << 9;
const int POW5 = 1953125;
const int M = POW2 * POW5;

array<int, 4> lift_solutions_mod_5_9(int n)
{
    if (n % 5 == 0)
    {
        return {0, POW5, 2 * POW5, 3 * POW5};
    }

    int n5 = n % 5;
    int p = 1;
    int base[4];
    for (int i = 0; i < 4; i++)
    {
        int u;
        if (i == 0)
            u = 1;
        else
        {
            p = (p * n5) % 5;
            u = p;
        }
        int t = (4 * (u - i)) % 5;
        if (t < 0)
            t += 5;
        base[i] = i + 4 * t;
    }

    array<int, 4> res;
    for (int idx = 0; idx < 4; idx++)
    {
        ll t = base[idx];
        int mod = 5;
        for (int step_i = 0; step_i < 8; step_i++)
        {
            mod *= 5;
            int nt = qpow(n, t, mod);
            ll diff = nt - (t % mod);
            diff %= mod;
            if (diff < 0)
                diff += mod;
            int step = mod / 5;
            int c = ((diff / step) % 5);
            int alpha = (c * 4) % 5;
            t += 4LL * step * alpha;
        }
        res[idx] = t;
    }
    return res;
}

int crt_best(int a, const array<int, 4> &b_list, ll inv128)
{
    int best = 0;
    for (int i = 0; i < 4; i++)
    {
        int b = b_list[i];
        if (((a - b) & 3) != 0)
            continue;
        ll delta = (ll)b - a;
        ll t = (delta / 4) % POW5;
        if (t < 0)
            t += POW5;
        t = (t * inv128) % POW5;
        int x = ((a + (ll)POW2 * t) % M);
        if (x > best)
            best = x;
    }
    return best;
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    int sol2[POW2];
    for (int n = 0; n < POW2; n++)
    {
        if ((n & 1) == 0)
        {
            sol2[n] = 0;
        }
        else
        {
            int cur = 1;
            for (int x = 0; x < POW2; x++)
            {
                if (cur == x)
                {
                    sol2[n] = x;
                    break;
                }
                cur = (cur * n) & (POW2 - 1);
            }
        }
    }

    ll inv128 = mod_inverse(128, POW5);

    ll total = 0;
    for (int n = 2; n <= N; n++)
    {
        if (n % 10 == 0)
            continue;
        int a = sol2[n & (POW2 - 1)];
        auto b_list = lift_solutions_mod_5_9(n);
        total += crt_best(a, b_list, inv128);
    }

    cout << total << "\n";
    return 0;
}