Description：

给一个序列 $a_i$ ，支持将第 $ l $ 个到第 $ r $ 个数 $a_l,a_{l+1},\dots,a_r$ 中的每一个数 $a_i$ 替换为 $c^{a_i}$ ，求第 $ l $ 个到第 $ r $ 个数的和 $\%p$ 的值。

$1 \leqslant n \leqslant 50000, 1 \leqslant m \leqslant 50000, 1 \leqslant p \leqslant 100000000$

Solution：

首先有一个叫做欧拉定理的东西： $$ a^{\varphi(p)}\equiv1(\text{mod }p) $$ 在 $\gcd(a,p)=1$ 的时候成立。

然后有一个叫做扩展欧拉定理的东西： $$ a^b\equiv a^{b\%\varphi(p)}\ \gcd(a,p)=1 $$

$$ a^b\equiv a^b\ b<\varphi(p) $$

$$ a^b\equiv a^{b\%\varphi(p)+\varphi(p)}\ \gcd(a,p)\ne 1 $$

均在 $\%p$ 意义下。

然后观察这道题，要求的其实就是： $$ \sum_{i=l}^rc^{c^{\cdot^{\cdot^{\cdot^{c^{a_i}}}}}} $$ 对于每个数 $c^{c^{a_i}}$ 用扩展欧拉定理展开一下就是 $c^{(c^{a_i\%\varphi(\varphi(p))+\varphi(\varphi(p))})\%\varphi(p)+\varphi(p)}\%p$ ，然后对于一个数，不停对他取 $\varphi$ 至多 $O(\log n)$ 次就会变成 $1$ ，原因是对偶数取 $\varphi$ 一定会除二，因为所有偶数都和他不互质，对于所有奇数下一次一定会变成偶数，因为设 $n=\prod p_i^{q_i}$ ，那么对于一个奇数所有的 $p_i$ 一定都是奇数，而 $\varphi(n)=n\prod \frac{p_i-1}{p_i}$ ， $p_i-1$ 必为偶数，所以一定不超过 $O(\log n)$ 次就变成 $1$ ，这说明对一个数 $a_i$ 不断用 $c^{a_i}$ 替换他至多 $O(\log n)$ 次任何 $a_i$ 最后得到的答案都一样了，因为 $a_i$ 还没来得及影响答案 $p$ 就已经被削成 $1$ 了，那么就联想到一个区间开根的题，这个题的区间操作也可以类似的用一棵线段树来维护，如果整个区间都不会再变了，那就打一个标记，如果不是，就递归下去暴力计算，每个点的计算次数是 $O(\log n)$ 的，每次计算快速幂的代价是 $O(\log n)$ 的，线段树的复杂度是 $O(n\log n)$ 的，总代价为 $O(n\log^3n)$ 的，既然底数相同，可以预处理 $c$ 的幂，但不能预处理 $p=10^8$ 个幂，可以预处理 $c$ 的 $0$ 到 $9999$ 次方分别对 $\varphi(p)$ ， $\varphi(\varphi(p))$ 等等取模的值，再预处理 $c^{10000}$ 的 $0$ 到 $10000$ 次方分别对 $\varphi(p)$ ， $\varphi(\varphi(p))$ 等等取模的值，查询时分两半计算就可以 $O(1)$ 了，于是复杂度少了一个 $\log$ 。

关于如何暴力计算某个数的 $c$ 的多少个某个次方的值，观察扩展欧拉定理， $\gcd(a,p)=1$ 的情况最优美，可以并入后面两个，也就是说第一个是后两个的特殊情况，那么我们只要区分 $2$ 和 $3$ 情况就行了，这两种情况只有质数和 $\varphi(p)$ 大小关系的差距。关于怎么比较这个，网上的题解都是用取模之后的值来判断，不知道这样是否合理

Code：

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<cstring>

using namespace std;

inline int read()

{

    register int res = 0,f = 1;

    register char c = getchar();

    while(!isdigit(c))

    {

        if(c == '-')f = -1;

        c = getchar();

    }

    while(isdigit(c))

    {

        res = (res << 1) + (res << 3) + c - '0';

        c = getchar();

    }

    return res;

}

int n,m,p,c;

#define MAXN 50010

int a[MAXN];

int phi(int x)

