[六省联考2017]相逢是问候（扩展欧拉定理+预处理幂）

https://www.luogu.com.cn/problem/P3747

已经记不请上一次遇到扩展欧拉定理是什么时候了。

(a^b~mod~p=)
(b<phi(p),a^b~mod~p)
(bge phi(p),a^{b~mod~phi(p)+phi(p)}~mod~p)

假设模数一直取(phi)的序列为(p[1],p[2],…,p[p0](p0 le 2log_2(mo)))

对于这题，可以对每个数记录一个(cnt)表示已经pow了多少次，当(cnt>p0)后，都和(cnt=p0+1)时的值一样

注意是和(cnt=p0+1)时的值一样，为了方便，多开一位：p[++p0]=1

用线段树维护区间(cnt)最小值，如果(le p0)，就下去修改就好了。

每个数要需要修改(log)次，每次修改要(log^2)的时间求新的答案（一个log来自快速幂），这样会T两个点。

注意到底数都是(c),所以预处理(c^{0..15000})和(c^{0..15000*15000})在(mod)每个(p[i])下的值，就可以(O(2))快速幂。

Code：

#include<bits/stdc++.h>
#define fo(i, x, y) for(int i = x, _b = y; i <= _b; i ++)
#define ff(i, x, y) for(int i = x, _b = y; i <  _b; i ++)
#define fd(i, x, y) for(int i = x, _b = y; i >= _b; i --)
#define ll long long
#define pp printf
#define hh pp("
")
using namespace std;

const int N = 50005;

int n, m, sp, c;
int v[N];

int phi(int n) {
	if(n == 1) return 1;
	int s = n;
	for(int x = 2; x * x <= n; x ++) if(n % x == 0) {
		s = s / x * (x - 1);
		for(; n % x == 0; n /= x);
	}
	if(n > 1) s = s / n * (n - 1);
	return s;
}

int p[N], p0;

namespace sub1 {
	struct P {
		ll x; int y;
		P(ll _x = 0, int _y = 0) {
			x = _x, y = _y;
		}
	};
	P mul(P a, P b, int mo) {
		a.y |= b.y;
		a.x *= b.x;
		if(a.x >= mo) a.x %= mo, a.y = 1;
		return a;
	}
	
	const int M = 15000;
	
	P t1[100][M], t2[100][M];
	
	void build() {
		fo(i, 1, p0) {
			t1[i][0] = t2[i][0] = P(1 % p[1], 1 >= p[1]);
			P w = P(c % p[i], c >= p[i]);
			fo(j, 1, M) t1[i][j] = mul(t1[i][j - 1], w, p[i]);
			w = t1[i][M];
			fo(j, 1, M) t2[i][j] = mul(t2[i][j - 1], w, p[i]);
		}
	}
	
	P b[N];
	ll calc(int *a, int a0) {
		a0 = min(a0, p0);
		b[a0] = P(a[a0] % p[a0], a[a0] >= p[a0]);
		fd(i, a0 - 1, 1) {
			int y = b[i + 1].x + (b[i + 1].y ? p[i + 1] : 0);
			int z1 = y % M, z2 = y / M;
			b[i] = mul(t1[i][z1], t2[i][z2], p[i]);
		}
		return b[1].x;
	}
}

using sub1 :: calc;

int a[N]; int a0;

#define i0 i + i
#define i1 i + i + 1
ll t[N * 4]; int mi[N * 4];

void bt(int i, int x, int y) {
	if(x == y) {
		t[i] = v[x] % p[1];
		return;
	}
	int m = x + y >> 1;
	bt(i0, x, m); bt(i1, m + 1, y);
	t[i] = (t[i0] + t[i1]) % p[1];
}
int pl, pr, px;
void add(int i, int x, int y) {
	if(mi[i] >= p0 - 1 || y < pl || x > pr) return;
	if(x == y) {
		mi[i] ++;
		a0 = mi[i] + 1;
		fo(j, 1, a0 - 1) a[j] = c; a[a0] = v[x];
		t[i] = calc(a, a0);
		return;
	}
	int m = x + y >> 1;
	add(i0, x, m); add(i1, m + 1, y);
	t[i] = (t[i0] + t[i1]) % p[1];
	mi[i] = min(mi[i0], mi[i1]);
}
void ft(int i, int x, int y) {
	if(y < pl || x > pr) return;
	if(x >= pl && y <= pr) {
		px = (px + t[i]) % p[1];
		return;
	}
	int m = x + y >> 1;
	ft(i0, x, m); ft(i1, m + 1, y);
}

int op;

int main() {
	scanf("%d %d %d %d", &n, &m, &sp, &c);
	p[p0 = 1] = sp;
	while(p[p0] != 1) p0 ++, p[p0] = phi(p[p0 - 1]);
	p[++ p0] = 1;
	fo(i, 1, n) scanf("%d", &v[i]);
	sub1 :: build();
	bt(1, 1, n);
	fo(ii, 1, m) {
		scanf("%d %d %d", &op, &pl, &pr);
		if(op == 0) {
			add(1, 1, n);
		} else {
			px = 0;
			ft(1, 1, n);
			pp("%d
", px);
		}
	}
}