炒冷饭
求逆元
exgcd(a,m,x,y); ans=(x%m+m)%m;
,万能。qpow(a,m-2,m)
,这个仅限于m是一个质数。inv[i]=(m-m/i)*inv[m%i]%m
,线性求逆元。
中国剩余定理(CRT)
现在已知一个同余方程组:
{
x
≡
a
1
m
o
d
m
1
x
≡
a
2
m
o
d
m
2
.
.
.
x
≡
a
n
m
o
d
m
n
\begin{cases} x\equiv a_1\mod m_1\\ x\equiv a_2\mod m_2\\ ...\\ x\equiv a_n\mod m_n \end{cases}
⎩⎪⎪⎪⎨⎪⎪⎪⎧x≡a1modm1x≡a2modm2...x≡anmodmn
满足
m
i
m_i
mi两两互质。
解法:
首先计算出下面的各个值:
M
=
∏
i
=
1
m
m
i
P
i
=
M
m
i
Q
i
≡
P
i
−
1
m
o
d
m
i
M=\prod\limits_{i=1}^m m_i\\ P_i={M\over m_i}\\ Q_i \equiv P_i^{-1}\mod m_i
M=i=1∏mmiPi=miMQi≡Pi−1modmi
那么最终的答案就是:
a
n
s
≡
∑
i
=
1
n
a
i
P
i
Q
i
m
o
d
M
ans\equiv\sum\limits_{i=1}^n a_iP_iQ_i \mod M
ans≡i=1∑naiPiQimodM
我不会证明
扩展中国剩余定理(EXCRT)
现在解决的是 m i m_i mi可能不两两互质的情况。
符号约定:
i n v ( a , m ) \mathrm{inv}(a,m) inv(a,m)表示 a a a在模 m m m下的逆元。
( a , b ) = gcd ( a , b ) (a,b)=\gcd(a,b) (a,b)=gcd(a,b)
考虑将两条方程进行合并。
现在有两条方程:
{
x
≡
a
1
m
o
d
m
1
x
≡
a
2
m
o
d
m
2
\begin{cases} x\equiv a_1\mod m_1\\ x\equiv a_2\mod m_2 \end{cases}
{x≡a1modm1x≡a2modm2
转化为不定方程的形式:
{
x
=
a
1
+
m
1
y
1
x
=
a
2
+
m
2
y
2
⇒
{
y
1
=
x
−
a
1
m
1
y
1
=
x
−
a
2
m
2
\begin{cases} x=a_1+m_1 y_1\\ x=a_2+m_2 y_2 \end{cases} \Rightarrow \begin{cases} y_1=\frac{x-a_1}{m_1}\\ y_1=\frac{x-a_2}{m_2} \end{cases}
{x=a1+m1y1x=a2+m2y2⇒{y1=m1x−a1y1=m2x−a2
那么可以得到:
a
1
+
m
1
y
1
=
a
2
+
m
2
y
2
m
1
y
1
−
m
2
y
2
=
a
2
−
a
1
\begin{aligned} a_1+m_1 y_1&=a_2+m_2 y_2\\ m_1 y_1-m_2 y_2&=a_2-a_1\\ \end{aligned}\\
a1+m1y1m1y1−m2y2=a2+m2y2=a2−a1
可以发现,这个方程有解的条件是
(
m
1
,
m
2
)
∣
(
a
2
−
a
1
)
(m_1, m_2)|(a_2-a_1)
(m1,m2)∣(a2−a1) ,方程两边同时除以
(
m
1
,
m
2
)
(m_1, m_2)
(m1,m2),然后化为同余式:
m
1
(
m
1
,
m
2
)
y
1
≡
a
2
−
a
1
(
m
1
,
m
2
)
m
o
d
m
2
(
m
1
,
m
2
)
\frac{m_1}{(m_1,m_2)}y_1\equiv\frac{a_2-a_1}{(m_1, m_2)}\mod{\frac{m_2}{(m_1,m_2)}}
(m1,m2)m1y1≡(m1,m2)a2−a1mod(m1,m2)m2
可以发现此时
(
m
1
(
m
1
,
m
2
)
,
m
2
(
m
1
,
m
2
)
)
=
1
\left(\frac{m_1}{(m_1, m_2)}, \frac{m_2}{(m_1,m_2)}\right)=1
((m1,m2)m1,(m1,m2)m2)=1,所以同余式两边同时乘上逆元:
⇒
y
1
≡
a
2
−
a
1
(
m
1
,
m
2
)
×
i
n
v
(
m
1
(
m
1
,
m
2
)
,
m
2
(
m
1
,
m
2
)
)
m
o
d
m
2
(
m
1
,
m
2
)
\Rightarrow y_1\equiv \frac{a_2-a_1}{(m_1, m_2)}\times \mathrm{inv}\left(\frac{m_1}{(m_1, m_2)},\frac{m_2}{(m_1,m_2)}\right) \mod \frac{m_2}{(m_1,m_2)}
⇒y1≡(m1,m2)a2−a1×inv((m1,m2)m1,(m1,m2)m2)mod(m1,m2)m2
接着将其再次转化为一个不定方程,然后带入
x
x
x:
⇒
y
1
=
a
2
−
a
1
(
m
1
,
m
2
)
×
i
n
v
(
m
1
(
m
1
,
m
2
)
,
m
2
(
m
1
,
m
2
)
)
%
m
2
(
m
1
,
m
2
)
+
k
×
m
2
(
m
1
,
m
2
)
⇒
x
−
a
1
m
1
=
a
2
−
a
1
(
m
1
,
m
2
)
×
i
n
v
(
m
1
(
m
1
,
m
2
)
,
m
2
(
m
1
,
m
2
)
)
%
m
2
(
m
1
,
m
2
)
+
k
×
m
2
(
m
1
,
m
2
)
⇒
x
=
a
2
−
a
1
(
m
1
,
m
2
)
×
i
n
v
(
m
1
(
m
1
,
m
2
)
,
m
2
(
m
1
,
m
2
)
)
%
m
2
(
m
1
,
m
2
)
×
m
1
+
a
1
+
k
×
m
1
m
2
(
m
1
,
m
2
)
\Rightarrow y_1= \frac{a_2-a_1}{(m_1, m_2)} \times \mathrm{inv}\left(\frac{m_1} {(m_1,m_2)},\frac{m_2}{(m_1,m_2)}\right) \% \frac{m_2}{(m_1,m_2)} + k \times \frac{m_2}{(m_1,m_2)}\\ \Rightarrow \frac{x-a_1}{m_1}= \frac{a_2-a_1}{(m_1, m_2)} \times \mathrm{inv}\left(\frac{m_1}{(m_1, m_2)},\frac{m_2}{(m_1,m_2)}\right) \% \frac{m_2}{(m_1,m_2)} + k \times \frac{m_2}{(m_1,m_2)}\\ \Rightarrow x= \frac{a_2-a_1} {(m_1, m_2)} \times \mathrm{inv}\left(\frac{m_1}{(m_1, m_2)},\frac{m_2}{(m_1,m_2)}\right) \% \frac{m_2}{(m_1,m_2)} \times m_1 + a_1 + k \times \frac{m_1 m_2}{(m_1,m_2)}
⇒y1=(m1,m2)a2−a1×inv((m1,m2)m1,(m1,m2)m2)%(m1,m2)m2+k×(m1,m2)m2⇒m1x−a1=(m1,m2)a2−a1×inv((m1,m2)m1,(m1,m2)m2)%(m1,m2)m2+k×(m1,m2)m2⇒x=(m1,m2)a2−a1×inv((m1,m2)m1,(m1,m2)m2)%(m1,m2)m2×m1+a1+k×(m1,m2)m1m2
然后转化为同余式:
⇒
x
≡
a
2
−
a
1
(
m
1
,
m
2
)
×
i
n
v
(
m
1
(
m
1
,
m
2
)
,
m
2
(
m
1
,
m
2
)
)
%
m
2
(
m
1
,
m
2
)
×
m
1
+
a
1
(
m
o
d
m
1
m
2
(
m
1
,
m
2
)
)
\Rightarrow x\equiv \frac{a_2-a_1} {(m_1, m_2)} \times \mathrm{inv}\left(\frac{m_1}{(m_1, m_2)},\frac{m_2}{(m_1,m_2)}\right) \% \frac{m_2}{(m_1,m_2)} \times m_1 + a_1 ~~~\left(\mathrm{mod} \frac{m_1 m_2}{(m_1,m_2)}\right)
⇒x≡(m1,m2)a2−a1×inv((m1,m2)m1,(m1,m2)m2)%(m1,m2)m2×m1+a1 (mod(m1,m2)m1m2)
这样就合并成功了。
同余式的转化
在实际问题中,可能给出的同余式是类似下面这个样子的:
a
x
≡
b
m
o
d
M
ax\equiv b\mod M
ax≡bmodM
我们需要将其转化为(EX)CRT的形式的形式。
将其转化为不定方程的形式:
a
x
+
M
y
=
b
ax+My=b
ax+My=b
可以发现,这个方程有解的条件是
(
a
,
M
)
∣
b
(a,M)|b
(a,M)∣b。
首先使用exgcd求出方程
a
x
+
M
y
=
(
a
,
M
)
ax+My=(a,M)
ax+My=(a,M)的特解,得到的一组解为
x
′
,
y
′
x',y'
x′,y′。
那么原方程的一组解为
S
x
=
b
(
a
,
M
)
x
′
S
y
=
b
(
a
,
M
)
y
′
Sx=\frac{b}{(a,M)}x'\\ Sy=\frac{b}{(a,M)}y'
Sx=(a,M)bx′Sy=(a,M)by′
那么原方程的通解就是
x
=
S
x
+
k
×
b
(
a
,
M
)
y
=
S
y
−
k
×
b
(
a
,
M
)
x=Sx+k\times \frac{b}{(a,M)}\\ y=Sy-k\times \frac{b}{(a,M)}
x=Sx+k×(a,M)by=Sy−k×(a,M)b
我们取x的通解表达式,将其转化为同余式:
x
≡
S
x
(
m
o
d
b
(
a
,
M
)
)
x\equiv Sx~~~~\left(\mathrm{mod}~~\frac{b}{(a,M)}\right)
x≡Sx (mod (a,M)b)
新内容
基本概念
阶:当 ( a , p ) = 1 (a,p)=1 (a,p)=1 时,满足 a x ≡ 1 ( m o d p ) a^x\equiv 1~~~~(\mathrm{mod}~~p) ax≡1 (mod p) 的最小 x x x ,计作 δ p ( a ) \delta_p(a) δp(a)。
原根:如果 δ p ( a ) = φ ( p ) \delta_p(a)=\varphi(p) δp(a)=φ(p) ,那么 a a a 就是模 p p p 下的原根。一个模数由原根,就必须先满足 p = 2 , 4 , p n , 2 p n p=2,4,p^n,2p^n p=2,4,pn,2pn 。
如果 p p p 是一个质数,那么可以用 a x % p a^x\%p ax%p 表示 [ 1 , p − 1 ] [1,p-1] [1,p−1] 中的所有数,其中 1 ≤ x ≤ p − 1 1\leq x\leq p-1 1≤x≤p−1
二次剩余
形式化的,就是找到满足下面的同余式的
x
x
x:
x
2
≡
a
(
m
o
d
p
)
x^2\equiv a~~~~(\mathrm{mod}~~p)
x2≡a (mod p)
通俗点来讲,就是将
a
a
a在模
p
p
p以一下进行开根。
如果 a a a 能够开根(上面的同余方程有解),那么我们就可以认为 “ a a a 是模 p p p 意义下的二次剩余”,否则认为“ a a a 不是模 p p p 意义下的二次剩余”。
上面同余方程要么有两个解(可能相同),要么没有解。
如果 p = 2 p=2 p=2 ,那么显然的, x x x 的取值只有 0 / 1 0/1 0/1 ,而且很容易找到,所以我们只讨论 p p p 属于奇素数的情况。
引理1:对于一个模数 p p p ,是二次剩余的数共有 p − 1 2 \frac{p-1}{2} 2p−1(不包括 0 0 0 ),非二次剩余的数也有 p − 1 2 \frac{p-1}{2} 2p−1 个,各占一半。
虽然我们很难在较短时间内准确知道是二次剩余的数的分布,但是没有关系。
引理2:对于一个数
a
a
a和一个模数
p
p
p,可以得到:
v
≡
a
p
−
1
2
(
m
o
d
p
)
{
v
=
1
⇔
a
是
模
p
的
二
次
剩
余
v
=
−
1
⇔
a
不
是
模
p
的
二
次
剩
余
v\equiv a^{\frac{p-1}{2}}~~~~(\mathrm{mod}~~p)\\ \begin{cases} v=1&\Leftrightarrow a是模p的二次剩余\\ v=-1&\Leftrightarrow a不是模p的二次剩余 \end{cases}
v≡a2p−1 (mod p){v=1v=−1⇔a是模p的二次剩余⇔a不是模p的二次剩余
接下来考虑算法(忽略各种特判)。
首先可以想到,第一步一定是判断一个数
a
a
a是否是二次剩余。
如果是一个二次剩余,那么继续执行算法,找到
x
x
x,否则直接返回一个表示不是二次剩余的结果。
接下来就是一系列神仙操作。
随机一个数
a
a
a,范围是
[
1
,
p
−
1
]
[1,p-1]
[1,p−1],满足
a
2
−
n
a^2-n
a2−n不是模
p
p
p的二次剩余。
接下来使用类似虚数这种拓展数域的方式,令
ω
=
a
2
−
n
\omega=\sqrt{a^2-n}
ω=a2−n,这个东西在模
p
p
p的数域下是不存在的。用二元组
(
a
,
b
)
(a,b)
(a,b)表示
a
+
b
ω
a+b\omega
a+bω,那么可以得到:
(
a
+
ω
)
p
−
1
=
(
a
+
ω
)
p
(
a
+
ω
)
=
(
∑
i
=
0
p
C
(
p
,
i
)
a
i
ω
p
−
i
)
(
a
+
ω
)
=
(
∑
i
=
1
p
C
(
⌊
p
/
p
⌋
,
⌊
i
/
p
⌋
)
C
(
p
%
p
,
i
%
p
)
a
i
ω
p
−
i
)
(
a
+
ω
)
\begin{aligned} (a+\omega)^{p-1}&=(a+\omega)^p(a+\omega)\\ &=\left(\sum\limits_{i=0}^p C(p,i)a^i\omega^{p-i}\right)(a+\omega)\\ &=\left(\sum\limits_{i=1}^p C(\lfloor{p/p}\rfloor,\lfloor{i/p}\rfloor) C(p\%p,i\%p)a^i\omega^{p-i}\right)(a+\omega)\\ \end{aligned}
(a+ω)p−1=(a+ω)p(a+ω)=(i=0∑pC(p,i)aiωp−i)(a+ω)=(i=1∑pC(⌊p/p⌋,⌊i/p⌋)C(p%p,i%p)aiωp−i)(a+ω)
很明显,
C
(
⌊
p
/
p
⌋
,
⌊
i
/
p
⌋
)
C(\lfloor{p/p}\rfloor,\lfloor{i/p}\rfloor)
C(⌊p/p⌋,⌊i/p⌋)的值永远为
1
1
1,而
C
(
p
%
p
,
i
%
p
)
C(p\%p,i\%p)
C(p%p,i%p)只有在
i
=
0
或
p
i=0或p
i=0或p的时候才会是
1
1
1,其他时刻都为
0
0
0,那么式子可以继续化简:
=
(
a
p
+
ω
p
)
(
a
+
ω
)
=
(
a
p
−
1
a
+
ω
p
−
1
ω
)
(
a
+
ω
)
=
(
a
p
−
1
a
+
(
a
2
+
n
)
p
−
1
2
ω
)
(
a
+
ω
)
\begin{aligned} &=(a^p+\omega^p)(a+\omega)\\ &=\left(a^{p-1}a+\omega^{p-1}\omega\right)(a+\omega)\\ &=\left(a^{p-1}a+\left(a^2+n\right)^{\frac{p-1}{2}}\omega\right)(a+\omega) \end{aligned}
=(ap+ωp)(a+ω)=(ap−1a+ωp−1ω)(a+ω)=(ap−1a+(a2+n)2p−1ω)(a+ω)
因为
a
p
−
1
≡
1
a^{p-1}\equiv 1
ap−1≡1,
(
a
2
+
n
)
p
−
1
2
≡
−
1
(a^2+n)^{\frac{p-1}{2}}\equiv -1
(a2+n)2p−1≡−1,所以
=
(
a
−
ω
)
(
a
+
ω
)
=
a
2
−
ω
2
=
n
\begin{aligned} &=(a-\omega)(a+\omega)\\ &=a^2-\omega^2\\ &=n \end{aligned}
=(a−ω)(a+ω)=a2−ω2=n
所以
(
a
+
ω
)
p
+
1
2
(a+\omega)^{\frac{p+1}{2}}
(a+ω)2p+1就是其中一个
x
x
x,用
p
−
x
p-x
p−x就能得到另一个
x
x
x。
struct Complex {
Complex(LL x, LL y) : x(x), y(y) {}
LL x, y;
};
inline Complex mul(Complex a, Complex b) {
return Complex(a.x*b.x%M + a.y*b.y%M*w%M, a.y*b.x*M+a.x*b.y%M);
}
LL qpow(LL a, LL b, LL M) {
LL ret = 1;
while (b) {
if (b & 1) ret = ret * a % M;
a = a * a % M, b >>= 1;
}
return ret;
}
Complex qpow(Complex a, LL b) {
Complex ret = Complex(1, 0);
while (b) {
if (b & 1) ret = mul(ret, a);
a = mul(a, a), b >>= 1;
}
return ret;
}
pair<LL, LL> solve(LL n, LL M) {
n %= M;
if (!n) return make_pair(1, 1);
if (M == 2) return make_pair(1, 1);
if (qpow(n, (M - 1) >> 1, M) == M - 1) return make_pair(-1,-1);
LL ans, ans1, ans2;
a = rand() % M;
while (!a || qpow((a * a - n + M) % M, (M - 1) >> 1, M)) a = rand() % M;
w = (a * a - n + M) % M;
ans1 = qpow(Complex(a, 1), (M + 1) >> 1).x;
ans2 = M - ans1;
if (ans1 > ans2) swap(ans1, ans2);
return make_pair(ans1, ans2);
}
大步小步(BSGS)
就是在模 p p p的情况下求出 l o g a b log_ab logab,非常暴力。
现在考虑的是 a , p a,p a,p互质的特殊情况,
首先使用Hash
表或者std::map
存储
b
×
a
i
(
0
≤
i
≤
⌈
b
⌉
)
b\times a^{i}(0\leq i\leq \lceil\sqrt{b}\rceil)
b×ai(0≤i≤⌈b⌉),以及对应的
i
i
i,然后枚举
a
i
⌈
p
⌉
(
0
≤
i
⌈
p
⌉
≤
p
)
a^{i\lceil{\sqrt{p}}\rceil}(0\leq i\lceil{\sqrt{p}}\rceil\leq p)
ai⌈p⌉(0≤i⌈p⌉≤p),如果得到一个值存在于表中,那么直接将指数相减就好了。
namespace BSGS {
map<LL, LL> val_map;
LL log(LL a, LL b, LL M) { // a^x = b (mod M)
val_map.clear();
LL sqrt_m = ceil(sqrt(M));
b %= M;
for (int i = 1; i <= sqrt_m; i++) {
b = b * a % M;
val_map[b] = i;
}
LL vec = qpow(a, sqrt_m, M); b = 1;
for (int i = 1; i <= sqrt_m; i++) {
b = b * vec % M;
if (val_map.count(b)) return (i * sqrt_m - val_map[b] + M) % M;
}
return -1;
}
}
拓展大步小步(EX_BSGS)
现在考虑 a , p a,p a,p不互质的情况。
设
g
=
(
a
,
p
)
a
x
≡
b
(
m
o
d
p
)
(
a
g
)
a
x
−
1
≡
b
g
(
m
o
d
p
g
)
设g=(a,p)\\ a^x\equiv b~~~~(\mathrm{mod}~~p)\\ \left(\frac{a}{g}\right)a^{x-1}\equiv \frac{b}{g}~~~~\left(\mathrm{mod}~~\frac{p}{g}\right)
设g=(a,p)ax≡b (mod p)(ga)ax−1≡gb (mod gp)
如果
b
b
b不能被
g
g
g整除,那么式子无解,因为
a
g
,
p
g
\frac{a}{g},\frac{p}{g}
ga,gp互质,所以一定有
a
g
\frac{a}{g}
ga在模
p
g
\frac{p}{g}
gp下的逆元,对式子进行变化:
a
x
−
1
≡
b
g
(
a
g
)
−
1
(
m
o
d
p
g
)
a^{x-1}\equiv \frac{b}{g}\left(\frac{a}{g}\right)^{-1}~~~~\left(\mathrm{mod}~~\frac{p}{g}\right)
ax−1≡gb(ga)−1 (mod gp)
这样子就得到另一个形式相同的式子,但是
a
a
a和
p
g
\frac{p}{g}
gp可能依然不互质,需要继续变化。如果最终互质了,就用BSGS进行求解。
namespace EX_BSGS {
LL log(LL a, LL b, LL M) {
if (b == 1 || M == 1) return 0;
LL g = gcd(a, M), cnt = 0, k = 1;
while (g > 1) {
if (b % g) return -1;
cnt++, b /= g, M /= g, k = k * (a / g) % M;
if (k == b) return cnt;
g = gcd(a, M);
}
LL ans2 = BSGS::log(a, b * inv(k, M) % M, M);
if (ans2 == -1) return -1;
else return ans2 + cnt;
}
}