分析
扩展欧拉定理：
$$ \displaystyle a^b\bmod p=\begin{cases}a^{b\bmod \varphi(p)+\varphi(p)}\bmod p&,b\geq \varphi(p)\\a^b\bmod p&,b<\varphi(p)\end{cases} $$
可以发现 $c^{c^{c^{c^\cdots}}}$ 至多 $\log$ 层后指数的模数就会变成 $1$ ，于是可以 $O(\log n)$ 计算这个东西。
于是修改直接暴力修改到叶子节点，如果修改次数已经超过最大值了就不用改直接返回。
具体的，稍微修改一下快速幂，使得它能够判断扩展欧拉定理的条件，然后返回对应的值。
这样子时间复杂度是 $O(n\log^3n)$ 的，Luogu 上只能拿到 90 分。
注意到 $p,\varphi(p),\varphi(\varphi(p)),\cdots,1$ 这些模数的个数很少，可以考虑对于每个模数暴力预处理出 $c^k$ 以及 $c^{10000k}$ ，然后就可以做到 $O(1)$ 快速幂了。当然要对应处理出 $c^k$ 和 $c^{10000k}$ 的值是否超过了模数，便于使用欧拉定理。
具体实现比较麻烦，建议先写一个 $O(n\log^3n)$ 的再改成 $O(n\log^2 n)$ 。
代码
// ===================================
// author: M_sea
// website: http://m-sea-blog.com/
// ===================================
#include <algorithm>
#include <iostream>
#include <cstdlib>
#include <cstring>
#include <cstdio>
#include <cmath>
#define re register
using namespace std;
inline int read() {
int X=0,w=1; char c=getchar();
while (c<'0'||c>'9') { if (c=='-') w=-1; c=getchar(); }
while (c>='0'&&c<='9') X=X*10+c-'0',c=getchar();
return X*w;
}
const int N=50000+10,M=50+10,LIM=10000+10;
int n,m,p,c,maxx;
int a[N],phi[M];
int pw1[LIM][M],pw2[LIM][M],flag1[LIM][M],flag2[LIM][M];
inline pair<int,int> qpow(int b,int id) {
int p=b%10000,q=b/10000;
int flag=0,res=1ll*pw1[p][id]*pw2[q][id]%phi[id];
if (1ll*pw1[p][id]*pw2[q][id]>=phi[id]) flag=1;
flag|=flag1[p][id]|flag2[q][id];
return make_pair(res,flag);
}
inline int getphi(int x) {
int res=x,m=x;
for (re int i=2;i*i<=x;++i) {
if (m%i) continue;
res=1ll*res*(i-1)/i;
while (m%i==0) m/=i;
}
if (m^1) res=1ll*res*(m-1)/m;
return res;
}
inline void init() {
int tmp=p; phi[0]=p;
while (tmp^1) tmp=getphi(tmp),phi[++maxx]=tmp;
phi[++maxx]=1;
for (re int i=0;i<=maxx;++i) {
pw1[0][i]=1;
for (re int j=1;j<=10000;++j) {
if (1ll*pw1[j-1][i]*c>=phi[i]) flag1[j][i]=1;
pw1[j][i]=1ll*pw1[j-1][i]*c%phi[i];
flag1[j][i]|=flag1[j-1][i];
}
}
for (re int i=0;i<=maxx;++i) {
pw2[0][i]=1,flag2[1][i]=flag1[10000][i];
for (re int j=1;j<=10000;++j) {
if (pw2[j-1][i]*pw1[10000][i]>=phi[i]) flag2[j][i]=1;
pw2[j][i]=1ll*pw2[j-1][i]*pw1[10000][i]%phi[i];
flag2[j][i]|=flag2[j-1][i];
}
}
}
inline int calc(int w,int now,int lim) {
if (now==lim) {
if (w>phi[now]) return w%phi[now]+phi[now];
else return w;
}
int t=calc(w,now+1,lim); pair<int,int> s=qpow(t,now);
if (s.second) return s.first+phi[now];
else return s.first;
}
#define ls (o<<1)
#define rs (o<<1|1)
int sumv[N<<2],addv[N<<2];
inline void build(int o,int l,int r) {
if (l==r) { sumv[o]=a[l]; return; }
int mid=(l+r)>>1;
build(ls,l,mid),build(rs,mid+1,r);
sumv[o]=(sumv[ls]+sumv[rs])%p;
}
inline void modify(int o,int l,int r,int ql,int qr) {
if (addv[o]>=maxx) return;
if (l==r) {
++addv[o],sumv[o]=calc(a[l],0,addv[o]);
return;
}
int mid=(l+r)>>1;
if (ql<=mid) modify(ls,l,mid,ql,qr);
if (qr>mid) modify(rs,mid+1,r,ql,qr);
sumv[o]=(sumv[ls]+sumv[rs])%p,addv[o]=min(addv[ls],addv[rs]);
}
inline int query(int o,int l,int r,int ql,int qr) {
if (ql<=l&&r<=qr) return sumv[o];
int mid=(l+r)>>1,res=0;
if (ql<=mid) res=(res+query(ls,l,mid,ql,qr))%p;
if (qr>mid) res=(res+query(rs,mid+1,r,ql,qr))%p;
return res;
}
#undef ls
#undef rs
int main() {
n=read(),m=read(),p=read(),c=read();
for (re int i=1;i<=n;++i) a[i]=read();
init(); build(1,1,n);
while (m--) {
int op=read(),l=read(),r=read();
if (!op) modify(1,1,n,l,r);
else printf("%d\n",query(1,1,n,l,r));
}
return 0;
}