逆元与快速幂学习
逆元与快速幂
逆元的数学原理与其在快速幂中的应用 原题为蓝桥杯国赛 小数的第n位
总结算法的学习过程：这是一个什么样的原理，我如何将其他问题转化为这样的算法
逆元
单位元
【定义1】在一个集合中，对于某种运算 ∗ * ∗（注意：这里代表通用运算的表示符号，并不是特指乘法），如果对于任何的集合元素 a a a ，和元素 e e e运算，得到还是集合元素 a a a本身，则称 e e e 为这个运算下的单位元
加法中的单位元为 a + e = a , e = 0 a + e = a, e = 0 a+e=a,e=0
乘法中的单位元为 a × e = a , e = 1 a \times e = a, e = 1 a×e=a,e=1
逆元
【定义2】在一个集合中，对于某种运算 ∗ * ∗ ，如果任意两个元素的运算结果等于单位元，则称这两个元素互为逆元。
在mod的意义下，不能直接除以一个数，而要除以它的逆元
Example：
a
∗
b
≡
1
(
m
o
d
p
)
⇔
x
/
b
=
x
∗
b
%
p
a * b \equiv 1(mod\ p)\\ \Leftrightarrow x / b = x * b \% p
a∗b≡1(mod p)⇔x/b=x∗b%p
加法运算中逆元 a + a − 1 = e = 0 , a − 1 = − a a + a^{-1} = e = 0, a^{-1} = -a a+a−1=e=0,a−1=−a
乘法运算中逆元 a × a − 1 = e = 1 , a − 1 = 1 a a \times a^{-1} = e = 1, a^{-1} = \frac{1}{a} a×a−1=e=1,a−1=a1
模乘中的单位元
对于 n n n的乘法，所有模 n n n和 a a a同余的数都可以表示成：
a ( m o d n ) = k n + a ( k ∈ Z ) a(mod\ n) = kn + a (k \in \Z) a(mod n)=kn+a(k∈Z)
另单位元为 e ( m o d n ) e(mod\ n) e(mod n)则根据单位元的概念可得
a ( m o d n ) × e ( m o d n ) = a ( m o d n ) ⇒ ( k 1 k 2 n + k 1 e + k 2 ) n + e a = k n + a ⇒ { k = k 1 k 2 n + k 1 e + k 2 e = 1 a(mod\ n) \times e(mod\ n) = a(mod\ n) \\ \Rightarrow (k_1 k_2 n + k_1 e + k_2)n + ea = kn + a\\ \Rightarrow \begin{cases} k = k_1 k_2 n + k_1 e + k_2 \\ e = 1 \end{cases} a(mod n)×e(mod n)=a(mod n)⇒(k1k2n+k1e+k2)n+ea=kn+a⇒{k=k1k2n+k1e+k2e=1
模乘中的逆元
因为对于任意模 n n n和 a a a同余的数的单位元都是 1 ( m o d n ) 1(mod\ n) 1(mod n)
要找 a ( m o d n ) a(mod\ n) a(mod n)的逆元 a − 1 ( m o d n ) a^{-1}(mod\ n) a−1(mod n)就要求两个元素满足
a ( m o d n ) × a − 1 ( m o d n ) = 1 ( m o d n ) a ( m o d n ) = k 1 n + a a − 1 ( m o d n ) = k 2 n + a − 1 a(mod\ n) \times a^{-1}(mod\ n) = 1(mod\ n)\\ a(mod\ n) = k_1 n + a \quad a^{-1}(mod\ n) = k_2 n + a^{-1} a(mod n)×a−1(mod n)=1(mod n)a(mod n)=k1n+aa−1(mod n)=k2n+a−1
也就是说要满足
( k 1 k 2 n + k 1 a − 1 + k 2 a ) n + a a − 1 = k 3 n + 1 ⇒ { k 1 k 2 n + k 1 a − 1 + k 2 a − k 3 = 0 a a − 1 − 1 = 0 Consequently: a a − 1 + k ⋅ n = 1 (k_1k_2n + k_1a^{-1} + k_2a)n + aa^{-1} = k_3n + 1 \\ \Rightarrow \begin{cases} k_1k_2n + k_1a^{-1} + k_2a -k_3 = 0 \\ aa^{-1} - 1= 0 \end{cases} \\ \text{Consequently:\quad} aa^{-1} + k \cdot n = 1 (k1k2n+k1a−1+k2a)n+aa−1=k3n+1⇒{k1k2n+k1a−1+k2a−k3=0aa−1−1=0Consequently:aa−1+k⋅n=1
其中 a , n a,n a,n已知，那么可以将求模的逆的问题转化成：
对于给定的正整数 a a a， b b b，求解下列方程中 x x x的最小整数解：（ x x x即为 a a a的逆元）
a x + b y = 1 ax + by = 1 ax+by=1
即欧几里德定理
思想： 整数域上的运算与一元多项式的运算是同构的，于是考虑的时候可以把一元多项式的相关结论拿出来用
求解过程：
(
1
)
a
x
+
b
y
=
1
(
2
)
=
gcd
(
a
,
b
)
(
3
)
=
gcd
(
b
,
a
(
m
o
d
b
)
)
(
4
)
=
b
x
′
+
(
a
m
o
d
b
)
y
′
(
5
)
=
b
x
′
+
(
a
−
b
⌊
a
b
⌋
)
y
′
(
6
)
=
a
y
′
+
b
(
x
′
−
⌊
a
b
⌋
)
\begin{aligned} (1)&&\quad ax + by &= 1 \\ (2)&&\quad &= \gcd(a,b) \\ (3)&&\quad &=\gcd(b, a(mod\ b)) \\ (4)&&\quad &=bx' + (a\ mod\ b)y' \\ (5)&&\quad &=bx' + (a - b\lfloor \frac{a}{b}\rfloor)y' \\ (6)&&\quad &=ay' + b(x' - \lfloor\frac{a}{b} \rfloor) \end{aligned}
(1)(2)(3)(4)(5)(6)ax+by=1=gcd(a,b)=gcd(b,a(mod b))=bx′+(a mod b)y′=bx′+(a−b⌊ba⌋)y′=ay′+b(x′−⌊ba⌋)
公式解释：
（1）题目给定等式
（2）a，b互素，则
1
=
gcd
(
a
,
b
)
1 = \gcd(a,b)
1=gcd(a,b)
（3）欧几里德定理，辗转相除法
（4）写出
b
b
b与
a
(
m
o
d
b
)
a(mod\ b)
a(mod b)的欧几里德定理形式
（5）根据取模的定义
a
m
o
d
b
=
a
−
b
⌊
a
b
⌋
a\mod b = a - b\lfloor\frac{a}{b} \rfloor
amodb=a−b⌊ba⌋（其中
⌊
a
b
⌋
\lfloor\frac{a}{b} \rfloor
⌊ba⌋表示向下取整）
（6）整理成关于a，b系数的形式
⇒
a
x
+
b
y
=
a
y
′
+
b
(
x
′
−
⌊
a
b
⌋
)
⇒
{
x
=
y
′
y
=
(
x
′
−
⌊
a
b
⌋
)
\Rightarrow ax+ by = ay' + b(x' - \lfloor \frac{a}{b} \rfloor) \\ \Rightarrow \begin{cases} x = y'\\ y = (x' - \lfloor \frac{a}{b} \rfloor) \end{cases}
⇒ax+by=ay′+b(x′−⌊ba⌋)⇒{x=y′y=(x′−⌊ba⌋)
得到了辗转相除法的递推公式，代码实现
编写成代码得：
typedef long long ll;
int EXgcd(ll a, ll b, ll &x, ll &y)
{
ll q, tmp
if(!b)
{
q = a;
x = 1;
y = 0;
}
else
{
q = EXgcd(b, a%b, x, y);
temp = x;
x = y;
y = temp - (a / b) * y;
}
return q;
}
求逆元题的标准形式
-
【例题1】给定正整数 a a a和 b b b，求满足等式 a x + b y = 1 ax + by = 1 ax+by=1的 x x x的最小正整数解。如果不存在返回 -1。
-
【例题2】给定正整数 a a a和 n n n，求满足等式 a x ≡ 1 ( m o d n ) ax \equiv 1(mod\ n) ax≡1(mod n)的 x x x的最小正整数解。如果不存在，则返回-1。
-
【例题3】给定素数 p p p和正整数 a a a，求满足 a ≡ 1 ( m o d p ) a\equiv 1(mod\ p) a≡1(mod p)的最小正整数 x x x，如果不存在返回-1。
对于例题3证明： 会用到二分快速幂
求逆元的方法
方法一：扩展欧几里得算法
公式原理：
a
∗
b
≡
1
(
m
o
d
p
)
a
∗
b
+
k
∗
p
=
1
a * b \equiv 1(mod\ p) \\ a *b +k * p = 1
a∗b≡1(mod p)a∗b+k∗p=1
其中 a a a为要求的逆元（其中b，p为已知），用辗转相除法求解。
辗转相除法（递归算法）递归推导过程见上
下以多项式证明（多项式空间与整数空间同构）
g
e
t
g
r
e
a
t
e
s
t
c
o
m
m
o
n
d
i
v
i
s
o
r
:
f
(
x
)
=
q
1
(
x
)
g
(
x
)
+
r
1
(
x
)
g
(
x
)
=
q
2
(
x
)
r
1
(
x
)
+
r
2
(
x
)
.
.
.
r
s
−
3
(
x
)
=
q
s
−
1
(
x
)
r
s
−
2
(
x
)
+
r
s
−
1
(
x
)
r
s
−
2
(
x
)
=
q
s
(
x
)
r
s
−
1
(
x
)
+
r
s
(
x
)
r
s
−
1
(
x
)
=
q
s
+
1
(
x
)
r
s
(
x
)
+
0
⇒
r
s
(
x
)
i
s
g
r
e
a
t
e
s
t
c
o
m
m
o
n
d
i
v
i
s
o
r
f
i
n
d
u
(
x
)
a
n
d
v
(
x
)
:
r
s
(
x
)
=
r
s
−
2
(
x
)
−
q
s
(
x
)
r
s
−
1
(
x
)
r
s
−
2
(
x
)
=
q
s
(
x
)
r
s
−
1
(
x
)
+
r
s
(
x
)
⇒
r
s
(
x
)
=
q
s
(
x
)
r
s
−
1
(
x
)
+
r
s
(
x
)
−
q
s
(
x
)
r
s
−
1
(
x
)
r
s
−
1
(
x
)
=
r
s
−
3
(
x
)
−
q
s
−
1
(
x
)
r
s
−
2
(
x
)
⇒
r
s
(
x
)
=
q
s
(
x
)
r
s
−
3
(
x
)
−
q
s
−
1
(
x
)
r
s
−
2
(
x
)
+
r
s
(
x
)
−
q
s
(
x
)
r
s
−
3
(
x
)
−
q
s
−
1
(
x
)
r
s
−
2
(
x
)
.
.
.
r
2
(
x
)
=
g
(
x
)
−
q
1
(
x
)
r
1
(
x
)
r
1
(
x
)
=
f
(
x
)
−
q
1
(
x
)
g
(
x
)
get \ greatest \ common \ divisor:\\ f(x) = q_1(x) g(x) + r_1(x) \\ g(x) = q_2(x) r_1(x) + r_2(x) \\ ... \\ r_{s-3}(x) = q_{s-1}(x) r_{s-2}(x) + r_{s-1}(x) \\ r_{s-2}(x) = q_{s}(x) r_{s-1}(x) + r_{s}(x) \\ r_{s-1}(x) = q_{s+1}(x) r_{s}(x) +0 \\ \Rightarrow r_{s}(x) \ is \ greatest \ common \ divisor \\ \ \\ find \ u(x) \ and \ v(x): \\ r_s(x) = r_{s-2}(x) - q_s(x)r_{s-1}(x) \\ r_{s-2}(x) = q_{s}(x) r_{s-1}(x) + r_{s}(x) \\ \Rightarrow r_s(x) = q_{s}(x) r_{s-1}(x) + r_{s}(x)- q_s(x)r_{s-1}(x) \\ r_{s-1}(x) = r_{s-3}(x) - q_{s-1}(x) r_{s-2}(x) \\ \Rightarrow r_s(x) = q_{s}(x) r_{s-3}(x) - q_{s-1}(x) r_{s-2}(x) + r_{s}(x)- q_s(x)r_{s-3}(x) - q_{s-1}(x) r_{s-2}(x)\\ ... \\ r_2(x) = g(x) - q_1(x)r_1(x) \\ r_1(x) = f(x) - q_1(x)g(x)
get greatest common divisor:f(x)=q1(x)g(x)+r1(x)g(x)=q2(x)r1(x)+r2(x)...rs−3(x)=qs−1(x)rs−2(x)+rs−1(x)rs−2(x)=qs(x)rs−1(x)+rs(x)rs−1(x)=qs+1(x)rs(x)+0⇒rs(x) is greatest common divisor find u(x) and v(x):rs(x)=rs−2(x)−qs(x)rs−1(x)rs−2(x)=qs(x)rs−1(x)+rs(x)⇒rs(x)=qs(x)rs−1(x)+rs(x)−qs(x)rs−1(x)rs−1(x)=rs−3(x)−qs−1(x)rs−2(x)⇒rs(x)=qs(x)rs−3(x)−qs−1(x)rs−2(x)+rs(x)−qs(x)rs−3(x)−qs−1(x)rs−2(x)...r2(x)=g(x)−q1(x)r1(x)r1(x)=f(x)−q1(x)g(x)
代码如下：
typedef long long ll;
ll ExpGcd(ll a, ll b, ll &x, ll &y)
{
ll q,tmp;
//边界条件
if(!b)
{
q = a;
x = 1;
y = 0;
}
// 被除数（作新a）
// 商（作新b）
// x - ay后的值（作新x）
// 余式（作新y）
else
{
q = ExpGcd(b, a % b, x, y);
tmp = x;
x = y;
y = tmp - (a / b) * y;
}
return q;
}
例如，对27a + 8k = 1方程（b = 27，p = 8）求解过程：
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-b3ZgpM9R-1680699079110)(./pic/divisor.png)]
费马小定理/欧拉定理
指数形式的逆元求法，一般式子中带有%
费马小定理原理：
if
p
is even, then
a
−
1
≡
1
(
m
o
d
p
)
⇒
a
p
−
2
≡
1
(
m
o
d
p
)
then
a
p
−
2
is the inverse element for a in mod p
\text{if }\ p \text{ is even, then } a^{-1} \equiv 1(mod\ p)\\ \Rightarrow a^{p - 2} \equiv 1(mod\ p) \\ \text{then } a^{p - 2} \text{ is the inverse element for a in mod p}
if p is even, then a−1≡1(mod p)⇒ap−2≡1(mod p)then ap−2 is the inverse element for a in mod p
欧拉定理原理：
if
a
,
p
are comprime then
a
φ
(
p
)
≡
1
(
m
o
d
p
)
Fermat’s Little Theorem
a
φ
(
p
)
−
1
∗
a
≡
1
(
m
o
d
p
)
a
φ
(
p
)
−
1
is the inverse element for a in mod p
\text{if }a,p \text{ are comprime then }a^{\varphi(p)} \equiv 1(mod\ p)\text{ Fermat's Little Theorem} \\ a^{\varphi(p) - 1} *a \equiv 1(mod\ p) \\ a^{\varphi(p) - 1}\text{ is the inverse element for a in mod p}
if a,p are comprime then aφ(p)≡1(mod p) Fermat’s Little Theoremaφ(p)−1∗a≡1(mod p)aφ(p)−1 is the inverse element for a in mod p
代码如下：
typedef long long ll;
//快速幂求解a^{p - 1}次幂
ll qkpow(ll a, ll p, ll mod)
{
ll t = 1, tt = a % mod;
while(p)
{
if(p & 1) t = t * tt % mod;
tt = tt* tt % mod;
p >>= 1;
}
return t;
}
llgetInv(ll a, ll mod)
{
return qkpow(a, mod - 2, mod);
}
递推求逆元
用于求很大范围内所有数的逆元（上述方法可能会超时）
原理：
set Mod
m
=
k
×
a
+
b
then
k
×
a
+
b
≡
0
(
m
o
d
m
)
multiplies
i
n
v
(
a
)
×
i
n
v
(
b
)
on both side
⇒
k
×
i
n
v
(
b
)
+
i
n
v
(
a
)
≡
0
(
m
o
d
m
)
∴
i
n
v
(
a
)
=
−
k
×
i
n
v
(
b
)
(
m
o
d
m
)
∵
k
=
m
/
a
,
b
=
m
%
a
⇔
i
n
v
(
a
)
=
−
(
m
/
a
)
×
(
m
%
a
)
(
m
o
d
m
)
\text{set Mod }m = k \times a + b \\ \text{then } k \times a + b \equiv 0(mod\ m) \\ \text{multiplies } inv(a) \times inv(b) \text{ on both side} \\ \Rightarrow k \times inv(b) +inv(a) \equiv 0 (mod\ m) \\ \therefore inv (a) = -k \times inv(b)(mod\ m) \\ \because k = m / a , b = m \% a \\ \Leftrightarrow inv(a) = -(m/a) \times (m \%a)(mod\ m)
set Mod m=k×a+bthen k×a+b≡0(mod m)multiplies inv(a)×inv(b) on both side⇒k×inv(b)+inv(a)≡0(mod m)∴inv(a)=−k×inv(b)(mod m)∵k=m/a,b=m%a⇔inv(a)=−(m/a)×(m%a)(mod m)
那么递推公式：
i
n
v
(
a
)
=
−
(
m
/
a
)
×
(
m
%
a
)
(
m
o
d
m
)
inv(a) = -(m/a) \times (m \%a)(mod\ m)
inv(a)=−(m/a)×(m%a)(mod m)
同余式（congruence）；
i
n
v
(
)
inv()
inv() 为求逆；
快速幂
二分递推的思想：
a
n
=
{
a
n
−
1
⋅
a
,
if n is odd
a
n
2
⋅
a
n
2
,
if n is even but not 0
1
,
if n = 0
a ^ n = \begin{cases} a^{n - 1} \cdot a, \ \text{if n is odd}\\ a^{\frac{n}{2}} \cdot a^{\frac{n}{2}}, \text{ if n is even but not 0} \\ 1, \text{if n = 0} \end{cases}
an=⎩
⎨
⎧an−1⋅a, if n is odda2n⋅a2n, if n is even but not 01,if n = 0
//递归快速幂
int qkpow(int a, int n)
{
if(n == 0)
{
return 1;
}
else if(n % 2 == 1)
{
return qkpow(a, n - 1) * a;
}
else
{
int tmp = qkpow(a, n / 2);
return tmp * tmp;
}
}
快速幂的优化
在实际问题中，题目常常会要求对一个大素数取模，这是因为计算结果可能会非常巨大，但是在这里考察高精度又没有必要。这时我们的快速幂也应当进行取模，此时应当注意，原则是步步取模，如果MOD较大，还应当开long long。
#define MOD 1000000007;
typedef long long ll;
ll qkpow (ll a, ll n)
{
if(n == 0)
return 1;
else if(n % 2 == 1)
return qkpow(a, n - 1) * a % MOD;
else
{
ll tmp = qkpow(a, n / 2) % MOD;
return tmp * tmp % MOD;
}
}
非递归快速幂
通过将指数写成二进制的形式，可以直接通过判断指数二进制的某一位（位运算）是否为1，来进行进行以上拆分的过程
Example：
算
7
1
0
7^10
710 将
10
10
10 转换成改写成十进制形式
(
1010
)
2
(1010)_2
(1010)2
那么接下来计算
7
(
1000
)
2
⋅
7
(
10
)
2
7^{(1000)_2} \cdot 7^{(10)_2}
7(1000)2⋅7(10)2
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-GP0g8wHe-1680699079110)(./pic/mi.webp)]
int qpow(int a, int n)
{
int ans = 1;
while(n)
{
if(n & 1)
ans *= a; //位上为1，进行乘法
a *= a; //确保每一位上乘的是该位的幂
n>>=1;
}
return ans;
}
快速幂的拓展
快速幂是一个模板，任何涉及到计算幂的计算都可以从其中抽象出来
//泛型的非递归快速幂
template <typename T>
T qpow(T a, ll n)
{
T ans = 1; // 赋值为乘法单位元，可能要根据构造函数修改
while (n)
{
if (n & 1)
ans = ans * a; // 这里就最好别用自乘了，不然重载完*还要重载*=，有点麻烦。
n >>= 1;
a = a * a;
}
return ans;
}
应用：矩阵快速幂求解斐波那契数列：
（洛谷P1962）斐波那契数列
已知斐波那契数列的递推关系式为
F n = { 1 ( n ≤ 2 ) F n − 1 + F n − 2 ( n ≥ 3 ) F_n = \begin{cases} 1&(n \leq 2) \\ F_{n - 1} + F_{n - 2} &(n \geq 3) \end{cases} Fn={1Fn−1+Fn−2(n≤2)(n≥3)
求 F n m o d 1 0 9 + 7 的值 F_n mod 10^9 + 7的值 Fnmod109+7的值
数学模型：
set matrix A = [ 0 1 1 1 ] , then A [ F n F n + 1 ] = [ F n + 1 F n + F n + 1 ] then [ F n F n + 1 ] = A [ F n − 1 F n ] = A 2 [ F n − 2 F n − 1 ] = . . . = A n − 1 [ F 1 F 2 ] = A n − 1 [ 1 1 ] \text{set matrix }A = \begin{bmatrix} 0 &1 \\ 1&1 \end{bmatrix}, \text{then } A\begin{bmatrix} F_n \\F_{n+1}\end{bmatrix} = \begin{bmatrix}F_{n+1}\\F_n + F_{n + 1}\end{bmatrix}\\ \text{then} \\ \begin{bmatrix} F_n \\ F_{n + 1} \end{bmatrix} = A\begin{bmatrix}F_{n - 1}\\ F_{n}\end{bmatrix} = A^2\begin{bmatrix}F_{n - 2}\\ F_{n - 1}\end{bmatrix} = ... = A^{n-1}\begin{bmatrix}F_{1}\\ F_{2}\end{bmatrix} = A^{n-1}\begin{bmatrix}1\\ 1\end{bmatrix} set matrix A=[0111],then A[FnFn+1]=[Fn+1Fn+Fn+1]then[FnFn+1]=A[Fn−1Fn]=A2[Fn−2Fn−1]=...=An−1[F1F2]=An−1[11]
于是斐波那契数列的大小就可以转化成求某个矩阵的幂的问题了。
代码如下：
#include <cstdio>
#define MOD 1000000007
typedef long long ll;
struct matrix
{
ll a1, a2, b1, b2;
matrix(ll a1, ll a2, ll b1, ll b2) : a1(a1), a2(a2), b1(b1), b2(b2) {}
matrix operator*(const matrix &y)
{
matrix ans((a1 * y.a1 + a2 * y.b1) % MOD,
(a1 * y.a2 + a2 * y.b2) % MOD,
(b1 * y.a1 + b2 * y.b1) % MOD,
(b1 * y.a2 + b2 * y.b2) % MOD);
return ans;
}
};
matrix qpow(matrix a, ll n)
{
matrix ans(1, 0, 0, 1); //单位矩阵
while (n)
{
if (n & 1)
ans = ans * a;
a = a * a;
n >>= 1;
}
return ans;
}
int main()
{
ll x;
matrix M(0, 1, 1, 1);
scanf("%lld", &x);
matrix ans = qpow(M, x - 1);
printf("%lld\n", (ans.a1 + ans.a2) % MOD);
return 0;
}
- 参考资料：
313
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
