Problem Statement
维护一个长度为 (n) 的序列，给定一个 (p) 和 (c) ，有 (m) 个操作：
0 l r将所有的 (iin [l, r]) ，进行 (a_i leftarrow c ^ {a_i}) .1 l r求出 (sum_{i = l} ^ r a_i) .
Constraints
(1le n, m le 5 imes 10 ^ 4) ，(0 < p le 10 ^ 8) ，(0 < c < p) ，(0le a_i < p) .
Editorial
由于是数据结构题，我们就会碰到这种东西：
[c ^ {c ^ {... ^ {c ^ {a_i}}}}
]
看到这种很鬼畜的一堆幂叠起来的操作，考虑扩展欧拉定理
[a^bequiv egin{cases} a^{b mod varphi(m)}, &gcd(a,m) = 1, \ a^b, &gcd(a,m)
e 1, b < varphi(m), \ a^{(b mod varphi(m)) + varphi(m)}, &gcd(a,m)
e 1, b ge varphi(m). end{cases} pmod m
]
这样的东西有一个很常见的暴力 dfs 算法，由于 (pleftarrow varphi(p)) 的操作会在 (mathcal O (log p)) 的计算次数后 (p) 变成 (1) ，所以直接递归求解的复杂度是可以接受的。
也就是说，随着 (c) 的幂次堆叠的越来越多，最后的答案在超过 (mathcal O (log p)) 次后永远是一样的。
那么，这个序列上的值最多只会发生 (mathcal O (nlog p)) 的修改，每次修改的值可以通过 dfs 预处理出来，区间求和用树状数组维护，需要修改的数字用并查集或 std::set 维护。
然后脑子和浆糊一样的我直接写完了，但是这个 dfs 递归层数为 (mathcal O (log p)) ，对于每个位置要求 (mathcal O (log p)) 个，每次快速幂要 (mathcal O (log p)) ，直接 (mathcal O (nlog ^3 p)) 。
预处理 TLE 了！！！
然后开始毒瘤，由于递归层数和总个数是不能变的，所以能优化的只有快速幂了。
然后你掏出来一个光速幂，由于真正需要的模数只有 (mathcal O (log p)) 个，成功的把预处理复杂度降为 (mathcal O (sqrt plog p + nlog ^ 2 p)) 。
unordered_map<int, int> Phi;
inline int phi(int p) {
if(Phi.count(p)) return Phi[p];
int pp = p; Phi[pp] = p;
for (int i = 2; i * i <= p; ++i) {
if(p % i == 0) {
Phi[pp] = Phi[pp] / i * (i - 1);
while(p % i == 0) p /= i;
}
}
if(p > 1) Phi[pp] = Phi[pp] / p * (p - 1);
return Phi[pp];
}
int PA[39][PS], PB[39][PS], mod[39], tot, c, B;
bool BA[39][PS], BB[39][PS];
inline void init(int p) {
int pp = p; mod[0] = p;
while(pp ^ 1) {
mod[++tot] = phi(pp), pp = mod[tot];
}
B = sqrt(2 * p); int SS = 2 * p / B + 3;
forn(k,0,tot) {
PA[k][0] = PB[k][0] = 1;
BA[k][0] = BB[k][0] = 0;
forn(i,1,B) {
BA[k][i] = BA[k][i - 1] | (1ll * PA[k][i - 1] * c >= mod[k]);
PA[k][i] = 1ll * PA[k][i - 1] * c % mod[k];
}
BB[k][1] = BA[k][B], PB[k][1] = PA[k][B];
forn(i,2,SS) {
BB[k][i] = BB[k][i - 1] | (1ll * PB[k][i - 1] * PB[k][1] >= mod[k]);
PB[k][i] = 1ll * PB[k][i - 1] * PB[k][1] % mod[k];
}
}
}
int n, m, p, a[N], b[40][N];
inline void Mod(int& A, const int& M) {(A >= M) && (A -= M);}
struct BIT {
int val[N], R, A;
inline void upd(int o, int k) {while(o <= R) Mod(val[o] += k, p), o += o & -o; }
inline int qry(int o) {A = 0; while(o) Mod(A += val[o], p), o -= o & -o; return A;}
} ZT;
// inline int euler_pow(int p, int k, int M) {
// if(k == 0) return (1 >= M) ? (1 % M + M) : 1;
// bool fl = (p >= M); int res = 1;
// (fl) && (p -= M);
// for(; k; k >>= 1, (1ll * p * p >= M) ? (fl = 1, p = 1ll * p * p % M) : (p = p * p))
// if(k & 1) (1ll * res * p >= M) ? (fl = 1, res = 1ll * p * res % M) : (res = res * p);
// return res + fl * M;
// }
// inline int solve(const int& num, int lim, int M) {
// if(lim == 0) return (num >= M) ? (num % M + M) : num;
// if(M == 1) return (c >= M) ? M : c;
// int subres = solve(num, lim - 1, phi(M));
// return euler_pow(c, subres, M);
// }
inline int FAST_POW(int k, int id) {
return 1ll * PA[id][k % B] * PB[id][k / B] % mod[id] + (BA[id][k % B] | BB[id][k / B] | (1ll * PA[id][k % B] * PB[id][k / B] >= mod[id])) * mod[id];
}
inline int FAST_solver(const int& num, int id, int lim) {
if(id == lim) return (num >= mod[id]) ? (num % mod[id] + mod[id]) : num;
if(mod[id] == 1) return (c >= mod[id]) ? mod[id] : c;
int subres = FAST_solver(num, id + 1, lim);
return FAST_POW(subres, id);
}
int lim[N], val[N]; set<int> P; vector<int> del;
inline void solve() {
Rdn(n, m, p, c), ZT.R = n, Phi[1] = 1;
init(p);
forn(i,1,n) {
Rdn(a[i]), ZT.upd(i, a[i]);
forn(j,0,27) b[j][i] = FAST_solver(a[i], 0, j) % p;
lim[i] = 27;
while(lim[i] && b[lim[i]][i] == b[lim[i] - 1][i]) lim[i] -- ;
P.insert(i);
}
while(m--) {
int opt, l, r;
Rdn(opt, l, r);
if(opt) Wtn((ZT.qry(r) - ZT.qry(l - 1) + p) % p, '
');
else {
del.clear();
auto ll = P.lower_bound(l);
auto rr = P.upper_bound(r);
for (auto it = ll; it != rr; it++) {
ZT.upd(*it, p - b[val[*it]][*it]);
ZT.upd(*it, b[val[*it] + 1][*it]);
if(++val[*it] >= lim[*it]) del.push_back(*it);
}
for (auto pos : del) P.erase(pos);
}
}
}