{

    int res = x,s = x;

    for(int i = 2;i * i <= x;++i)

    {

        if(s % i == 0)

        {

            res = res / i * (i - 1);

            while(s % i == 0)s /= i;

        }

    }

    if(s != 1)res = res / s * (s - 1);

    return res;

}

#define MAX 10001

#define MAXPHI_TIMES 100

int phi_times;

int phi_val[MAXPHI_TIMES];

int p1[MAXPHI_TIMES][MAX];

int p2[MAXPHI_TIMES][MAX];

int f[MAXN][MAXPHI_TIMES];

int power_c(int times,int mod)

{

    return 1ll * p2[mod][times / 10000] * p1[mod][times % 10000] % phi_val[mod];

}

int calc(int base,int cur,int times,int top)

{

    if(cur == times)return top;

    int k = calc(base,cur + 1,times,top);

    k = power_c(k,cur);

    if(k == 0)k = phi_val[cur];

    return k;                                      

}

void pre_work()

{

    phi_val[0] = p;

    for(phi_times = 1;;++phi_times)

    {

        phi_val[phi_times] = phi(phi_val[phi_times - 1]);

        if(phi_val[phi_times - 1] == 1)break;

    }

    for(int i = 0;i <= phi_times;++i)

    {

        p1[i][0] = 1;

        for(int j = 1;j < MAX;++j)

        {

            p1[i][j] = 1ll * p1[i][j - 1] * c % phi_val[i];

        }

    }

    for(int i = 0;i <= phi_times;++i)

    {

        int base = p1[i][10000];

        p2[i][0] = 1;

        for(int j = 1;j < MAX;++j)

        {

            p2[i][j] = 1ll * p2[i][j - 1] * base % phi_val[i];

        }

    }

    for(int i = 1;i <= n;++i)

        for(int k = 0;k <= phi_times;++k)

            f[i][k] = calc(c,0,k,a[i]);

    return;

}

struct node

{

    int lc,rc;

    int sum,tag;

    node(){sum = tag = 0;}

}t[MAXN << 1];

int ptr = 0;

int newnode(){return ++ptr;}

int root;

#define mid ((l + r) >> 1)

void build(int &rt,int l,int r)

{

    rt = newnode();

    if(l == r)

    {

        t[rt].sum = a[l] % p;t[rt].tag = 0;

        return;

    }

    build(t[rt].lc,l,mid);

    build(t[rt].rc,mid + 1,r);

    t[rt].sum = (t[t[rt].lc].sum + t[t[rt].rc].sum) % p;

    t[rt].tag = min(t[t[rt].lc].tag,t[t[rt].rc].tag);

    return;

}

void change(int rt,int L,int R,int l,int r)

{

    if(t[rt].tag >= phi_times)return;

    if(l == r)

    {

        t[rt].sum = calc(c,0,++t[rt].tag,a[l]);

        return;

    }

    if(L <= mid)change(t[rt].lc,L,R,l,mid);

    if(R > mid)change(t[rt].rc,L,R,mid + 1,r);

    t[rt].sum = (t[t[rt].lc].sum + t[t[rt].rc].sum) % p;

    t[rt].tag = min(t[t[rt].lc].tag,t[t[rt].rc].tag);

    return;

}

int query(int rt,int L,int R,int l,int r)

{

    if(L <= l && r <= R)return t[rt].sum;

    int res = 0;

    if(L <= mid)res += query(t[rt].lc,L,R,l,mid);

    if(R > mid)res += query(t[rt].rc,L,R,mid + 1,r);

    return res % p;

}

void pt(int rt,int l,int r)

{

    if(l == r)

    {

        cout << t[rt].sum << " ";if(r == n)cout << endl;

        return;

    }

    pt(t[rt].lc,l,mid);

    pt(t[rt].rc,mid + 1,r);

    return;

}

int main()

{

    scanf("%d%d%d%d",&n,&m,&p,&c);

    for(int i = 1;i <= n;++i)scanf("%d",&a[i]);

    pre_work();

    build(root,1,n);

    int opt,l,r;

    for(int i = 1;i <= m;++i)

    {//pt(root,1,n);

        opt = read();l = read();r = read();

        if(opt == 0)

        {

            change(root,l,r,1,n);

        }

        else

        {

            printf("%d\n",query(root,l,r,1,n));

        }

    }

    return 0;

}