luogu P3747 [六省联考 2017] 相逢是问候

题面传送门
首先你需要知道扩展欧拉定理：
当\(b\geq \varphi(p)\)时，有\(a^b\equiv a^{b\bmod \varphi(p)+\varphi(p)}\pmod p\)。
然后你需要做过两道题：P4139 上帝与集合的正确用法与P4145 上帝造题的七分钟 2 / 花神游历各国，知道扩展欧拉定理和线段树势能分析。
首先这个形式是扩展欧拉定理的形式，直接暴力计算可以做到\(O(nm^2)\)的复杂度。可以拿到50分。
考虑优化，观察样例二发现一个数很快衰减到不会变的情况。也就是说在扩展欧拉定理的时候，\(\varphi(p)\)会经过\(O(\log p)\)次后衰减为\(1\)，当其为\(1\)时显然为\(0\)，可以直接返回。
然后写一个线段树，如果一个节点里面全部不会改变了就返回，否则继续递归修改，对于每次修改暴力跑答案，时间复杂度\(O(n\log n\log ^2p)\)看上去不太能过的亚子。
考虑优化，我们发现只有\(O(\log p)\)个模数，因此可以处理光速幂，就可以做到\(O(n\log n\log p+\sqrt p\log p)\)的复杂度了。
如果真的是这样的话就非常美好与优美，但是这个题和上帝那题的区别在于：上帝那题\(2^{2^{2^{\dots}}}\)一定大于\(\varphi (p)\)，但是这道题不是，所以要特判\(b< \varphi(p)\)，这样的话讨论就会很多，非常屎山。
代码感觉只能拿来对拍了。
code:

#include<bits/stdc++.h>
#define I inline
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
#define abs(x) ((x)>0?(x):-(x))
#define ll long long
#define db double
#define lb long db
#define N (50000+5)
#define M (6000000+5)
#define K (350)
#define mod 998244353
#define Mod (mod-1)
#define eps (1e-9)
#define ull unsigned ll
#define it iterator
#define Gc() getchar() 
#define Me(x,y) memset(x,y,sizeof(x))
#define Mc(x,y) memcpy(x,y,sizeof(x))
#define d(x,y) ((k+1)*(x)+(y))
#define R(n) (1ll*rand()*rand()%(n)+1)
#define Pc(x) putchar(x)
#define LB lower_bound
#define UB upper_bound
#define PB push_back
using namespace std;
int n,m,p,c,op,x,y,A[N],TL[N][105],P[105],Ph,k[105];ll P1[105][20005],P2[105][20005];
I int GP(int p){int ToT=1;for(int i=2;i*i<=p;i++) if(p%i==0){ToT=ToT*(i-1);p/=i;while(p%i==0) ToT=ToT*i,p/=i;}p^1&&(ToT*=p-1);return ToT;}
I ll mpow(int y,int d){return P1[d][y%k[d]]*P2[d][y/k[d]]%P[d];}
I ll calc(int p,int cnt,int d,int &Fl,int Id){if(p==1){Fl=1;return 0;}if(!cnt) return A[Id]%p;return (~TL[Id][cnt-1]&&TL[Id][cnt-1]<P[d])?mpow(TL[Id][cnt-1],d-1):mpow(calc(P[d],cnt-1,d+1,Fl,Id)+P[d],d-1);}
namespace Tree{
	#define ls v<<1
	#define rs v<<1|1
	ll sum[N<<2];int Fl[N<<2],Ct[N];I void Up(int v){Fl[v]=Fl[ls]&Fl[rs];sum[v]=sum[ls]+sum[rs];}
	I void BD(int l=1,int r=n,int v=1){if(l==r) {sum[v]=A[l];return;}int m=l+r>>1;BD(l,m,ls);BD(m+1,r,rs);Up(v);}
	I ll Qry(int x,int y,int l=1,int r=n,int v=1){if(x<=l&&r<=y) return sum[v];int m=l+r>>1;return (x<=m?Qry(x,y,l,m,ls):0)+(y>m?Qry(x,y,m+1,r,rs):0);}
	I void Ins(int x,int y,int l=1,int r=n,int v=1){if(Fl[v]) return;if(l==r) return Ct[l]++,c^1?(sum[v]=calc(p,Ct[l],1,Fl[v],l)):(sum[v]=1,Fl[v]=1),void();int m=l+r>>1;x<=m&&(Ins(x,y,l,m,ls),0);y>m&&(Ins(x,y,m+1,r,rs),0);Up(v);}
	#undef ls
	#undef rs
}
int main(){
	freopen("1.in","r",stdin);
	int i,j;scanf("%d%d%d%d",&n,&m,&p,&c);P[Ph=0]=p;for(i=1;P[Ph]^1;i++) P[i]=GP(P[Ph++]);for(i=0;i<=Ph;i++){k[i]=sqrt(P[i]*2);for(P1[i][0]=j=1;j<=k[i];j++) P1[i][j]=P1[i][j-1]*c%P[i];for(P2[i][0]=j=1;j<=P[i]*2/k[i];j++) P2[i][j]=P2[i][j-1]*P1[i][k[i]]%P[i];}
	for(i=1;i<=n;i++){
		scanf("%d",&A[i]);TL[i][0]=A[i];for(j=1;j<=100;j++){if(TL[i][j-1]==-1||c==1){TL[i][j]=-1;continue;}TL[i][j]=mpow(TL[i][j-1],0);if(fabs(log(TL[i][j])/log(c)-TL[i][j-1])>eps)TL[i][j]=-1;}
	} Tree::BD();while(m--)scanf("%d%d%d",&op,&x,&y),op?printf("%lld\n",Tree::Qry(x,y)%p):(Tree::Ins(x,y),0);
}