思维难度作为一道省选题还是有待商榷，但是代码确实挺恶心的。

记一下这种有关无穷层幂嵌套（无穷幂塔）的套路。

思路

扩展欧拉定理 + 线段树。

首先看到不断嵌套幂并且模数较大，优先考虑扩展欧拉定理。

首先 \(\gcd(a, p) = 1\) 的情况因为 \(a^{\varphi(p)} \equiv 1 \pmod p\) ，可以按照 \(\gcd(a, p) \neq 1\) 的情况处理。

现在的问题就是不断对一个区间嵌套幂，然后区间求和。

• 结论：对于正整数 \(x\)，不断令 \(x = \varphi(x)\)，嵌套的次数不超过 \(O(\log x)\) 次。

考虑对 \(a_i\) 进行操作，等价于令 \(a_i = c^{a_i} = c^{a_i \bmod \varphi(p) + \varphi(p)}\)，也就是递归计算 \(c^{a_i \bmod \varphi(p)} = c^{a_i \bmod \varphi(\varphi(p)) + \varphi(\varphi(p))}\)。显然递归一次之后 \(a_i \leq \varphi(p)\)，则根据结论知这样的总层数是 \(O(\log)\) 级别的。

考虑最后一层递归有 \(\varphi(p) = 1\)，所以回带 \(O(\log)\) 层之后得到的一定是定值。换言之，对于任意正整数，对其进行 \(O(\log)\) 次修改之后，再任意修改得到的都是定值。

于是每次修改时判断一下修改次数，不足就暴力修改即可。均摊下来时间复杂度是 \(O(n \log^3 n)\).

考虑用光速幂，最终的时间复杂度就是 \(O(n \log^2 n)\).

代码

#include <cstdio>
#include <algorithm>
using namespace std;

typedef long long ll;

const int maxn = 1e5 + 5;
const int blk_sz = 1e4 + 5;
const int blk = 1e4;
const int sgt_sz = maxn << 2;
const int lg_sz = 60;

int n, m, len;
int a[maxn];
ll p, c;
ll phi[lg_sz], pw1[blk_sz][lg_sz], pw2[blk_sz][lg_sz];
bool chk_th, vis1[blk_sz][lg_sz], vis2[blk_sz][lg_sz];

ll get_phi(ll x)
{
    ll res = x, tmp = x;
    for (int i = 2; i * i <= x; i++)
    {
        if (tmp % i == 0)
        {
            res = res / i * (i - 1);
            while (tmp % i == 0) tmp /= i;
        }
    }
    if (tmp > 1) res = res / tmp * (tmp - 1);
    return res;
}

void init()
{
    ll cur = p;
    phi[len = 0] = p;
    while (cur != 1) cur = get_phi(cur), phi[++len] = cur;
    phi[++len] = 1;
    for (int i = 0; i <= len; i++)
    {
        pw1[0][i] = 1;
        for (int j = 1; j <= blk; j++)
        {
            pw1[j][i] = pw1[j - 1][i] * c;
            if (pw1[j][i] >= phi[i]) pw1[j][i] %= phi[i], vis1[j][i] = true;
            vis1[j][i] |= vis1[j - 1][i];
        }
    }
    for (int i = 0; i <= len; i++)
    {
        pw2[0][i] = 1;
        vis2[1][i] = vis1[blk][i];
        for (int j = 1; j <= blk; j++)
        {
            pw2[j][i] = pw2[j - 1][i] * pw1[blk][i];
            if (pw2[j][i] >= phi[i]) pw2[j][i] %= phi[i], vis2[j][i] = true;
            vis2[j][i] |= vis2[j - 1][i];
        }
    }
}

ll calc(ll pw, int idx)
{
    chk_th = false;
    int v1 = pw % blk, v2 = pw / blk;
    ll res = pw1[v1][idx] * pw2[v2][idx];
    if (res >= phi[idx]) res = res % phi[idx], chk_th = true;
    chk_th |= (vis1[v1][idx] | vis2[v2][idx]);
    return res;
}

ll dfs(ll v, int dep, int lim)
{
    chk_th = false;
    if (dep == lim)
    {
        if (v >= phi[dep]) v %= phi[dep], chk_th = true;
        return v;
    }
    ll c = dfs(v, dep + 1, lim);
    return calc(chk_th ? c + phi[dep + 1] : c, dep);
}

namespace SGT
{
    #define ls (k << 1)
    #define rs (k << 1 | 1)

    int cnt[sgt_sz];
    ll sum[sgt_sz];

    void push_up(int k)
    {
        sum[k] = sum[ls] + sum[rs];
        if (sum[k] >= p) sum[k] -= p;
        cnt[k] = min(cnt[ls], cnt[rs]);
    }

    void build(int k, int l, int r)
    {
        if (l == r) return sum[k] = a[l], cnt[k] = 0, void();
        int mid = (l + r) >> 1;
        build(ls, l, mid), build(rs, mid + 1, r), push_up(k);
    }

    void update(int k, int l, int r, int ql, int qr)
    {
        if (cnt[k] >= len) return;
        if (l == r) return cnt[k]++, sum[k] = dfs(a[l], 0, cnt[k]), void();
        int mid = (l + r) >> 1;
        if (ql <= mid) update(ls, l, mid, ql, qr);
        if (qr > mid) update(rs, mid + 1, r, ql, qr);
        push_up(k);
    }

    ll query(int k, int l, int r, int ql, int qr)
    {
        if ((l >= ql) && (r <= qr)) return sum[k];
        int mid = (l + r) >> 1; ll res = 0;
        if (ql <= mid) res = query(ls, l, mid, ql, qr);
        if (qr > mid) res += query(rs, mid + 1, r, ql, qr);
        if (res > p) res -= p;
        return res;
    }
}

int main()
{
//     freopen("1.in", "r", stdin);
//     freopen("1.out", "w", stdout);
    scanf("%d%d%d%d", &n, &m, &p, &c);
    for (int i = 1; i <= n; i++) scanf("%d", &a[i]);
    init();
    SGT::build(1, 1, n);
    while (m--)
    {
        int opt, l, r;
        scanf("%d%d%d", &opt, &l, &r);
        if (!opt) SGT::update(1, 1, n, l, r);
        else printf("%lld\n", SGT::query(1, 1, n, l, r) % p);
    }
    return 0;
}