可以先考虑对一个位置修改,如果修改了(k)次,那权值就是(c^{c^{cdots^{c^{a_i}}}}mod {p}),其中里面有(k)个(c),这里不妨设为(f_{i,k}).根据拓展欧拉定理,可得(f_{i,k}=c^{f_{i,k-1} \% varphi(p)+[f_{i,k-1}ge varphi(p)]*varphi(p)}\% p),然后(f_{i,k-1})就可以递归求.可以发现对于一个(p),求大约(O(log p))次(varphi)就会变成(1),所以对于一个位置修改若干次后他的值就会确定下来
现在考虑区间修改,我们维护区间修改次数的最小值,如果被修改区间的修改次数的最小值还没到对应阈值就暴力递归到叶子更新,由于每个位置只会更新(O(log n ))次,所以花在暴力递归上的复杂度为(O(nlog^2n))
这里的复杂度瓶颈在于每个位置的求值,每次的复杂度为两个log(递归+每次要快速幂),所以考虑预处理模数为每一种可能的模数的幂数组,例如(c^imod p_j,c^{10000i}mod p_j),就可以做到线性了
//下面是瞎写的
#include<bits/stdc++.h>
#define LL long long
#define db double
using namespace std;
const int N=50000+10;
const LL gg=1e9;
int rd()
{
int x=0,w=1;char ch=0;
while(ch<'0'||ch>'9'){if(ch=='-') w=-1;ch=getchar();}
while(ch>='0'&&ch<='9'){x=x*10+(ch^48);ch=getchar();}
return x*w;
}
int n,q,p,c,a[N],va[N][110],pp[110],pz[110][N],an[110][N],t;
int s[N<<3],mi[N<<3];
void psup(int o){s[o]=(s[o<<1]+s[o<<1|1])%p,mi[o]=min(mi[o<<1],mi[o<<1|1]);}
int fpow(int a,int b,int ii)
{
int mod=pp[ii];
if(mod<=N/2&&b<N) return an[ii][b];
int na=1;
for(int i=0;b;++i)
{
if(b&1) na=1ll*na*pz[ii][i]%mod;
b>>=1;
}
return na;
}
int pw(LL a,LL b)
{
if(b==gg) return b;
if(a==1||!b) return 1;
LL an=1;
while(b){if(b&1) an=min(an*a,gg);a=min(a*a,gg),b>>=1;}
return an;
}
int cal(int x,int y)
{
int an=x%pp[y],zk=x;
for(int i=y-1;~i;--i)
{
int nb=an+(zk>=pp[i+1]?pp[i+1]:0);
zk=pw(c,zk),an=fpow(c,nb,i);
}
return an;
}
void modif(int o,int l,int r,int ll,int rr)
{
if(mi[o]>=t) return;
if(l==r){++mi[o],s[o]=va[l][mi[o]];return;}
int mid=(l+r)>>1;
if(ll<=mid) modif(o<<1,l,mid,ll,rr);
if(rr>mid) modif(o<<1|1,mid+1,r,ll,rr);
psup(o);
}
int quer(int o,int l,int r,int ll,int rr)
{
if(ll<=l&&r<=rr) return s[o];
int an=0,mid=(l+r)>>1;
if(ll<=mid) an=(an+quer(o<<1,l,mid,ll,rr))%p;
if(rr>mid) an=(an+quer(o<<1|1,mid+1,r,ll,rr))%p;
return an;
}
void bui(int o,int l,int r)
{
if(l==r){s[o]=a[l]%p,mi[o]=0;return;}
int mid=(l+r)>>1;
bui(o<<1,l,mid),bui(o<<1|1,mid+1,r);
psup(o);
}
int main()
{
n=rd(),q=rd(),p=rd(),c=rd();
pp[t=0]=p;
while(pp[t]>1)
{
pp[++t]=1;
int x=pp[t-1],sqt=sqrt(x);
for(int j=2;j<=sqt;++j)
if(x%j==0)
{
pp[t]*=j-1,x/=j;
while(x%j==0) pp[t]*=j,x/=j;
}
if(x>1) pp[t]*=x-1;
}
pp[++t]=1;
for(int i=0;i<=t;++i)
{
if(pp[i]<=N/2)
{
an[i][0]=1%pp[i];
for(int j=1;j<N;++j) an[i][j]=1ll*an[i][j-1]*c%pp[i];
}
else
{
pz[i][0]=c%pp[i];
for(int j=1;j<40;++j) pz[i][j]=1ll*pz[i][j-1]*pz[i][j-1]%pp[i];
}
}
for(int i=1;i<=n;++i) a[i]=rd();
for(int j=0;j<=t;++j)
for(int i=1;i<=n;++i)
va[i][j]=cal(a[i],j);
bui(1,1,n);
while(q--)
{
int op=rd(),l=rd(),r=rd();
if(!op) modif(1,1,n,l,r);
else printf("%d
",quer(1,1,n,l,r));
}
return 0;
}