抱歉,您的浏览器无法访问本站

本页面需要浏览器支持(启用)JavaScript


了解详情 >

七つの海

今日、海を見た。もう怖くない

首先简单复习一下辗转相除法(欧几里得算法)和扩展欧几里得算法:

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
// 最大公约数 & 最小公倍数
int gcd(int a, int b)
{
if(a<b) return gcd(b,a); // 确保 a>b
else if(!b) return a;
else return gcd(b,a%b);
}

int lcm(int a, int b)
{
return a/gcd(a,b)*b; // 如果是 a*b/gcd(a,b),中间结果可能会爆
}

// 拓展欧几里得: 解不定方程 ax+by=gcd(a,b)
void exgcd(int a, int b, int& d, int&x, int& y)
{
if(a<b) return exgcd(b,a,d,y,x);
else if(!b) // 基础情况: 1*a-0*0 = a = gcd(a,0)
{
d = a; // gcd(a,b)
x = 1; // a 的系数
y = 0; // b 的系数
}
else
{
exgcd(b,a%b,d,y,x);
y -= x*(a/b);
}
}

根据写法流派的不同,exgcd的参数d(最大公约数)可以作为返回值返回。

扩展欧几里得

说这个之前先说一下贝祖定理

贝祖定理:如果a、b是整数,那么一定存在整数x、y使得ax+by=gcd(a,b)。

也就是说,如果ax+by=m有解,那么m一定是gcd(a,b)的若干倍。特别地,如果ax+by=1有解,那么gcd(a,b) = 1。直接的应用就是知道不定方程的情况下,判断这样的等式有没有整数解。

但是这样并不能得到这个不定方程的解,这种时候就需要使用扩展欧几里得算法。

扩展欧几里得算法可以在求得a、b的最大公约数的同时,能找到整数x、y(其中一个很可能是负数),使它们满足贝祖等式 ax+by=gcd(a,b)。如果a是负数,可以把问题转化成 |a|(-x')+by=gcd(a,b),然后令 x=(-x')。

和欧几里得算法一样,扩展欧几里得也是使用递归求解,有一样的递归边界:确保a>b的前提下,当b = 0的时候,GCD的值是a,此时对于不定方程 a*x + b*y = gcd(a,b) 有显而易见的解(x,y) = (1,0)。

现在我们可以从这个基本情况进行反推,求出一般a和b的可能解(x,y)。假设其中的某一阶段是这样的:已知 b、a%b 和 gcd(b,a%b)=gcd(a,b)=d,已经求解得到了 b*x' + (a%b)*y' = d。接下来要继续反推求出 ax + by = d 的x和y:

显然,a%b = a - (a/b)*b。将这个式子待入到上述状态的不等方程中可以得到:

d = b * x' + (a%b) * y' d = b * x' + (a - (a/b) * b) * y' = b * x' + a * y' - (a/b) * b * y'
d = a * y' + b * (x' - (a/b) * y')

小声bb:我马上就去学习怎么输入数学公式

综上所述,可以得到:若 b*x' + (a%b)*y' = d,那么解为 x = y';y = x' - (a/b)*b,即 a*y' + b*(x'-(a/b)*b) = d。递归返回的过程中根据这个更新值就可以。

更新方法:x和y是引用值,每层递归要使用下一层计算出的y作为自己的x,使用下一层计算的x作为自己的y后再进行本层值的处理。所以有了y -= x*(a/b);

扩展欧几里得算法可以用来计算模反元素(也叫模逆元),而模反元素在RSA加密算法中有重要的作用。

模逆元

再说这个之前先简单的复习一下同余和模算数的相关知识:

同余: a ≡ b (mod n),含义是 a, b除以 n的余数相同。它的充要条件是 a-b 可以被 n整除。若有整数 y满足 a ≡ y (mod n),那么 y ≡ b (mod n)。所以实际上同余方程的一个解指的是一整个同余等价类。

模线性方程组: 形如 ax ≡ b (mod n) 的方程。由上述等价条件可知 ax - b = ny,再进行移项可得 ax - ny = b,正是可以被扩展欧几里得求解的不等方程的形式。先对 a, n 进行 exgcd,若 d不能被 b整除,那么无解;否则将求出的 x, y按照比例放大即可。

模算数: 同余也有四则运算
- 加法: (a + b) % n = ((a%n) + (b%n)) % n - 减法: (a - b) % n = ((a%n) - (b%n) + n) % n - 乘法: ab % n = (a%n)(b%n) % n - 除法:<无> - 同余式: c≠0, ac ≡ bc (mod n) => a ≡ b (mod n/gcd(c,n))

特别注意,最后一项不能简单理解为除法。最后一项的特点是都可以整除。

费马小定理: 若p是质数,且a、p互质,那么a^(p-1) mod p = 1。若p是质数,有 a^p ≡ a (mod p)。两式在 p|a(p整除a)时不等价。

那么什么是模逆元?就是上述模线性方程组的 b=1 的特殊情况:当 ax ≡ 1 (mod n) 时,解出的 x 就是 a 关于模 n 的逆。在这种情况下,a 的逆存在的必要条件就是该模线性方程组(或其等价形式 ax - ny = 1)有解。而它有解的必要条件是 1|gcd(a,n),也就是 a 和 n 互素,gcd(a,n)=1。此时,该方程有唯一解。(唯一解是刘汝佳书上的原话,但是不是这样的)

此时,可以利用模逆元来求模算数的除法。虽然模除法相当于对一个分数有理数取模,意义不明,但是它依然满足同余的性质。若已有 x 使得 bx ≡ 1 (mod n),又由模乘法可以得到: (a/b) % m = (a/b)*1 % m = (a/b)*bx % m = ax % m,即 (a/b) ≡ ax (mod m),式中 x 是 b 关于模 m 的逆。

