[bzoj4869] [loj#2142] [Shoi2017] 相逢是问候

题意简述

题面

Informatik verbindet dich und mich.
信息将你我连结。

维护一个长度为 (n) 的数组 (a[])，支持两个操作:

1. 对 (i in [l,r]) ，进行替换 (a[i] gets c^{a[i]})
2. 求 (sumlimits_{i=l}^r a[i]) ((mod) (P))

其中 (c,P) 为给定常数。
共操作 (m) 次。
(n,mleq 50000,P leq 10^8)

想法

看到 (c^{a[i]}\%P) 可以想到扩展欧拉定理：

[egin{equation*} c^xequiv egin{cases} c^x (mod P)& x<varphi(P)\ c^{x\%varphi(P)+varphi(P)} (mod P)& x geq varphi(P) end{cases} end{equation*} ]

递推一下有：
(c^{c^x}equiv c^{varphi(P)+c^x\% P} equiv c^{varphi(P)+c^{varphi(varphi(P))+x\%varphi(varphi(P))}\% varphi(P)})

指数上的 (varphi(varphi(varphi(...)))) 在经过 (Theta(log_2 P)) 后就会恒定为 1
（证明：若 (P) 为偶数，则 (varphi(P)leq P/2) ；若 (P) 为奇数，则 (varphi(P)) 为偶数）
由此可知，对每个 (a[i])，操作1执行 (Theta(log_2 P)) 遍后其值就恒定不变了。

可以预处理出第二次 (varphi(varphi(varphi(...)))=1) 时 (varphi()) 的个数 (w)。
（为何是第二次？见大佬博客）

建立线段树，每个节点记录 (sum) 与它所对应区间中每个值是否都修改过 (w) 次以上（即是否都已恒定不变）。
询问操作就是常规求和。
暴力进行修改操作，若某节点中所有值都恒定不变则不用再修改。
对于需要修改的叶子节点，通过扩展欧拉定理计算修改后的值。

这样的话找到待修改的点的复杂度为 (O(nlogn)) ，修改的复杂度是 (O(log^2n)) (快速幂需要 (O(logn)))
总复杂度是 (O(log^3n)) ，可能会超时。

考虑优化快速幂。
发现我们要计算的底数都是 (c) ，于是用类似大步小步法优化。
假设我们要算的是 (c^x (mod P)) ，令 (s=sqrt{P},x=p imes s+r, p,r<P)，则 (c^x=(c^s)^p imes c^r)
预处理出 (c^s) 的幂与 (c) 的幂即可。

最后注意一大堆乱七八糟的细节。

总结

技巧

1. 像不断开方或此题这种奇怪的不好维护的操作，可以思考会不会很少的修改次数内成为定值，这样就可以暴力修改。

2. 类大步小步法优化同底数快速幂。

3. 对扩展欧拉定理的应用：判断指数与 (varphi(P)) 大小要在各种有取模操作的位置判断。
   引用 (Sengxian) 大佬之言：

指数循环节公式只在 (xgeq varphi(n)) 时成立，在 (UVa 10692) 中，用试乘来判断是否 (geq varphi(n))，我们在试乘的时候，是以上一层返回的取模后结果作为幂试乘，这样并不准确，应该使用上一层的答案进行试乘。但是放心，经过验证，这样做没有任何问题。因为如果 (xge varphi(n))，那么 (a^xgeq n) 只在 (n = 6) 时不成立，经过验证，这个带来的一系列后续影响不会造成答案的错误，所以大可放心使用。

误区

1. (a^{b^c}=a^{(b^c)} eq (a^b)^c=a^{bc})

2. 扩展欧拉定理 (c^xequiv c^{x\%varphi(P)+varphi(P)} (mod P)) 的使用条件是 (x geq varphi(P)) ！
   并不是普适！要注意判断是否要加上 (varphi(P))

手残

1. 线段树建树到叶子节点时，不要弄混改点的存储位置与它代表的有实际意义的点的坐标（(x) 与 (l)）。

2. 暴力找质因子是边界条件为 (i imes i leq x) ，别忘了可取等！

代码

#include<cstdio>
#include<iostream>
#include<algorithm>

using namespace std;

int read(){
	int x=0;
	char ch=getchar();
	while(!isdigit(ch)) ch=getchar();
	while(isdigit(ch)) x=x*10+ch-'0',ch=getchar();
	return x;
}

const int N = 50005;

int n,m,P,c,a[N],w;
int phi[36],mul[20005][36],bml[20005][36],is[20005][36],bis[200005][36];

int getphi(int x){
	int y=x,ret=x;
	for(int i=2;i*i<=x;i++) /**/
		if(y%i==0){
			ret=ret/i*(i-1);
			while(y%i==0) y/=i;
 		}
	if(y!=1) ret=ret/y*(y-1);
	return ret;
}
int flag;
int Pow(int x,int y) { /**/
	flag=(1ll*mul[x%20000][y]*bml[x/20000][y]>=phi[y]);
	flag=max(flag,max(is[x%20000][y],bis[x/20000][y])); /**/
	return 1ll*mul[x%20000][y]*bml[x/20000][y]%phi[y];
} 

int root,cnt,sum[N*2],ch[N*2][2],mn[N*2];
void update(int x){
	sum[x]=(sum[ch[x][0]]+sum[ch[x][1]])%P;
	mn[x]=min(mn[ch[x][0]],mn[ch[x][1]]);
}
void build(int x,int l,int r){
	if(l==r) { sum[x]=a[l]%P;/**/ mn[x]=0; return; }
	int mid=(l+r)>>1;
	build(ch[x][0]=++cnt,l,mid);
	build(ch[x][1]=++cnt,mid+1,r);
	update(x);
}
void modify(int x,int l,int r,int L,int R){
	if(l==r){
		mn[x]++;
		if(mn[x]>w) return;
		flag=a[l]>=phi[mn[x]];/**/
		sum[x]=a[l]%phi[mn[x]]; /**/
		for(int i=mn[x];i>=1;i--){
			if(flag) sum[x]=Pow(sum[x]+phi[i],i-1);
			else sum[x]=Pow(sum[x],i-1);
		}
		return;
	}
	if(mn[x]>=w) return;
	int mid=(l+r)>>1;
	if(L<=l && r<=R){
		modify(ch[x][0],l,mid,L,R);
		modify(ch[x][1],mid+1,r,L,R);
	}
	else{
		if(L<=mid) modify(ch[x][0],l,mid,L,R);
		if(R>mid) modify(ch[x][1],mid+1,r,L,R);
	}
	update(x);
}
int ask(int x,int l,int r,int L,int R){
	if(L<=l && r<=R) return sum[x];
	int ret=0,mid=(l+r)>>1;
	if(L<=mid) ret=(ret+ask(ch[x][0],l,mid,L,R))%P;
	if(R>mid) ret=(ret+ask(ch[x][1],mid+1,r,L,R))%P;
	return ret;
}

int main()
{
	n=read(); m=read(); P=read(); c=read();
	for(int i=1;i<=n;i++) a[i]=read();
	
	phi[0]=P;
	for(int i=1;i<36;i++){
		phi[i]=getphi(phi[i-1]);
		if(phi[i]==1 && phi[i-1]==1) { w=i; break; }
	}
	mul[0][0]=1; is[0][0]=mul[0][0]>=phi[0];
	for(int i=1;i<=20000;i++) {
		is[i][0]=max(is[i-1][0],1ll*mul[i-1][0]*c>=phi[0]?1:0);
		mul[i][0]=1ll*mul[i-1][0]*c%P;
	}
	bml[0][0]=1; bis[0][0]=bml[0][0]>=phi[0];
	for(int i=1;i<=20000;i++) {
		bis[i][0]=max(bis[i-1][0],1ll*bis[i-1][0]*mul[20000][0]>=phi[0]?1:0);
		bml[i][0]=1ll*bml[i-1][0]*mul[20000][0]%P;
	}
	for(int i=1;i<=w;i++){
		mul[0][i]=1; is[0][i]=mul[0][i]>=phi[i];
		for(int j=1;j<=20000;j++) {
			is[j][i]=max(is[j-1][i],1ll*mul[j-1][i]*c>=phi[i]?1:0);
			mul[j][i]=1ll*mul[j-1][i]*c%phi[i];
		}
		bml[0][i]=1; bis[0][i]=bml[0][i]>=phi[i];
		for(int j=1;j<=20000;j++) {
			bis[j][i]=max(bis[j-1][i],1ll*bis[j-1][i]*mul[20000][i]>=phi[i]?1:0);
			bml[j][i]=1ll*bml[j-1][i]*mul[20000][i]%phi[i];
		}
	}
	
	build(root=++cnt,1,n);
	
	int opt,l,r;
	while(m--){
		opt=read(); l=read(); r=read();
		if(opt==0) modify(root,1,n,l,r);
		else printf("%d
",ask(root,1,n,l,r)%P);
	}
	
	return 0;
}