转载自:https://blog.csdn.net/whereisherofrom/article/details/78922798
二、数论基础知识
1、欧几里德算法(辗转相除法)
2、扩展欧几里德定理
a.线性同余
b.同余方程求解
c.逆元
3、中国剩余定理(孙子定理)
4、欧拉函数
a.互素
b.筛选法求解欧拉函数
c.欧拉定理和费马小定理
5、容斥原理
二、数论基础知识
1、欧几里德定理(辗转相除法)
定理:gcd(a, b) = gcd(b, a % b)。
证明:a = kb + r = kb + a%b,则a % b = a - kb。令d为a和b的公约数,则d|a且d|b 根据整除的组合性原则,有d|(a-kb),即d|(a%b)。
这就说明如果d是a和b的公约数,那么d也一定是b和a%b的公约数,即两者的公约数是一样的,所以最大公约数也必定相等。
这个定理可以直接用递归实现,代码如下:
int gcd(int a, int b) {
return b ? gcd(b, a%b) : a;
}
return b ? gcd(b, a%b) : a;
}
这个函数揭示了一个约定俗成的概念,即任何非零整数和零的最大公约数为它本身。
【例题8】f[0] = 0, 当n>1时,f[n] = (f[n-1]+a) % b,给定a和b,问是否存在一个自然数k (0 <= k< b),是f[n]永远都取不到的。
永远有多远?并不是本题的范畴。
但是可以发现的是这里的f[...]一定是有循环节的,如果在某个循环节内都无法找到那个自然数k,那么必定是永远都找不到了。
求出f[n]的通项公式,为f[n] = an % b,令an = kb + r,那么这里的r = f[n],如果t = gcd(a, b),r = an-kb = t ( (a/t)n - (b/t)k ),则有t|r,要满足所有的r使得t|r,只有当t = 1的时候,于是这个问题的解也就出来了,只要求a和b的gcd,如果gcd(a, b) > 1,则存在一个k使得f[n]永远都取不到,直观的理解是当gcd(a, b) > 1,那么f[n]不可能是素数。
2、扩展欧几里德定理
a、线性同余
线性同余方程(也可以叫模线性方程)是最基本的同余方程,即ax≡b (mod n),其中a、b、n都为常量,x是未知数,这个方程可以进行一定的转化,得到:ax = kn + b,这里的k为任意整数,于是我们可以得到更加一般的形式即:ax + by + c = 0,这个方程就是二维空间中的直线方程,但是x和y的取值为整数,所以这个方程的解是一些排列成直线的点集。
b、同余方程求解
求解同余方程第一步是转化成一般式:ax + by = c,这个方程的求解步骤如下:
i) 首先求出a和b的最大公约数d = gcd(a, b),那么原方程可以转化成d(ax/d + by/d) = c,容易知道(ax/d + by/d)为整数,如若d不能整除b,方程必然无解,算法结束;否则进入ii)。
ii) 由i)可以得知,方程有解则一定可以表示成 ax + by = c = gcd(a, b)*c',那么我们先来看如何求解d = gcd(a, b) = ax + by,根据欧几里德定理,有:
d = gcd(a, b) = gcd(b, a%b) = bx' + (a%b)y' = bx' + [a-b*(a/b)]y' = ay' + b[x' - (a/b)y']
于是有x = y', y = x' - (a/b)y'。
由于gcd(a, b)是一个递归的计算,所以在求解(x, y)时,(x', y')其实已经利用递归计算出来了,递归出口为b == 0的时候(对比辗转相除,也是b == 0的时候递归结束),那么这时方程的解x0 = 1, y0 = 0。代码如下:
#define LL __int64
LL Extend_Euclid(LL a, LL b, LL &X, LL &Y) {
LL q, temp;
if( !b ) {
X = 1; Y = 0;
return a;
}else {
q = Extend_Euclid(b, a % b, X, Y);
temp = X;
X = Y;
Y = temp - (a / b) * Y;
return q;
}
}
LL Extend_Euclid(LL a, LL b, LL &X, LL &Y) {
LL q, temp;
if( !b ) {
X = 1; Y = 0;
return a;
}else {
q = Extend_Euclid(b, a % b, X, Y);
temp = X;
X = Y;
Y = temp - (a / b) * Y;
return q;
}
}
扩展欧几里德算法和欧几里德算法的返回值一致,都是gcd(a, b),传参多了两个未知数X, Y,采用引用的形式进行传递,对应上文提到的x, y,递归出口为b == 0,这时返回值为当前的a,因为gcd(a, 0) = a,(X, Y)初值为(1, 0),然后经过回溯不断计算新的(X, Y),这个计算是利用了之前的(X, Y)进行迭代计算的,直到回溯到最上层算法终止。最后得到的(X, Y)就是方程gcd(a, b) = ax + by的解。
通过扩展欧几里德求的是ax + by = gcd(a, b)的解,令解为(x0, y0),代入原方程,得:ax0 + by0 = gcd(a, b),如果要求ax + by = c = gcd(a, b)*c',可以将上式代入,得:ax + by = c = (ax0 + by0)c',则x = x0c', y = y0c',这里的(x, y)只是这个方程的其中一组解,x的通解为 { x0c' + kb/gcd(a, b) | k为任意整数 },y的通解可以通过x通解的代入得出。
【例题9】有两只青蛙,青蛙A和青蛙B,它们在一个首尾相接的数轴上。设青蛙A的出发点坐标是x,青蛙B的出发点坐标是y。青蛙A一次能跳m米,青蛙B一次能跳n米,两只青蛙跳一次所花费的时间相同。数轴总长L米。要求它们至少跳了几次以后才会碰面。
假设跳了t次后相遇,则可以列出方程:(x + mt) % L = (y + nt) % L
将未知数t移到等式左边,常数移到等式右边,得到模线性方程:(m-n)t%L = (y-x)%L (即 ax≡b (mod n) 的形式)
利用扩展欧几里德定理可以求得t的通解{ t0 + kd | k为任意整数 },由于这里需要求t的最小正整数,而t0不一定是最小的正整数,甚至有可能是负数,我们发现t的通解是关于d同余的,所以最后的解可以做如下处理:ans = (t0 % d + d) % d。
c、逆元
模逆元的最通俗含义可以效仿乘法,a*x = 1,则称x为a在乘法域上的逆(倒数);同样,如果ax≡1 (mod n),则称b为a模n的逆,简称逆元。求a模n的逆元,就是模线性方程ax≡b (mod n)中b等于1的特殊形式,可以用扩展欧几里德求解。并且在gcd(a, n) > 1时逆不存在。
3、中国剩余定理
上文提到了模线性方程的求解,再来介绍一种模线性方程组的求解,模线性方程组如图二-3-1所示,其中(ai, mi)都是已知量,求最小的x满足以下n个等式:
![](http://www.cppblog.com/images/cppblog_com/menjitianya/ysrj_sl_8.png)
图二-3-1
将模数保存在mod数组中,余数保存在rem数组中,则上面的问题可以表示成以下几个式子,我们的目的是要求出一个最小的正整数K满足所有等式:
K = mod[0] * x[0] + rem[0] (0)
K = mod[1] * x[1] + rem[1] (1)
K = mod[2] * x[2] + rem[2] (2)
K = mod[3] * x[3] + rem[3] (3)
... ...
这里给出我的算法,大体的思想就是每次合并两个方程,经过n-1次合并后剩下一个方程,方程的自变量取0时得到最小正整数解。算法描述如下:
i) 迭代器i = 0
ii) x[i] = (newMod[i]*k + newRem[i]) (k为任意整数)
iii) 合并(i)和(i+1),得 mod[i] * x[i] - mod[i+1] * x[i+1] = rem[i+1] - rem[i]
将x[i]代入上式,有 newMod[i]*mod[i]*k - mod[i+1] * x[i+1] = rem[i+1] - rem[i] - newRem[i]*mod[i]
iv) 那么产生了一个形如 a*k + b*x[i+1] = c的同余方程,
其中a = newMod[i]*mod[i], b = - mod[i+1], c = rem[i+1] - rem[i] - newRem[i]*mod[i]
求解同余方程,如果a和b的gcd不能整除c,则整个同余方程组无解,算法结束;
否则,利用扩展欧几里德求解x[i+1]的通解,通解可以表示成 x[i+1] = (newMod[i+1]*k + newRem[i+1]) (k为任意整数)
v) 迭代器i++,如果i == n算法结束,最后答案为 newRem[n-1] * mod[n-1] + rem[n-1];否则跳转到ii)继续迭代计算。
4、欧拉函数
a、互素
两个数a和b互素的定义为:gcd(a, b) = 1,那么如何求不大于n且与n互素的数的个数呢?
朴素算法,枚举i从1到n,当gcd(i, n)=1时计数器++,算法时间复杂度O(n)。
这里引入一个新的概念:用φ(n)表示不大于n且与n互素的数的个数,该函数以欧拉的名字命名,称为欧拉函数。
如果n是一个素数,即n = p,那么φ(n) = p-1(所有小于n的都互素);
如果n是素数的k次幂,即n = p^k,那么φ(n) = p^k - p^(k-1) (除了p的倍数其它都互素);
如果m和n互素,那么φ(mn) = φ(m)φ(n)(可以利用上面两个性质进行推导)。
将n分解成如图二-4-1的素因子形式,那么利用上面的定理可得φ(n)如图二-4-2所示:
![](http://www.cppblog.com/images/cppblog_com/menjitianya/ysrj_sl_9.png)
图二-4-1
![](http://www.cppblog.com/images/cppblog_com/menjitianya/ysrj_sl_10.png)
图二-4-2
前面已经讲到n的因子分解复杂度为O(k),所以欧拉函数的求解就是O(k)。
b、筛选法求解欧拉函数
由于欧拉函数的表示法和整数的素数拆分表示法很类似,都可以表示成一些素数的函数的乘积,所以同样可以利用筛选法进行求解。伪代码如下:
#define MAXP 2000010
#define LL __int64
void Eratosthenes_Phi() {
notprime[1] = true;
for(int i = 1; i < MAXP; i++) phi[i] = 1;
for(int i = 2; i < MAXP; i++) {
if( !notprime[i] ) {
phi[i] *= i - 1;
// 和传统素数筛法的区别在于这个i+i
for(int j = i+i; j < MAXP; j += i) {
notprime[j] = true;
int n = j / i;
phi[j] *= (i - 1);
while(n % i == 0) n /= i, phi[j] *= i;
}
}
}
}
#define LL __int64
void Eratosthenes_Phi() {
notprime[1] = true;
for(int i = 1; i < MAXP; i++) phi[i] = 1;
for(int i = 2; i < MAXP; i++) {
if( !notprime[i] ) {
phi[i] *= i - 1;
// 和传统素数筛法的区别在于这个i+i
for(int j = i+i; j < MAXP; j += i) {
notprime[j] = true;
int n = j / i;
phi[j] *= (i - 1);
while(n % i == 0) n /= i, phi[j] *= i;
}
}
}
}
这里的phi[i]保存了i这个数的欧拉函数,还是利用素数筛选将所有素数筛选出来,然后针对每个素因子计算它的倍数含有该素因子的个数,利用欧拉公式计算该素因子带来的欧拉函数分量,整个筛选过程可以参考素数筛选。
c、欧拉定理和费马小定理
欧拉定理:若n,a为正整数,且n,a互素,则:
。
![](http://www.cppblog.com/images/cppblog_com/menjitianya/ysrj_sl_11.png)
费马小定理:若p为素数,a为正整数且和p互素,则:
。
![](http://www.cppblog.com/images/cppblog_com/menjitianya/ysrj_sl_12.png)
由于当n为素数时φ(n) = p-1,可见费马小定理是欧拉定理的特殊形式。
证明随处可见,这里讲一下应用。
【例题10】整数a和n互素,求a的k次幂模n,其中k = X^Y, 正整数a,n,X,Y(X,Y<=10^9)为给定值。
问题要求的是a^(X^Y) % n,指数上还是存在指数,需要将指数化简,注意到a和n互素,所以可以利用欧拉定理,令X^Y = kφ(n) + r,那么kφ(n)部分并不需要考虑,问题转化成求r = X^Y % φ(n),可以采用快速幂取模,二分求解,得到r后再采用快速幂取模求解a^r % n。
5、容斥原理
容斥原理是应用在集合上的,来看图二-5-1,要求图中两个圆的并面积,我们的做法是先将两个圆的面积相加,然后发现相交的部分多加了一次,予以减去;对于图二-5-2的三个圆的并面积,则是先将三个圆的面积相加,然后减去两两相交的部分,而三个圆相交的部分被多减了一次,予以加回。
![](http://www.cppblog.com/images/cppblog_com/menjitianya/ysrj_sl_13.png)
图二-5-1
![](http://www.cppblog.com/images/cppblog_com/menjitianya/ysrj_sl_14.png)
图二-5-2
这里的“加”就是“容”,“减”就是“斥”,并且“容”和“斥”总是交替进行的(一个的加上,两个的减去,三个的加上,四个的减去),而且可以推广到n个元素的情况。
【例题11】求小于等于m(m < 2^31)并且与n(n < 2^31)互素的数的个数。
当m等于n,就是一个简单的欧拉函数求解。
但是一般情况m都是不等于n的,所以可以直接摈弃欧拉函数的思路了。
考虑将n分解成素数幂的乘积,来看一种最简单的情况,当n为素数的幂即n = p^k时,显然答案等于m - m/p(m/p表示的是p的倍数,去掉p的倍数,则都是和n互素的数了);然后再来讨论n是两个素数的幂的乘积的情况,即n = p1^k1 * p2^k2,那么我们需要做的就是找到p1的倍数和p2的倍数,并且要减去p1和p2的公公倍数,这个思想其实已经是容斥了,所以这种情况下答案为:m - ( m/p1 + m/p2 - m/(p1*p2) )。
类比两个素因子,如果n分解成s个素因子,也同样可以用容斥原理求解。
容斥原理其实是枚举子集的过程,常见的枚举方法为dfs,也可以采用二进制法(0表示取,1表示不取)。这里给出一版dfs版本的容斥原理的伪代码,用于求解小于等于m且与n互素的数的个数。
#define LL __int64
void IncludeExclude(int depth, LL m, LL mul, int op, int* p, LL &ans) {
if(m < mul) return ;
if(depth == p[0]) {
ans += (op ? -1 : 1) * (m / mul);
return ;
}
for(int i = 0; i < 2; i++) {
// 0 表示不取, 1表示取
IncludeExclude( depth+1, m, mul * (i?p[depth+1]:1), op^i, p, ans );
}
}
void IncludeExclude(int depth, LL m, LL mul, int op, int* p, LL &ans) {
if(m < mul) return ;
if(depth == p[0]) {
ans += (op ? -1 : 1) * (m / mul);
return ;
}
for(int i = 0; i < 2; i++) {
// 0 表示不取, 1表示取
IncludeExclude( depth+1, m, mul * (i?p[depth+1]:1), op^i, p, ans );
}
}
p[ 1 : p[0] ]存储的是n的所有素因子,p[0]表示数组长度,mul表示该次的素因子子集的乘积,op表示子集的奇偶性,ans存储最后的答案。
例如求[1, 9]中和6互素的数的个数,这时p = [2, 2, 3] (注意p[0]是存素数的个数的,6分解的素因子为2和3)。
ans = 9/1 - (9/2 + 9/3) + 9/6 = 3,ans分为三部分,0个数的组合,1个数的组合,2个数的组合。