Luogu 3934 Nephren Ruq Insania

和Ynoi2016 炸脖龙重题了。

BZOJ 5394。

首先是扩展欧拉定理：

一开始傻掉了……递归的层数和区间长度无关……也就是说我们每一次直接暴力递归求解子问题一定不会超过$logP$层，因为当模数变成$1$了的时候，不管怎么乘方都不会再变化了。

然后要注意$b < phi (P)$的情况，在快速幂的时候判断一下，如果超过了$phi (P)$，那么就再加上一个$P$。

对于区间修改，树状数组或者线段树均可，每一层直接查询。

一共有不超过$logP$层，每一层有$logn$的开销查询元素，所以最后的复杂度是$O(Maxn + qlogP * logn)$。

另外，BZOJ上要把$pri$开成一半才不会爆空间……

Code：

#include <cstdio>
#include <cstring>
using namespace std;
typedef long long ll;

const int N = 5e5 + 5;
const int M = 2e7 + 5;
const int Maxn = 2e7;

int n, qn, pCnt = 0, pri[M / 2];
ll a[N], phi[M];
bool np[M];

template <typename T>
inline void read(T &X) {
    X = 0; char ch = 0; T op = 1;
    for(; ch > '9' || ch < '0'; ch = getchar())
        if(ch == '-') op = -1;
    for(; ch >= '0' && ch <= '9'; ch = getchar())
        X = (X << 3) + (X << 1) + ch - 48;
    X *= op;
}    

inline void sieve() {
    phi[1] = 1;
    for(int i = 2; i <= Maxn; i++) {
        if(!np[i]) pri[++pCnt] = i, phi[i] = i - 1;
        for(int j = 1; j <= pCnt && pri[j] * i <= Maxn; j++) {
            np[i * pri[j]] = 1;
            if(i % pri[j] == 0) {
                phi[i * pri[j]] = phi[i] * pri[j];
                break;
            }
            phi[i * pri[j]] = phi[i] * (pri[j] - 1);
        }
    }
}

namespace Bit {
    ll s[N];
    
    #define lowbit(p) (p & (-p))
    
    inline void modify(int p, ll v) {
        for(; p <= n; p += lowbit(p))
            s[p] += v;
    }
    
    inline ll query(int p) {
        ll res = 0;
        for(; p > 0; p -= lowbit(p))
            res += s[p];
        return res;
    }
    
} using namespace Bit;

inline ll fpow(ll x, ll y, ll P) {
    ll res = 1; bool tag = 0, tagx = 0;
    for(; y > 0; y >>= 1) {
        if(y & 1) {
            res = res * x;
            tag |= tagx;
            if(res >= P) res %= P, tag = 1;
        }
        if(x >= P) x %= P, tagx = 1;
        x = x * x;
        if(x >= P) x %= P, tagx = 1;
    }
    return res + (tag ? P : 0);
}

ll solve(int x, int y, ll P) {
    if(P == 1) return 1;
    if(x > y) return 1;
    ll val = a[x] + query(x), cur = solve(x + 1, y, phi[P]);
    return fpow(val, cur, P);
}

int main() {
    sieve();
    read(n), read(qn);
    for(int i = 1; i <= n; i++) read(a[i]);
    for(int op, x, y; qn--; ) {
        read(op), read(x), read(y);
        if(op == 1) {
            ll v; read(v);
            modify(x, v), modify(y + 1, -v);
        } else {
            ll P; read(P);
            printf("%lld
", solve(x, y, P) % P);
        }
    }
    return 0;
}