辗转相除法和有理数取模

七海

算法数学

发布于：2020年2月12日

次浏览

首先简单复习一下辗转相除法（欧几里得算法）和扩展欧几里得算法：

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

更多练习题

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

参考资料

• 《算法竞赛入门经典（第2版）》 - 刘汝佳编著，清华大学出版社
  
• https://www.cnblogs.com/zhanhonhao/p/11329772.html
• https://www.luogu.com.cn/problemnew/solution/P2613

更新于：2020年2月12日

算法

ACM

数论