[BZOJ4869][Shoi2017]相逢是问候 扩展欧拉定理+势能线段树

根据扩展欧拉定理
当$x\geq \phi(p)$时有$a^x\equiv a^{x \% \phi(p)+\phi(p)}(mod p)$
我们发现在至多进行$t$次操作后，$a \equiv c^a(modp)$，t为进行$p=\phi(p)$多少次之后$p=1$，$t$的大小为$O(log_2p)$级别的
对此，每个数至多操作$log_2p$次就不变了，对于长度为$L$的区间，至多操作$L*log_2p$次修改
线段树维护区间势能，就是这个区间还能被修改多少次，不修改势能为0的节点
时间代价为$T(N)=2T(\frac{N}{2})+Nlog_2N$，$T(N)=Nlog_2^2N$

#include <bits/stdc++.h>
#define mid ((l+r)>>1)
#define ls l,mid,t<<1
#define rs mid+1,r,t<<1^1
#define ___ int l,int r,int t

#define N 50050
#define tp 50000

using namespace std;
typedef long long LL;
struct Node{int sum,tot;}tr[4*N];
int C[tp+5][25],CC[tp+5][25],a[N],vis[N],n,m,ll,rr,p,c,ans;
int F[N][25],phi[N],cnt;
inline int rd() {int r;scanf("%d",&r);return r;}
inline void inc(int &x,int y) {x=(x+y)%p;}

inline int qp(int a,int b,int p) {
    int ret = 1;
    while (b) {
        if (b&1) ret = 1LL * ret * a % p;
        b >>= 1, a = 1LL * a * a % p;
    }
    return ret;
}

Node operator+(Node p1, Node p2) {
    Node tmp;
    tmp.sum = (p1.sum + p2.sum) % p;
    tmp.tot = p1.tot + p2.tot;
    return tmp;
}

Node build(___) {
    return tr[t] = l==r ? (Node){a[l],vis[l]<cnt} : build(ls) + build(rs);
}

void update(___) {
    if (!tr[t].tot) return ;
    if (l==r) {
        vis[l]<cnt ? ++vis[l]:0;
        tr[t].sum = F[l][vis[l]] % p;
        tr[t].tot = (int)(vis[l]<cnt);
        return ;
    }
    if (ll>=l && r<=rr && !tr[t].tot) return ;
    if (ll<=mid) update(ls);
    if (rr>mid) update(rs);
    tr[t]=tr[t<<1]+tr[t<<1^1];
}

void query(___) {
    if (l>=ll&&r<=rr) {inc(ans, tr[t].sum);return ;}
    if (ll <= mid) query(ls);
    if (rr > mid) query(rs);
}

int get(int x) {
    int ret = x;;
    for (int i=2;i*i<=x;i++) if (x%i == 0) {
        ret = 1LL * ret * (i-1) / i;
        while (x%i == 0) x /= i;
    }
    if (x-1) return 1LL * ret * (x-1) / x;
    return ret;
}

inline int nqp(int t,int p) {
    if (t<=tp)
        return C[t][p];
    else
        return 1LL * CC[t/tp][p] * C[t%tp][p] % phi[p];
}

void pre() {

    int cur = p;
    while (cur != 1) {
        phi[cnt] = cur;
        ++cnt;
        cur = get(cur);
    }
    phi[cnt] = 1;
    phi[++cnt] = 1;

}

int main() {
    n = rd(), m = rd(), p = rd(), c = rd();
    for (int i=1;i<=n;i++) a[i] = rd(), vis[i] = 0;
    pre();

    for (int x=0;x<=cnt;x++) {
        C[0][x] = 1 % phi[x];
        for (int i=1;i<=tp;i++) C[i][x] = 1LL * c * C[i-1][x] % phi[x];
        int g = qp(c, tp, phi[x]);
        CC[0][x] = 1 % phi[x];
        for (int i=1;i<=tp;i++) CC[i][x] = 1LL * g * CC[i-1][x] % phi[x];
    }

    for (int i=1;i<=n;i++) {
        F[i][0] = a[i];
        for (int j=1;j<=cnt;j++) {
            int cur = a[i];
            if (cur >= phi[j]) cur = cur % phi[j] + phi[j];
            for (int k=j;k>=1;k--) {
                cur = nqp(cur, k-1);
                if (!cur) cur += phi[k-1];
            }
            F[i][j] = cur % phi[0];
        }
    }

    build(1,n,1);
    while (m--) {
        int cmd = 0;
        cmd = rd(), ll = rd(), rr = rd();
        if (cmd == 0) {
            update(1,n,1);
        } else {
            ans = 0;
            query(1,n,1);
            printf("%d\n",ans);
        }
    }
    return 0;
}