十三、数论基础

模运算

模运算是大数运算中的常用操作。如果一个数太大，无法直接输出，或者不需要直接输出，可以把它取模后缩小数值再输出。定义取模运算为a除以m的余数,记为：
$$amodm=a\%m$$
取模的结果满足$0\le a\%m\le m-1$，题目用给定的$m$限制计算结果的范围。例如$m=10$,就是取计算结果的个位数。取模操作满足以下性质:
$$加：(a+b)\%m=((a\%m)+(b\%m))\%m$$
$$减：(a-b)\%m=((a\%m)-(b\%m))\%m$$
$$乘：(a\times b)\%m=((a\%m)\times (b\%m))\%m$$
然而，对除法取模进行下面的类似操作是错误的：
$$(a/b)\%m=((a\%m)/(b\%m))\%m$$
如：$(100/50)\%20=2，(100\%20)/(50\%20)\%20=0$,两者不相等。除法的取模需要用到逆元。

快速幂

快速幂以及扩展的矩阵快速幂，由于应用场景比较常见，也是竞赛中常见的题型。幂运算$a^n$即$n$个$a$相乘。快速幂就是高效地算出$a^n$。当$n$很大时，例如：$n=10^9$,计算$a^n$这样大的数，一是数字太大，二是计算时间很长。
下面先考虑如何缩短计算时间，如果用暴力的方法直接算$a^n$,即逐个做乘法，复杂度是$O(n)$,即使能算出来，也会超时。
读者很容易想到快速幂的办法：先算$a^2$,然后继续算平方$(a^2{)}^2$,一直算到$n$次幂。这
是分治法的思想，复杂度为$O(lo{g}_{2}n)$。下面是代码，请读者自己理解：

int fastPow(int a, int n) {
    if (n == 1) return a;
    int temp = fastPow(a, n / 2);   //分治
    if (n % 2 == 1) {               //奇数个a,此处也可以写为if(n &1)
        return temp * temp * a;
    } 
    else {                          //偶数个a
        return temp * temp;
    }
}

程序中的递归，层数只有$lo{g}_{2}n$，不用担心溢出的问题。

矩阵快速幂

给定一个$m\times m$的矩阵$A$,求它的$n$次幂$A^n$,这也是常见的计算。同样有矩阵快速幂的算法，原理是把矩阵当作变量来操作，程序和上面的很相似。
首先需要定义矩阵的结构体，并且定义矩阵相乘的操作。注意矩阵相乘也需要取模：

const int MAXN = 2;                     //根据题目要求定义矩阵的阶，本例中是2
const int MOD = 1e4;                    //根据题目要求定义模
struct Matrix {                         //定义矩阵
    int m[MAXN][MAXN];
    Matrix() {                          //初始化矩阵
        memset(m, 0, sizeof(m));
    }
};

Matrix Multi(Matrix a, Matrix b) {      //矩阵的乘法
    Matrix res;
    for (int i = 0; i < MAXN; i++)
        for (int j = 0; j < MAXN; j++)
            for (int k = 0; k < MAXN; k++)
                res.m[i][j] = (res.m[i][j] + a.m[i][k] * b.m[k][j]) % MOD;
    return res;
}

下面是矩阵快速幂的程序代码，和前面单变量的快速幂的代码非常相似：

Matrix fastm(Matrix a, int n) {
    Matrix res;
    for (int i = 0; i < MAXN; i++) res.m[i][i] = 1;
			//初始化为单位矩阵，相当于前面程序中的"int res =1;"
    while (n) {
        if (n & 1)res = Multi(res, a);
        a = Multi(a, a);
        n /= 2;
    }
    return res;
}

最大公约数，最小公倍数

最大公约数GCD和最小公倍数LCM是竞赛中常见的知识点，虽然这两个知识点很容易理解，但往往会与其他知识点结合起来出综合题。

int GCD(int a, int b){
    if(a % b == 0)return b;
    return GCD(b, a % b);
}

int LCM(int a, int b){
    return a * b / GCD(a, b);
}

扩展欧几里得算法与二元一次方程的整数解

给出整数$a、b、n$,问方程$ax+by=n$什么时候有整数解？如何求所有的整数解？
有解的充分必要条件是：$gcd(a,b)$可以整除$n$。简单证明如下：

令：a = gcd(a, b)a'、b = gcd(a, b)b';
有：ax + by = gcd(a, b) (a'x + b'y) = n;
如果x、y、a'、b'都是整数，
那么：n必须是gcd(a,b)的倍数才有整数解。

