Solution
这道题的话。。长得就是线段树的样子qwq
如果做过的话。。可能会联想到bzoj3211（没写博qwq晚点再说吧哈哈。。）
首先大胆猜一波结论：这题跟3211一样也是修改到了一定程度就不会再有变化了！那然后写起来就是线段树暴力修改然后如果整个区间达到了修改上限的话那就不走了
然而这个上限是啥呢。。
这里还是要用到扩展欧拉定理
（然而实际上在这题里面第一条并没有什么用毕竟gcd(a,p)=1这个限制有点令人难受qwq）
然后接着我们来看一下这个修改操作
经过多次修改之后应该是长成一堆(c)指数翻上去，然后最后那个指数是(c^a_i)这个样子（有没有长得很像上帝和集合（bzoj3884）那道题的那一堆(2)?)
然后我们就用类似bzoj3884中的处理方法，大力降幂，跟那题一样，当(varphi(p))降到(1)了之后(b%varphi(1)+varphi(1)=1)，也就是说接下来无论再怎么搞都不会再变了（每个数最多被修改(log)次就不会再变了，具体为啥的话在bzoj3884里面有写就不提了）
与那题不同的是，这里需要判断指数与模数的大小关系来判断能否降成(b\% varphi(p)+varphi(p))
然而这题写起来其实。。有点恶心
具体一点的话就是：
首先我们预处理出(p)被取(varphi)多少次之后会变成(1)，记这个次数为(Mx)，并且将中间每一步的结果记录在一个数组里面，方便我们后面求值的时候直接调用作为模数
然后我们开一棵线段树，对于修改操作，我们在每个区间记录一个(sum[x])和一个(mn[x])，前者表示答案，后者表示这个区间中被修改最少的那个位置被修改了多少次
在修改的时候，如果说当前区间的(mn[x]>Mx)那么说明这个区间所有的位置都已经达到上限可以直接跳过不用管了，否则递归进去处理，一直到叶子节点。叶子节点的话，每次修改的时候能直接改的只是这个节点的(mn[x])，而单点的(sum[x])则要重新再算一次
至于怎么计算(sum[x])的话，可以考虑从降幂的最后一层层推回来，然而因为每次我们推出来的都是下一层的指数部分，如果用最暴力的快速幂求解的话，加上层数最多是(log)以及线段树的一个(log)，总的复杂度会变成(O(nlog^3n))，在TLE的边缘试探。。
然后我们注意到每次快速幂底数其实都是(c)，所以我们可以考虑将(c)的一些次方预处理出来，这样就可以直接调用了
然而这里有一个小问题，次方数。。最大可以去到(10^8)，所以这里要用一个小trick：
我们预处理出(c)的(1-10000)次方以及(c^{10000})的(1-10000)次方，这样对于一个指数，大于(10^4)的部分和小于(10^4)的部分分开来调用就好了，然后就可以十分愉快地省掉一个(log)，(O(nlog^2n))相对来说就优秀很多啦
此外还有一个比较麻烦的地方是，我们需要判断指数和(p)的大小，注意，这里的指数是不能取模的，那就有点糟糕了qwq
这里的解决方案是，注意到由于我们是倒过来推的，这里的指数在用扩展欧拉定理降幂之后会变成(c)的若干次方，所以我们可以先预处理出(c)的若干次方的值（处理到这个值(>10^9)就可以了，也可以更加小一点，反正(>10^8)就行），对于(>10^8)的部分我们直接打上一个(tag)（代码里面我是直接将表示当前幂指数的那个变量赋成(-1)），然后根据是否有(tag)判断是否能够用扩展欧拉定理的第三条就好了
细节什么的比较。。嗯qwq（调试调了不知道多久最后终于调出来了发现是判断的时候下标忘记+1也是很开心。。。）
代码大概长这个样子
#include<iostream>
#include<cstdio>
#include<cstring>
#define ll long long
using namespace std;
const int MAXN=50010,SEG=MAXN*4;
int mod[MAXN],prime[MAXN];
int pwc[10010][30],pwc1[10010][30],pc[60];
bool vis[MAXN];
int n,m,P,c,Mx,Up=10000,ans;
int a[MAXN];
namespace Seg{/*{{{*/
int ch[SEG][2],sum[SEG],mi[SEG];
int n,tot;
void pushup(int x){
sum[x]=(1LL*sum[ch[x][0]]+sum[ch[x][1]])%P;
mi[x]=min(mi[ch[x][0]],mi[ch[x][1]]);
}
void _build(int x,int l,int r){
sum[x]=0; mi[x]=0;
if (l==r) {sum[x]=a[l];return;}
int mid=l+r>>1;
ch[x][0]=++tot; _build(ch[x][0],l,mid);
ch[x][1]=++tot; _build(ch[x][1],mid+1,r);
pushup(x);
}
void build(int _n){n=_n;tot=1;_build(1,1,n);}
int powc(int Pw,int Mod){return 1LL*pwc1[Pw/Up][Mod]*pwc[Pw%Up][Mod]%mod[Mod];}
int calc(int loc,int Pw){
int nowmi=a[loc],ret=a[loc];
for (int i=Pw-1;i>=0;--i){
if (nowmi==-1||nowmi>=mod[i+1]){
ret=powc(ret%mod[i+1]+mod[i+1],i);
nowmi=-1;
}
else{
ret=powc(ret,i);
if (nowmi<60) nowmi=pc[nowmi];
else nowmi=-1;
}
}
return ret;
}
void _update(int x,int l,int r,int lx,int rx){
if (mi[x]>=Mx) return;
if (lx==rx){
++mi[x];
sum[x]=calc(lx,mi[x]);
return;
}
int mid=lx+rx>>1;
if (l<=mid) _update(ch[x][0],l,r,lx,mid);
if (r>mid) _update(ch[x][1],l,r,mid+1,rx);
pushup(x);
}
void update(int l,int r){_update(1,l,r,1,n);}
int _query(int x,int l,int r,int lx,int rx){
if (l<=lx&&rx<=r) return sum[x];
int mid=lx+rx>>1,ret=0;
if (l<=mid) ret=(1LL*ret+_query(ch[x][0],l,r,lx,mid))%P;
if (r>mid) ret=(1LL*ret+_query(ch[x][1],l,r,mid+1,rx))%P;
return ret;
}
int query(int l,int r){return _query(1,l,r,1,n);}
}/*}}}*/
void prework(int n);
int phi(int x);
int ksm(int x,int y,int p);
int main(){
#ifndef ONLINE_JUDGE
freopen("a.in","r",stdin);
#endif
int op,l,r;
scanf("%d%d%d%d",&n,&m,&P,&c);
for (int i=1;i<=n;++i) scanf("%d",a+i);
prework(MAXN);
Seg::build(n);
for (int i=1;i<=m;++i){
scanf("%d%d%d",&op,&l,&r);
if (op==0)
Seg::update(l,r);
else{
ans=Seg::query(l,r);
printf("%d
",ans);
}
}
}
void prework(int n){
int cnt=0;
for (int i=2;i<n;++i){
if (!vis[i])
prime[++cnt]=i;
for (int j=1;j<=cnt&&prime[j]*i<n;++j){
vis[i*prime[j]]=1;
if (i%prime[j]==0) break;
}
}
mod[0]=P;Mx=0;
while (mod[Mx]!=1) ++Mx,mod[Mx]=phi(mod[Mx-1]);
mod[++Mx]=1;
for (int i=0;i<=Up;++i)
for (int j=0;j<=Mx;++j)
pwc[i][j]=ksm(c,i,mod[j]);
for (int i=0;i<=Up;++i)
for (int j=0;j<=Mx;++j)
pwc1[i][j]=ksm(ksm(c,Up,mod[j]),i,mod[j]);
//...?
pc[0]=1;
int flag=100000;
for (int i=1;i<60;++i){
if (pc[i-1]==-1) pc[i]=-1;
if (1LL*pc[i-1]*c<=1000000000)
pc[i]=pc[i-1]*c;
else
pc[i]=-1,flag=min(flag,i);
}
}
int phi(int x){
int ret=x;
for (int i=1;prime[i]*prime[i]<=x;++i){
if (x%prime[i]) continue;
ret=ret-ret/prime[i];
while (x%prime[i]==0) x/=prime[i];
}
if (x>1) ret=ret-ret/x;
return ret;
}
int ksm(int x,int y,int p){
int ret=1,base=x;
for (;y;y>>=1,base=1LL*base*base%p)
if (y&1) ret=1LL*ret*base%p;
return ret;
}