和《花神游历各国》有异曲同工之妙。
首先能想到扩展欧拉定理：
\[a^b\equiv \begin{cases} a^{b\bmod \varphi(p)+\varphi(p)}&\text{if }b\geq\varphi(p)\\ a^b&\text{if }b< \varphi(p) \end{cases} \pmod p \]观察到一个数经过若干次操作后肯定是 \(c^{c^{\cdots^{a_i}}}\) 的形态，所以我们肯定是不断递归取对这个数取模。
结论：如果对一个数 \(p\) 不断取 \(\varphi\)，那么最多取到 \(O(\log p)\) 次这个数就会变成 \(1\)。
证明：
- 若 \(p\) 为偶数，那么 \(\leq p\) 的数中，和 \(p\) 互质的数肯定只会是奇数，故 \(\varphi(p)\leq \dfrac{p}{2}\)。
- 若 \(p\) 为奇数，那么一次取 \(\varphi\) 后 \(p\) 就会变成偶数，转回第一种情况。
所以对一个数 \(p\) 不断取 \(\varphi\) 的次数是 \(\log p\) 级别的。
所以如果一个数 \(a_i\) 经过一定次数的操作，它 \(\bmod p\) 的值就会不变。
准确地说，如果 \(p\) 取 \(cntp\) 次 \(\varphi\) 就能变成 \(1\)，那么 \(a_i\) 经过 \(cntp\) 次操作后再操作都一直是一个定值。
所以每个数的前 \(cntp\) 次操作我们暴力算即可，这个过程中我们需要实现求 \(c^x\bmod p\)。
我们可以类似分块一样，设置一个值 \(B=\sqrt{p}\)，并预处理出 \(c^{0},c^1,\cdots,c^{B}\) 模 \(p\) 的值，以及 \(c^{B},c^{2B},\cdots,c^{B\times B}\) 模 \(p\) 的值，这样就可以实现 \(O(1)\) 求 \(c^x\bmod p\) 了。
注意 \(a=0\) 需要特判，因为它要经过 \(cntp+1\) 次操作后才是一个定值。
注意暴力递归算 \(c^{c^{\cdots^{a_i}}}\) 的时候需要注意使用扩展欧拉定理时是按第一种情况还是按第二种情况。
代码如下：
#include<bits/stdc++.h>
#define LP 70
#define N 50010
#define ll long long
using namespace std;
namespace modular
{
inline int add(int x,int y,int mod){return x+y>=mod?x+y-mod:x+y;}
inline int dec(int x,int y,int mod){return x-y<0?x-y+mod:x-y;}
inline int mul(int x,int y,int mod){return 1ll*x*y%mod;}
}using namespace modular;
inline int read()
{
int x=0,f=1;
char ch=getchar();
while(ch<'0'||ch>'9')
{
if(ch=='-') f=-1;
ch=getchar();
}
while(ch>='0'&&ch<='9')
{
x=(x<<1)+(x<<3)+(ch^'0');
ch=getchar();
}
return x*f;
}
int n,m,P,c,a[N];
int cntp,p[LP];
const int B=10000;
int pow1[LP][10010],pow2[LP][20010];
bool f1[LP][10010],f2[LP][20010];
//pow1[i]:c^i mod P,pow2[i]:c^(Bi) mod P
int minn[N<<2],sum[N<<2];
inline int getphi(int n)
{
int ans=n;
for(int i=2,t=sqrt(n);i<=t;i++)
{
if(!(n%i))
{
ans=ans/i*(i-1);
while(!(n%i)) n/=i;
}
}
if(n!=1) ans=ans/n*(n-1);
return ans;
}
void init()
{
p[0]=P;
while(P!=1) p[++cntp]=P=getphi(P);
P=p[0];
for(int i=0;i<cntp;i++)
{
int np=p[i];
pow1[i][0]=1;
for(int j=1;j<=B;j++)
{
pow1[i][j]=mul(pow1[i][j-1],c,np);
f1[i][j]=(f1[i][j-1]||(1ll*pow1[i][j-1]*c>=np));
}
pow2[i][0]=1;
for(int j=1;j*B<=(np<<1);j++)
{
pow2[i][j]=mul(pow2[i][j-1],pow1[i][B],np);
f2[i][j]=(f2[i][j-1]||f1[i][B]||(1ll*pow2[i][j-1]*pow1[i][B]>=np));
}
}
}
bool flag;
inline int poww(int x,int id)
{
// assert(x<=(p[id]<<1));
ll tmp=1ll*pow1[id][x%B]*pow2[id][x/B];
flag|=f1[id][x%B]|f2[id][x/B];
flag|=(tmp>=p[id]);
return tmp%p[id];
}
int calc(int a,int numc)
{
flag=0;
int now=a;
if(numc==cntp+1)
{
now=114514;
numc--;
}
if(now>=p[numc]) flag=1,now=now%p[numc];
for(int i=numc-1;i>=0;i--)
{
if(flag) now+=p[i+1];
flag=0;
now=poww(now,i);
}
return now;
}
inline void up(int k)
{
sum[k]=add(sum[k<<1],sum[k<<1|1],P);
minn[k]=min(minn[k<<1],minn[k<<1|1]);
}
void build(int k,int l,int r)
{
if(l==r)
{
sum[k]=a[l]=read();
return;
}
int mid=(l+r)>>1;
build(k<<1,l,mid);
build(k<<1|1,mid+1,r);
up(k);
}
void update(int k,int l,int r,int ql,int qr)
{
if(minn[k]>cntp) return;
if(l==r)
{
minn[k]++;
sum[k]=calc(a[l],minn[k]);
return;
}
int mid=(l+r)>>1;
if(ql<=mid) update(k<<1,l,mid,ql,qr);
if(qr>mid) update(k<<1|1,mid+1,r,ql,qr);
up(k);
}
int query(int k,int l,int r,int ql,int qr)
{
if(ql<=l&&r<=qr) return sum[k];
int mid=(l+r)>>1,ans=0;
if(ql<=mid) ans=add(ans,query(k<<1,l,mid,ql,qr),P);
if(qr>mid) ans=add(ans,query(k<<1|1,mid+1,r,ql,qr),P);
return ans;
}
int main()
{
n=read(),m=read(),P=read(),c=read();
init();
build(1,1,n);
while(m--)
{
int opt=read(),l=read(),r=read();
if(!opt) update(1,1,n,l,r);
else printf("%d\n",query(1,1,n,l,r));
}
return 0;
}
/*
1 6 4 2
0
0 1 1
1 1 1
0 1 1
1 1 1
0 1 1
1 1 1
*/