使用扩展欧几里得求模逆元

求解这个逆元就要使用到前面说的扩展欧几里得算法。bx ≡ 1 (mod n) 等价于 bx - ny = 1。若 gcd(b,n)=1,那么可以使用扩展欧几里得求出逆元 x。因为实际上线性不定方程组有无穷多解,这里只求正的最小的逆元。代码如下:

1
2
3
4
5
6
7
8
9
10
// 使用扩展欧几里得求 a 关于模 m 的逆元
int inverse(int a, int m)
{
int x,y,d; // 等价: ax - ny = 1
exgcd(a,m,d,x,y);
if(1%d) return -1; // gcd(a,m)≠1 时,无解
x *= 1/gcd; // 解不定方程通用:放大倍数
if(x<0) x+=m; // 求出负数时,将它加为正数
return x%m;
}

上述代码并没有处理 m 是负数的情况。如果 m 是负数,那么在 exgcd 之后将 m 取绝对值就好了。

使用费马小定理求模逆元

当模数 m 是素数的时候,求模逆元也可以使用费马小定理进行。很多题目在使用模数时都会给一个很大的素数,与任何比它小的正数互质,就可以使用。设这个很大的质数是 p,那么根据费马小定理,有 a^(p-1) ≡ 1 (mod p),a 是任何比 p 小的正数。那么根据模逆元的定义,显然有 a*a^(p-2) ≡ 1 (mod p),故 a 关于模 p 的逆元是 a^(p-2)。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// 快速幂: 计算 a^n (n>=0),非取模
int fastpow(int a, int n)
{
if(!n) return a?1:0;
int ans = 1;
while(n)
{
if(n&1) ans*=a;
a*=a;
n>>=1;
}
return ans;
}

// 使用费马小定理求 a 关于模 p 的逆元
int inverse(int a, int p)
{
if(a>p||a<=0) return -1;
return fastpow(a,p-2)%p; // 返回模值
}

特别注意:上述代码中的快速幂,如果要取模,即使结果不超过int范围,考虑到中间结果的可能性,参数和中间变量的类型应当改为long long

这样的话,就可以很方便的解决一些带模的除法问题了,而不是使用高精度计算了。比如求 x = (a/c) % p,可以简单的转化为 a*c^(-1) % p。又因为费马小定理,1 ≡ c^(p-1) (mod p),原式可以化为 x = a*c^(p-2) % p。这样就可以完成分数的求模了。

分数取模

首先看定义:

分数取模运算的定义: x = a/b (mod m) <=> xb = a (mod m),仅当 gcd(b,m)=1 时成立

当上述恒等式成立时,求已知分数的模相当于求解等价的不等方程: xb - my = a。这可以使用扩展欧几里得求解。当 m 是质数的时候,也可以使用费马小定理快速求解。

但是现在有一些题目会给分数取模的整数作为输入数据。这种题目如果是恰好需要计算概率的话,就很容易陷入“将取模后分数转化为分数”的误区中。但是实际上将取模后分数转化为分数并不是分数取模的逆运算——分数的模只是不定方程的一个解,另一个解被舍弃了。然而题目并不会告诉逆被舍弃的解。但是因为分数取模和分数满足同余等式,可以直接使用分数的模当作分数进行计算,只是要注意负数的情况即可。

练习题

一些相关的题目。写了我就放上来()

洛谷 P2613

题目链接:洛谷 P2613【模板】有理数取余

真正的模板题,直接告诉了要对分数 a/b 取模。但是数据范围着实吓人,不像是可以用一个变量装下的。如果要是打高精度的话风险太大,也不是这个题目的正确解法。因为分数取模最后可以通过逆元转化为模乘法,所以运算依然是线性的。可以修改C快读模板,在C快读模板每读入一位数字时进行取模操作,最后就可以得到待输入数字的关于给定数字的模。然而这里取模对题目是不影响的。

检查一下,发现题目给的19260817是一个素数。这太好了,直接费马小定理就完事了。当然也可以使用扩展欧几里得解不定方程解出。因为实际上是分母和大质数求不定方程,有解仅当它们互质,那么易得无解的情况:分母是0,根本就没有和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
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
#include <iostream>
#include <cstdio>
#include <cctype>

using namespace std;

typedef long long longs;

const longs MOD = 19260817LL; // 是素数
const int Power = (int)MOD-2;
const char No[] = "Angry!";

longs a,b;

longs getlongs() // 未处理负数的带模快读
{
longs ll = 0ll;
char ch;
ch = getchar();
while(!isdigit(ch) && ch != EOF)
ch = getchar();
while(isdigit(ch))
{
ll = (ll<<1)+(ll<<3)+ch-'0';
ll %= MOD; // 立即取模
ch = getchar();
}
return ll;
}

int mod(longs a, longs b) // 费马小定理法
{
int n = Power;
longs ans = 1;
while(n)
{
if(n&1) ans = (ans*b)%MOD;
b = (b*b)%MOD;
n >>= 1;
}
ans = (ans*a)%MOD;
return (int)ans;
}

int main()
{
// 要使用C快读,必须关闭C++快读
// ios::sync_with_stdio(false);
// cin.tie(nullptr);

a = getlongs();
b = getlongs();

if(!b) cout<<No<<endl;
else cout<<mod(a,b)<<endl;

return 0;
}

HDU - 4828

题目链接:HDU - 4828 Grids

更多练习题

如果遇到好题目以后会补在这里。如果你们有推荐的题目也欢迎在评论区留言~

参考资料

评论