RT.目测复杂度是对的（？）

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N=5e4+5;
const int SZ=1e4;
const int LG=60+5;
// his[i][k]表示ai的所有历史版本,其中初始k=0 cnt[i]表示第i个位置被操作了多少次
LL a[N],his[N][LG],pw1[SZ+5][LG],pw2[SZ+5][LG],bit[N];
// DIV[i]表示第i格 第一个>=p的历史版本 其中ai为第0个历史版本
int fa[N],cnt[N],DIV[N],nPhi,n,q;
LL MOD,c,phi[N];
inline int LSB(int x){return x&-x;}
void add(int u,LL x)
{
    x=(x%MOD+MOD)%MOD;
    while(u<=n){(bit[u]+=x)%=MOD;u+=LSB(u);}
}
LL psq(int x)
{
    LL ans=0;
    while(x){(ans+=bit[x])%=MOD;x-=LSB(x);}
    return ans;
}
int find(int u){return u==fa[u]?u:fa[u]=find(fa[u]);}
LL PHI(LL x)
{
    LL ans=x;
    for(int i=2;i*i<=x;++i)
        if(x%i==0)
        {
            ans=ans*(i-1)/i;
            while(x%i==0) x/=i;
        }
    if(x>1) ans=ans*(x-1)/x;
    return ans;
}
void prework_phi()
{
    phi[0]=MOD;
    nPhi=0;
    while(phi[nPhi]>1)
    {
        ++nPhi;
        phi[nPhi]=PHI(phi[nPhi-1]);
    }
    phi[nPhi+1]=1;
}
void prework_pw()
{
    // pw1表示c^1...c^10000
    // pw2表示c^10000 ... c^1e8
    for(int i=0;i<=nPhi;++i)
    {
        pw1[0][i]=1;
        for(int j=1;j<=SZ;++j)
            pw1[j][i]=pw1[j-1][i]*c%phi[i];
        pw2[0][i]=1;
        for(int j=1;j<=SZ;++j)
            pw2[j][i]=pw2[j-1][i]*pw1[SZ][i]%phi[i];
    }
}
// 返回 c^n mod phi[id]
LL qpow(LL n,int id){return pw1[n%SZ][id]*pw2[n/SZ][id]%phi[id];}
inline LL read()
{
    LL x=0; bool f=false;
    char c=getchar();
    while(!isdigit(c)){if(c=='-')f=true;c=getchar();}
    while(isdigit(c)){x=x*10+c-'0';c=getchar();}
    return f?-x:x;
}
int main()
{
    ios::sync_with_stdio(false);
    n=read(); q=read(); MOD=read(); c=read();
    prework_phi();
    prework_pw();
    LL MID=0,X=c;
    while(X<MOD) ++MID,X*=c; // MID = log(p)/log(c)
    for(int i=1;i<=n;++i)
    {
        a[i]=read(); his[i][0]=a[i];
        add(i,a[i]);
        fa[i]=i;
    }
    fa[n+1]=n+1;
    while(q--)
    {
        int op,l,r;
        op=read();l=read();r=read();
        if(op==0)
        {
            for(int i=find(l);i<=r;i=find(i+1))
            {
                if(DIV[i]==0) // 值<=p
                {
                    if(his[i][cnt[i]]>=MID) DIV[i]=cnt[i]+1; // 从cnt+1个历史版本开始有ai>p
                    // 求解c^a[i] mod p 即 a[i] mod phi(p)
                    int j;
                    for(j=cnt[i];j&&cnt[i]-j+1<=nPhi;--j) if(his[i][j]<=phi[cnt[i]-j+1]) break;
                    LL x=qpow(his[i][j],cnt[i]-j);
                    ++j;
                    for(;j<=cnt[i];++j)
                        x=qpow(x+phi[cnt[i]-j+1],cnt[i]-j);
                    add(i,x-his[i][cnt[i]]);
                    ++cnt[i];
                    his[i][cnt[i]]=x;
                }
                else //从<=p的值开始找
                {
                    int j;
                    for(j=DIV[i]-1;j&&cnt[i]-j+1<=nPhi;--j) if(his[i][j]<=phi[cnt[i]-j+1]) break;
                    if(cnt[i]-j+1>nPhi&&his[i][j]!=0){fa[i]=i+1;continue;}
                    LL x=qpow(his[i][j],cnt[i]-j);
                    ++j;
                    for(;j<=cnt[i];++j)
                        x=qpow(x+phi[cnt[i]-j+1],cnt[i]-j);
                    add(i,x-his[i][cnt[i]]);
                    ++cnt[i];
                    his[i][cnt[i]]=x;
                    
                }
            }
            
        }
        else 
            cout << ((psq(r) - psq(l-1))%MOD+MOD)%MOD << "\n";
    }
    return 0;
}