根据扩展欧拉定理
a
b
%
p
≡
a
b
%
ϕ
(
p
)
+
ϕ
(
p
)
%
p
,
b
≥
ϕ
(
p
)
a^b\% p ≡ a^{b \%\phi(p)+\phi(p)}\% p_,b≥\phi(p)
ab%p≡ab%ϕ(p)+ϕ(p)%p,b≥ϕ(p)
a
b
%
p
≡
a
b
%
ϕ
(
p
)
%
p
,
b
<
ϕ
(
p
)
a^b\% p ≡ a^{b \%\phi(p)}\% p,b<\phi(p)
ab%p≡ab%ϕ(p)%p,b<ϕ(p)
当
p
p
p 是偶数时,
p
p
p 变
ϕ
(
p
)
\phi(p)
ϕ(p) 至少除以
2
2
2;当
p
p
p 是奇数时,根据
ϕ
(
p
)
\phi(p)
ϕ(p) 的计算公式 (
q
[
i
]
q[i]
q[i] 为
p
p
p 的所有质因子,且互不相同)
p
∗
Π
i
=
1
k
(
q
[
i
]
−
1
)
÷
q
[
i
]
p*\Pi_{i=1}^k(q[i]-1)\div q[i]
p∗Πi=1k(q[i]−1)÷q[i]
q
[
i
]
−
1
q[i]-1
q[i]−1必有偶数,那么
ϕ
(
p
)
\phi(p)
ϕ(p)是偶数
可知
p
p
p 会在
O
(
log
p
)
O(\log p)
O(logp) 内次变为
1
1
1,结束递归
而上题中,对扩展欧拉定理的应用为:
2
2
2
2
.
.
.
%
p
≡
2
2
2
2
.
.
.
%
ϕ
(
p
)
+
ϕ
(
p
)
2^{2^{2^{2^{...}}}}\% p ≡ 2^{2^{2^{2^{...}}} \%\phi(p)+\phi(p)}
2222...%p≡2222...%ϕ(p)+ϕ(p)
假设
2
2
2 的个数有限,等式左边有
x
x
x 个
2
2
2, 那么等式右边的指数,也就是递归处理的部分,有
x
−
1
x-1
x−1 个
2
2
2
也就是说,每次递归减少了一个
2
2
2
那么减少
O
(
log
p
)
O(\log p)
O(logp) 个
2
2
2 后必定结束递归,即当
2
2
2 的个数超过
log
p
\log p
logp 时,答案全部一样
回到本题,可以得出一个数的修改次数上限为
O
(
log
p
)
O(\log p)
O(logp)
那么维护线段树,如果整个区间修改次数都到上限了,就
r
e
t
u
r
n
return
return
易得一共会碰到叶子节点
O
(
n
log
p
)
O (n \log p)
O(nlogp) 次
为了方(偷)便(懒) ,记
b
[
i
]
[
j
]
b[i][j]
b[i][j] 为
c
c
c
.
.
.
a
i
c^{c^{c^{...^{a_{i}}}}}
ccc...ai (共
j
j
j 个
c
c
c,注意没有取模),记
f
(
i
,
j
,
p
)
f(i,j,p)
f(i,j,p) 为
c
c
c
.
.
.
a
i
%
p
c^{c^{c^{...^{a_{i}}}}} \%p
ccc...ai%p (共
j
j
j 个
c
c
c),那么有递推式
f
(
i
,
j
,
p
)
=
c
f
(
i
,
j
−
1
,
ϕ
(
p
)
)
+
(
b
[
i
]
[
j
−
1
]
<
ϕ
(
p
)
?
0
:
ϕ
(
p
)
)
f(i,j,p)=c^{f(i,j-1,\phi(p))+(b[i][j-1]<\phi(p)?0:\phi(p))}
f(i,j,p)=cf(i,j−1,ϕ(p))+(b[i][j−1]<ϕ(p)?0:ϕ(p))
预处理几个值在
1
e
8
1e8
1e8 内的
b
[
i
]
[
j
]
b[i][j]
b[i][j] 就好了,显然
j
j
j 都比较小,
j
j
j 稍大一点就是天文数字了,肯定比
ϕ
(
p
)
\phi(p)
ϕ(p) 大,当然,注意特判
c
=
1
c=1
c=1 和
j
j
j 较大的情况
p
p
p 有
1
e
8
1e8
1e8,不能用数组把所有的
ϕ
(
x
)
,
x
≤
p
\phi(x),x≤p
ϕ(x),x≤p 都存下来,发现用到的
x
x
x 只有
p
,
ϕ
(
p
)
,
ϕ
(
ϕ
(
p
)
)
,
.
.
.
,
1
p,\phi(p),\phi(\phi(p)),...,1
p,ϕ(p),ϕ(ϕ(p)),...,1,考虑预处理出所有这样的
x
x
x,并存下它们的
ϕ
\phi
ϕ,记这些
x
x
x 分别为
d
[
m
.
.
1
]
d[m..1]
d[m..1](即
d
[
k
−
1
]
=
ϕ
(
d
[
k
]
)
d[k-1] = \phi(d[k])
d[k−1]=ϕ(d[k])),那么上式改写为
f
(
i
,
j
,
k
)
=
c
f
(
i
,
j
−
1
,
k
−
1
)
+
(
b
[
i
]
[
j
−
1
]
<
d
[
k
−
1
]
?
0
:
d
[
k
−
1
]
)
f(i,j,k)=c^{f(i,j-1,k-1)+(b[i][j-1]<d[k-1]?0:d[k-1])}
f(i,j,k)=cf(i,j−1,k−1)+(b[i][j−1]<d[k−1]?0:d[k−1])
注意递归边界:
j
=
0
j=0
j=0 或
k
=
1
k=1
k=1
由于用到快速幂,上式时间复杂度
O
(
log
2
p
)
O (\log^2p)
O(log2p)
发现
f
(
i
,
j
,
m
)
f(i,j,m)
f(i,j,m) 不能直接由
f
(
i
,
j
−
1
,
m
)
f(i,j-1,m)
f(i,j−1,m) 得来,也就是说,每次到叶子节点都要重新计算,那么总复杂度
O
(
n
log
n
log
2
p
)
O (n \log n \log^2p)
O(nlognlog2p),直接爆炸,考虑优化掉快速幂
对于这
m
m
m个模数,分别预处理
c
1
[
i
]
c1[i]
c1[i] 为
c
i
∗
20000
c^{i*20000}
ci∗20000,
c
2
[
i
]
c2[i]
c2[i] 为
c
i
c^i
ci,那么
c
y
=
c
1
[
y
/
20000
]
∗
c
2
[
y
%
20000
]
c^y=c1[y/20000]*c2[y\%20000]
cy=c1[y/20000]∗c2[y%20000]
然后就可以
O
(
1
)
O(1)
O(1) 快速幂了
总时间复杂度
O
(
n
(
log
n
log
p
+
log
2
p
)
)
O(n (\log n \log p + \log^2p))
O(n(lognlogp+log2p))
Code
#include<bits/stdc++.h>usingnamespace std;#define ll long long#define p2 p << 1#define p3 p << 1 | 1template<classt>inlinevoidread(t & res){char ch;while(ch =getchar(),!isdigit(ch));
res = ch ^48;while(ch =getchar(),isdigit(ch))
res = res *10+(ch ^48);}constint e =5e4+5, o =5e4+5, bl =2e4, z =105;int a[e], cnt[o *4], c1[z][20005], c2[z][20005], n, m, c, mod, phi[e];int num[e], d[e], c4, c3[e], b[e][20], sum[o *4];inlinevoidupt(int&x,int y){
x = y;if(x >= mod) x -= mod;}inlineintsolve(int x){int s =sqrt(x), res = x, i;for(i =2; i <= s; i++)if(x % i ==0){while(x % i ==0) x /= i;
res /= i;
res *= i -1;}if(x >1){
res /= x;
res *= x -1;}return res;}inlinevoidinit(){
d[d[0]=1]= mod;for(;;){
d[0]++;
d[d[0]]=solve(d[d[0]-1]);if(d[d[0]]==1)break;}reverse(d +1, d + d[0]+1);int i, j;
ll y =1;
c3[0]=1;for(i =1; i < bl; i++){if(y <=1e8){
y = y * c;
c3[i]= y;if(y <=1e8) c4 = i;}elsebreak;}for(i =1; i <= d[0]; i++){
c1[i][0]= c2[i][0]=1;for(j =1; j < bl; j++) c2[i][j]=(ll)c2[i][j -1]* c % d[i];
c1[i][1]=(ll)c2[i][bl -1]* c % d[i];for(j =2; j < bl; j++) c1[i][j]=(ll)c1[i][j -1]* c1[i][1]% d[i];}for(i =1; i <= n; i++){
b[i][0]= a[i];for(j =1; j <=6; j++){int x = b[i][j -1], p1 = x / bl, p0 = x % bl;if(c3[p0]>1e8|| p1 || p0 > c4)break;
b[i][j]= c3[p0];
num[i]= j;}}}inlineintksm(int y,int id){return(ll)c1[id][y / bl]* c2[id][y % bl]% d[id];}inlineintf(int i,int j,int mod){if(mod ==1)return0;if(j ==0)return a[i]% d[mod];if(c ==1||(c !=1&& j -1<= num[i]&& b[i][j -1]< d[mod -1]))returnksm(f(i, j -1, mod -1), mod);elsereturnksm(f(i, j -1, mod -1)+ d[mod -1], mod);}inlinevoidcollect(int p){
cnt[p]=min(cnt[p2], cnt[p3]);upt(sum[p], sum[p2]+ sum[p3]);}inlinevoidbuild(int l,int r,int p){if(l == r){
sum[p]= a[l];return;}int mid = l + r >>1;build(l, mid, p2);build(mid +1, r, p3);collect(p);}inlinevoidupdate(int l,int r,int s,int t,int p){if(cnt[p]> d[0])return;if(l == r){
cnt[p]++;
sum[p]=f(l, cnt[p], d[0]);return;}int mid = l + r >>1;if(s <= mid)update(l, mid, s, t, p2);if(t > mid)update(mid +1, r, s, t, p3);collect(p);}inlineintquery(int l,int r,int s,int t,int p){if(l == s && r == t)return sum[p];int mid = l + r >>1, res =0;if(t <= mid) res =query(l, mid, s, t, p2);elseif(s > mid) res =query(mid +1, r, s, t, p3);elseupt(res,query(l, mid, s, mid, p2)+query(mid +1, r, mid +1, t, p3));return res;}intmain(){read(n);read(m);read(mod);read(c);int i, opt, l, r;for(i =1; i <= n;++i)read(a[i]);init();build(1, n,1);while(m--){read(opt);read(l);read(r);if(!opt)update(1, n, l, r,1);elseprintf("%d\n",query(1, n, l, r,1));}return0;}