幂塔与欧拉降幂 - Semi-Puzzle: Brain Storm - 牛客多校Day9

Semi-Puzzle: Brain Storm

前往 算竞题解 分类 || 前往 算竞模板 分类

原题链接

B-Semi-Puzzle: Brain Storm 2023牛客暑期多校训练营9

简明题意

给定 $a(1 \le a \le 10^9)$ 和 $m(1 \le m \le 10^9)$，求 $u$ 使得 $a^u \equiv u \pmod{m}$，要求 $1 \le u \le 10^{18}$。

多组测试数据，$T \le 1000$。

前置芝士

在开始之前你必须要知道以下三个知识点，否则就别想按着题解思路做了。会的应该也不需要此题解了。

我也是边学边写的这道题，过程可以说很艰辛，所以作此文整理一下。

扩展欧拉定理

这个东西通常用于给快速幂降幂，但是确是下一个知识点的前置。

欧拉函数

​定义： 在区间 $[1, x)$ 中与 $x$ 互质的数的个数，记作 $\phi(x)$。

​计算： 设分解质因数 $x=\prod{p_i^{\alpha_i}}$，则：

$$\phi(x)=x\prod\frac{p_i-1}{p_i}$$

扩展欧拉定理

对于任意数 $a$ 和 $p$，如果 $gcd(a,p)=1$，满足欧拉定理：$a^{\phi(p)} \equiv 1 \pmod{p}$。

如对任意数 $a$ 和 $p$，可得扩展欧拉定理，又称欧拉降幂：

