数论及数学相关

数论

Gcd 以及 EX_Gcd

学习了一下gcd和拓展gcd。之前一直以为是十分难得东西, 其实好像没有那么复杂。
Gcd(最小公倍数) 我们使用欧几里得算法来求解最小公倍数问题, 代码如下:

1
2
3
int gcd(int x, int y ){
return (b==0) ? a : gcd(b,a%b);
}

这个一行的代码可以快速求出两数的gcd, 时间复杂度为O(logn),效率还算不错。

lcm 如果我们已经求出了两个数的gcd,那么只需要一个式子就可以求出两数的lcm了:

1
2
3
inline int lcm(int a, int b) {
return a/gcd(a, b)*b;
}//lcm(a, b) = a*b/gcd(a, b);

引理

有一个叫做裴蜀等式或者贝祖等式的东西, 可以方便我们来求拓展gcd。 大概内容就是 ax+by=z 有整数解当且仅当 gcd(a,b) | z时存在整数解。并且当裴蜀等式有解时该方程必然有无穷多个整数解,每组解x、y都称为裴蜀数,可用扩展欧几里得算法(Extended Euclidean algorithm)求得。

然后就是拓展欧几里得算法了, 这个算法的目的是求出类似于ax+by=z 形势的一个二元一次方程的可行解,因为二元一次方程是不定方程, 可能存在多组解的情况。而拓展gcd就是求出一组可行解。证明如下, 之后再加上代码:

从网上摘了一个证明:

1
2
3
4
5
6
7
8
9
10
11
设:a>b。
  推理1,显然当 b=0,gcd(a,b)=a。此时 x=1,y=0;//推理1
  推理2,ab!=0
  设 ax1+by1=gcd(a,b);
  bx2+(a mod b)y2=gcd(b,a mod b);
  根据朴素的欧几里德原理有 gcd(a,b)=gcd(b,a mod b);
  则:ax1+by1=bx2+(a mod b)y2;
  即:ax1+by1=bx2+(a-(a/b)*b)y2=ay2+bx2-(a/b)*by2;
  根据恒等定理得:x1=y2 ,y1=x2-(a/b)*y2;//推理2
  这样我们就得到了求解 x1,y1 的方法:x1,y1 的值基于 x2,y2.
上面的思想是以递归定义的,因为 gcd 不断的递归求解一定会有个时候 b=0,所以递归可以结束。

这样我们就找到了递推关系:

递推的终止条件再上面也已经给出:

Code:

1
2
3
4
5
6
7
8
9
10
11
12
void Ex_gcd(int a, int b, int& d, int &x, int &y) {
if (b == 0) {
d = a;
x = 1;
y = 0;
return;
}
int x1, y1;
Ex_gcd(b, a%b, d, x1, y1);
x = y1;
y = x1 - (a/b)*y1;
}// 此时x 与 y中储存的即为这个方程的一组解, d为顺带求得一下 a 与 b 的最小公倍数

附带一个很良心的手推过程:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
exgcd(47, 30, x, y)
{
r = exgcd(30, 17,x, y)
{
r = exgcd(17, 13, x, y)
{
r = exgcd(13, 4, x, y)
{
r = exgcd(4, 1, x, y)
{
r = exgcd(1, 0, x, y)
{
x = 1;
y = 0;
return 1;
}
t = y = 0;
y = x - (4/1) * y = 1;
x = t = 0;
return r = 1;
}
t = 1;
y = 0 - (13/4) * 1 = -3;
x = 1;
return 1;
}
t = -3;
y = 1 - (17/13) * (-3) = 4;
x = -3;
return 1;
}
t = 4;
y = -3 - (30/17) * 4 = -7;
x = 4;
return 1;
}
t = -7;
y = 4 - (47/30) * (-7) = 11;
x = -7;
return 1;
}
最后的结果:
r = exgcd(47,30,x,y) = 1;
x = -7;
y = 11;

乘法逆元

下面来说一个和拓展欧几里得定理十分有关的一个东西就是乘法逆元,我们一般在模算术中出现过膜的加法, 膜的减法, 膜的乘法等。我们把一些数被一个数膜的集合叫做这个数的完全剩余系,写为:

我们还定义如果n的完全剩余系中的数字还和n互素, 那么这些数组成的集合叫做n的是简化剩余系又叫做简系

但是我们在做题的会出现一些除法取模的问题, 这时候我们发现在数论的领域里没有分数这种东西, 这时候我们会发现我们需要寻找ab = 1的b这个数,也可以叫做a^-1。 而这个b也就叫做a在膜n下的逆, 也叫作乘法逆元。

1. 拓展欧几里得

乘法逆元可以有很多种做法, 先写一下比较一般的写法:

这个式子的推导就把a的逆元x用不定方程表示了出来, 直接把这个方程带进exgcd去解即可:

1
2
3
4
5
int inv(int a, int p) {
int x, y, d;
exgcd(a, p, d, x, y);
return d == 1 ? (a+p)%p : -1;
}

如果返回负一则这个数没有逆元 复杂度O(log n) 。

2. 费马小定理

费马小定理的式子是这样的:

所以当模数是一个质数的时候,可以用费马小定理求解,即

直接快速幂: power_mod(a, p-2, mod);(一定要记清这个不是一般情况的求法, 仅限于膜数是一个质数的时候)

1
2
3
4
5
6
7
8
9
10
11
12
int power(int a, int b, int p) {
int res = 1 % p;
while(b) {
if (b&1) res = (ll) res * a % p ;
a = (ll) a * a % p;
b >>= 1;
}
return res;
}
inline int niyuan(int a, int b) {
return power(a, b-2, b);
}

3.线性递推

O(n)O(n)的时间可以处理出1~n 在mod p 意义下的逆元。那么我们就可以做到线性递推 :

1
inv[i] = (mod-mod/i) * inv[mod % i] % mod;

线性递推的证明过程:

1
2
3
4
5
6
7
8
设 x = p % a,y = p / a;  (y就是p/a的商, x就是p/a后的余数) 
于是: x + y * a = p;
等式两边同时对p取余: (x + y * a)% p = 0;
去括号并移项:x % p = - y * a % p;
将a移至左边: inv (a) * x % p = - y % p;
       inv(a) * x % p = (-y + p ) % p;
将x移至右边 inv(a) =( p - y ) * inv(x) % p;
将证明开始设的x,y值代入上面式子中,则有 inv(a) = ( p - p/a )*inv( p % a ) % p;

网上扒了一段代码粘上来: 几种求逆元的方法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
#define LL long long
int inv[1000010];
LL ksm(LL a,LL b,LL mod) {
int ans=1;
while(b) {
if(b&1) ans=(ans*a)%mod;
a=(a*a)%mod;
b>>=1;
}
return ans;
}
LL exgcd(LL a,LL b,LL &x,LL &y) {
if(!b) {
x=1;
y=0;
return a;
}
LL GCD=exgcd(b,a%b,x,y);
LL tmp=x;
x=y;
y=tmp-a/b*y;
return GCD;
}
LL inv1(LL a,LL mod) {//扩展欧几里得求逆元
LL x,y;
LL d=exgcd(a,mod,x,y);
if(d==1) return (x%mod+mod)%mod;
return -1;
}
LL inv2(LL a,LL mod) {//费马小定理
return ksm(a,mod-2,mod);
}
void inv3(LL mod) {//线性递推求逆元
inv[1]=1;
for(int i=2;i<=mod-1;i++)
{
inv[i]=(mod-mod/i)*inv[mod%i]%mod;
cout<<inv[i]<<" ";
}
}

中国剩余定理

​ 本篇博客正在编写中…