例如$4x+6y=8$、$2x+3y=4$有整数解，$4x+6y=7$则没有整数解。如果确定有解，一种解题方法是先找到一个解$(x_0,y_0)$,那么通解公式如下：
$$\{\begin{array}{l}x=x_0+bt\\ y=y_0-att为任意整数\end{array}$$
所以，问题转化为如何求$(x_0,y_0)$。利用扩展欧几里得算法可以求出这个特解。

扩展欧几里得算法

当方程符合$ax+by=gcd(a,b)$时，可以用扩展欧几里得算法求$(x_0,y_0)$。程序如下：

void ex_gcd(int a, int b, int &x, int &y) { //x, y变化后形成特解（x0, y0）
    if (b == 0) {
        x = 1, y = 0;
        return;
    }
    ex_gcd(b, a % b, x, y);
    int tmp = x;
    x = y;
    y = tmp - (a / b) * y;
}

有时候为了简化描述，在$ax+by=gcd(a,b)$两边除以$gcd(a,b)$，得到$cx+dy=1$。
其中：$c=a/gcd(c,b),d=b/gcd(a,b)$。
很明显，$c、d$是互质的。$cx+dy=1$的通解如下：
$$\{\begin{array}{l}x=x_0+dt\\ y=y_0-ctt为任意整数\end{array}$$

求任意方程ax+by=n的一个整数解

用扩展欧几里得算法求解$ax+by=gcd(a,b)$后，利用它可以进一步解任意方程$ax+by=n$，得到一个整数解。其步骤如下：

1. 判断方程$ax+by=n$是否有整数解，有解的条件是$gcd(a,b)$可以整除$n$
2. 用扩展欧几里得算法求$ax+by=gcd(a,b)$的一个解$(x_0,y_0)$
3. 在$a{x}_0+b{y}_0=gcd(a,b)$两边同时乘以$n/gcd(a,b)$得:$a{x}_{0}n/gcd(a,b)+b{y}_{0}n/gcd(a,b)=n$
4. 对照$ax+by=n$,得到它的一个解$(x‘,y‘)$是：
   $$\{\begin{array}{l}x‘=x_0\ast n/gcd(a,b)\\ y‘=y_0\ast n/gcd(a,b)\end{array}$$

应用场合

扩展欧几里得算法是一个很有用的工具，在竞赛题目中常用于以下场合：

1. 求解不定方程
2. 求解模的逆元
3. 求解同余方程
   虽然用扩展欧几里得算法可以算$ax+by=gcd(a,b)$的通解，不过一般没有这个需求，而是用于求某些特殊的解，例如求解逆元，逆元是除法取模操作常用的工具。

同余与逆元

同余在数论中非常有用，它用类似处理等式的方式来处理整除关系，非常简便。

同余

两个整数$a、b$和一个正整数$m$,如果$a$除以$m$所得的余数和$b$除以$m$所得的余数相等。
即：$a\%m=b\%m$，称$a$和$b$对于$m$同余，$m$称为同余的模。
同余的概念也可以这样理解：$m|(a-b)$。
即：$a-b$是$m$的整数倍。例如$6|(23-5)$，$23$和$5$对模$6$同余。
同余的符号记为$a\equiv b(modm)$

一元线性同余方程

$ax\equiv b(modm)$,即$ax$除以$m,b$除以$m$的余数相同，这里$a、b、m$都是整数，求解$x$的值。
方程也可以这样理解：$ax-b$是$m$的整数倍。设$y$是倍数，那么$ax-b=my$。
移项得到：$ax-my=b$。因为$y$可以是负数，改写为：$ax+my=b$。
这就是在扩展欧几里得算法中的二元一次不定方程。
当且仅当$gcd(a,m)$能整除$b$时有整数解。例如$15x+6y=9$，有整数解$x=1,y=-1$。当$gcd(a,m)=b$时，可以直接用扩展欧几里得算法求解$ax+my=b$。
如果不满足$gcd(a,m)=b$，还能用扩展欧几里得算法求解$ax+my=b$吗？答案是肯定的，但是需要结合逆元。

逆元

给出$a$和$m$，求解方程$ax\equiv 1(modm)$，即$ax$除以$m$余数是$1$。
根据前面的讨论，有解的条件是$gcd(a,m)=1$，即$a、m$互质。该问题等价于求解$ax+my=1,可以用上一节的扩展欧几里得算法求解。例如$8x\equiv 1(mod31)$,等价于求解$8x+31y=1$，用扩展欧几里得算法求得一个特解是$x=4,y=-1$。
方程$ax\equiv 1(modm)$的一个解$x$，称$x$为$a$模$m$的逆。注意，这样的$x$有很多，把它们都称为逆。
求逆元的代码如下：

	int inverse_mod(int a, int m){  
		int x, y;  
		ex_gcd(a, m, x, y);  
		return(m + x % m ) % m;  
}

试除法判断质数

问题：输入一个很大的数$n$，判断它是不是素数。
素数定义：一个数$n$，如果不能被$[2,n-1]$内的所有数整除，$n$就是素数。当然，并不需
要把$[2,n-1]$内的数都试一遍，这个范围可以缩小到$[2,\sqrt{n}]$。
给定$n$,如果它不能整除$[2,\sqrt{n}]$内的所有数，它就是素数。证明如下：

设$n=a\times b$，有$min(a,b)\le \sqrt{n}$，令$a\le b$。只要检查$[2,\sqrt{n}]$内的数，如果$n$不是素数，就能找到一个$a$。如果不存在这个$a$，那么$[\sqrt{n},n-1]$内也不存在$b$。

以上判断的范围可以再缩小一点：$[2,\sqrt{n}]$内所有的素数。其原理很简单，读者在学过下文提到的埃式筛法之后更容易理解。
用试除法判断素数，复杂度是$O(\sqrt{n})$，对于$n\le 10^{12}$的数是没有问题的。下面是试除法判断素数的代码：

bool is_pri(long long n){  
	if(n <= 1)return false;  
	for(int i = 2; i * i <= n; i++){  
		if(n % i == 0)return false;  
	}  
	return ture;  
}