$$a^x \equiv \left\{ \begin{array}{ll} a^{x\%\phi(p)} & ,gcd(a,p)=1 \\ a^x & ,gcd(a,p)\neq1 & ,x<\phi(p) \\ a^{x\%\phi(p)+\phi(p)} & ,gcd(a,p)\neq1 & ,x\ge\phi(p) \\ \end{array} \right. \pmod{p}$$

我们可以归纳出一个粗糙且更实用的公式：

$$设~x~\text{exmod}~p=(x-p)\%p+p~, \\ a^x \equiv a^{x~\text{exmod}~\phi(p)} \pmod{p}$$

其中 $\phi(x)$ 表示欧拉函数。

带模幂塔问题

这算是一个经典的扩展欧拉定理的问题，证明了幂塔的值只与前有限个项有关。

这个问题的经典例题是：Codeforces Round 906 - D: Power Tower

带模幂塔问题

​定义： 将幂作为指数不断迭代得到的表达式，形如：

$$x = a_1^{\left(a_2^{\left(a_3^{\cdot^{\cdot^\cdot}}\right)}\right)}$$

​带模幂塔问题与欧拉降幂：

因为幂塔的增长非常快，我们通常处理的是带模幂塔问题，必须使用欧拉降幂，而这其中又有很好的性质。

$$\begin{array}{l} 设~x=x_1,~x_i \equiv a_i^{x_{i+1}} \pmod{p}~, \\ 设~x~\text{exmod}~p=(x-p)\%p+p~, \\ 设~\phi_1(x)=\phi(x),~\phi_i(x)=\phi(\phi_{i-1}(x))~, \\ \end{array}$$

$$\begin{array}{rll} & x_1 \equiv a_1^{x_2} & (\text{mod}~~{p}) \\ ~~~~\Longrightarrow~~~~ & x_1 \equiv a_1^{x_2~\text{exmod}~\phi(p)} & (\text{mod}~~{p}) \\ ~~~~\Longrightarrow~~~~ & x_2 \equiv a_2^{x_3~\text{exmod}~\phi_2(p)} & (\text{exmod}~~\phi(p)) \\ ~~~~\Longrightarrow~~~~ & x_3 \equiv a_3^{x_4~\text{exmod}~\phi_3(p)} & (\text{exmod}~~\phi_2(p)) \\ ~~~~\Longrightarrow~~~~ & \cdots \\ \end{array}$$

我们发现通过不断递归，最后 $\phi_m(p)=1$，那么再往后的递归都是同余的，即：

$$\forall m \le n,~\phi_m(p)=0,~a_1^{\cdot^{\cdot^{\cdot^{a_n}}}} \equiv a_1^{\cdot^{\cdot^{\cdot^{a_{m-1}}}}} \pmod{p}$$

根据欧拉函数定义，$m$ 一定是 $\ln{p}$ 级别的，使用递归求解即可，注意将全过程模重载为 $\text{exmod}$。

Miller Rabin 素性测试与 Pollard Rho 算法

这是两个关于素数问题的随机算法。Miller Rabin 可以在 $O(k\log_2n)$ 的时间内判断 $n$ 是否为素数，其中 $k$ 表示测试次数，测试次数越多，算法正确率越高。Pollard Rho 在 $n$ 不为素数时，可在 $O(n^{1/4})$ 的时间内找到 $n$ 的一个因数（$1$ 除外），注意不一定是质因数。

这两种算法都有随机性，但实际表现均远高于理论性能，前往 Luogu 自行搜索 Pollard Rho 学习即可。这里还想谴责一下，出题人在题解丝毫没有提及就在标程随手贴了一套没学过的 Rho 算法，否则会 TLE 飞，对本萌新不是很友好。

题解思路

根据带模幂塔问题，令 $A_i \equiv a^{A_{i-1}}$，最后有 $A_{i-1} \equiv A_i \equiv a^{A_{i-1}}$，所以只需不断递归 $A_i$ 在 $log_2m$ 步级别可以得到 $A_i \equiv a^{A_i}$。

在这个式子中我们需要特别注意 $A_i \equiv a^{A_i}$ 中是不能直接将模 $m$ 后的 $A_i$ 拿来用的，因为根据欧拉降幂，等式应该是 $A_i~\text{mod}~m=a^{A_i~\text{exmod}~\phi(m)}~\text{mod}~m$，其中 $x~\text{exmod}~m$ 表示 $(x-m)\%m+m$。但是设 $d = \text{lcm}(m,\phi(m))$，我们发现全过程模 $d$ 等式依然成立，注意由于欧拉降幂这里的模建议都使用 $\text{exmod}$。于是乎我们只需要求出第一个 $A_i \equiv a^{A_i} \pmod{\text{lcm}(m,\phi(m))}$ 即可。

特别注意，在求幂塔的过程中，由于 $m$ 很大，涉及乘法需要使用 $int128$，求欧拉函数需要使用 Pollard Rho 算法。

最后要有一个特判，因为要求输出 $1 \le u \le 10^{18}$ 而我们最后得到的 $\text{exmod}$ 后可以达到 $2 \times 10^{18}$，因此我们特判 $\text{mod}$ 后是否成立即可。题解提到了这样做的正确性，但没有给出证明，有懂哥可以在评论区留言。

代码

请自行跳转 156 行，有效代码 20 行。

// Semi-Puzzle: Brain Storm
// Module Version 230418
// Powered by GreatLiangpi
#ifdef ONLINE_JUDGE
#define TEST 1 // Test Switch
#else
// #define TEST 1
#define LOCAL 1 // Debug Switch
#endif
#ifndef LOCAL
#define TEST 1
#endif
#include <bits/stdc++.h>
#define int int64_t
#ifdef TEST
#pragma GCC optimize(3, "Ofast", "inline")  // Optimization Switch
#define endl '\n'
#endif
#ifdef LOCAL
#define endl " **\n"
#endif
#define tomax(x, y) x = max(x, (y))
#define tomin(x, y) x = min(x, (y))
using namespace std;
#ifdef TEST
const int MAX = 200005;
#else
#ifdef LOCAL
const int MAX = 105;
#endif
#endif
const int INF = 0x3f3f3f3f3f3f3f3fll;
const int MOD = 1000000007;
int call_count = 0;

namespace Prime_Tests {
// 快速幂, 该快速幂使用 mod
int quick_pow(__int128_t b, int e, int m) {
    __int128_t res = 1 % m;
    b = b % m;
    while (e) {
        if (e & 1) res = res * b % m;
        b = b * b % m;
        e >>= 1;
    }
    return res;
}

// Miller Rabin 素性测试
// 154590409516822759, 2305843009213693951, 4384957924686954497
int test[] = {
    2,  3,  5,  7,  11,
    13, 17, 19, 23, 29,
    31, 37, 41, 43, 47,
    79, 83, 89, 97, 233,
};  // 测试集: 优先使用小质数
bool Miller_Rabin(int n) {
    if (n < 2 or n > 2 and n % 2 == 0)
        return false;
    int s = n - 1;
    while (not(s & 1)) s >>= 1;
    for (int p : test) {
        if (n == p)
            return true;
        __int128_t t = s, m = quick_pow(p, s, n);
        while (t != n - 1 and m != 1 and m != n - 1) {
            m = m * m % n;
            t <<= 1;
        }
        if (m != n - 1 and not(t & 1))
            return false;
    }
    return true;
}

// 随机数生成器
mt19937 rnd((unsigned int)chrono::steady_clock::now().time_since_epoch().count());
int randint(int l, int r = 0) {
    if (l > r) swap(l, r);
    return rnd() % (r - l + 1) + l;
}

// Pollard Rho 因数测试
int Pollard_Rho(int n) {
    if (n % 2 == 0)  // Pollard 不能判 4
        return 2;
    if (Miller_Rabin(n))
        return n;
    while (true) {
        int c = randint(1, n - 1);
        auto nxt = [&](__int128_t x) -> int { return (x * x + c) % n; };
        int turtle = 1, rabbit = -1;
        __int128_t mul = 1;
        while (turtle != rabbit) {
            for (int i = 0; i < 128; i++) {
                turtle = nxt(turtle);
                rabbit = nxt(nxt(rabbit));
                __int128_t tmp = mul * abs(turtle - rabbit) % n;
                // 龟兔相遇或积将为 0 退出
                if (turtle == rabbit or tmp == 0)
                    break;
                mul = tmp;
            }
            int g = __gcd((int)mul, n);
            if (g > 1)
                return g;
        }
    }
}

// 使用 Pollard Rho 的分解质因数
int divisors[128];
int prime_divisors(int x) {
    int cnt = 0, p = 0;
    divisors[cnt++] = x;
    while (p < cnt) {
        int x = divisors[p];
        int tmp = Pollard_Rho(x);
        if (tmp == x) {
            p++;
        } else {
            divisors[p] = tmp;
            divisors[cnt++] = x / tmp;
        }
    }
    sort(divisors, divisors + cnt);
    cnt = unique(divisors, divisors + cnt) - divisors;
    return cnt;
}
}

// 使用 Pollard Rho 的欧拉函数
int euler(int x) {
    if (x < 2) return 0;
    int cnt = Prime_Tests::prime_divisors(x);
    int res = x;
    for (int i = 0; i < cnt; i++) {
        int p = Prime_Tests::divisors[i];
        res = res / p * (p - 1);
    }
    return res;
}

// 其实这样定义 exmod 更加清晰, 重点是跑的还快一点...
#define exmod(x, m) ((x) < (m) ? (x) : (x) % (m) + (m))

// 快速幂, 该快速幂使用 exmod
int quick_pow(__int128_t b, int e, int m) {
    __int128_t res = exmod(1, m);
    b = exmod(b, m);
    while (e) {
        if (e & 1) res = exmod(res * b, m);
        b = exmod(b * b, m);
        e >>= 1;
    }
    return res;
}

// 计算幂塔
int calc(int a, int m) {
    if (m == 1) return exmod(a, m);
    return quick_pow(a, calc(a, euler(m)), m);
}

void solve() {
    int a, m;
    cin >> a >> m;
    int phi = euler(m);
    int lcm = m / __gcd(m, phi) * phi;
    int res = calc(a, lcm);
    if (Prime_Tests::quick_pow(a, res % lcm, m) == res % lcm % m) {
        cout << res % lcm << endl;
    } else {
        cout << res << endl;
    }
}

int32_t main() {
#ifdef TEST
    // cin.tie(0);
    // cout.tie(0);
    // ios::sync_with_stdio(false); 
#endif
#ifdef LOCAL
    // freopen("B-Data\\3.in", "r", stdin);
    // freopen("out.txt", "w", stdout);
    clock_t start_time = clock();
#endif
    cout << fixed << setprecision(2);


    int t = 1;
    cin >> t;
    while (t--)
        solve();


#ifdef LOCAL
    cout << "Used " << call_count << " Function Calls." << endl;
    cout << "Used " << (int)(clock() - start_time) << " Microseconds of Time." << endl;
#endif 
    return 0;
}