数论基础算法

最新推荐文章于 2022-12-28 16:54:40 发布

斯文刘

最新推荐文章于 2022-12-28 16:54:40 发布

阅读量2.3k

点赞数 1

文章标签：算法加密 c algorithm n2 behavior

本文链接：https://blog.csdn.net/lzj888777/article/details/4446752

版权

数论基础算法

一、引言

数论曾经被视为数学领域中华而不实的一个分支，然而现在它已经得到了广泛的应用，这其中的部分原因应归结为以大素数为基础的密码体系的建立。这种体系的可行性在于我们可以轻松地找到一些大素数，而体系的安全性则在于将大素数的乘积重新分解因数往往十分困难。本文将介绍数论中比较基本的一些算法，面向的读者应具有代数结构的基础知识。

二、求最大公约数

我们首先要介绍的是能够有效计算两个整数的最大公约数的欧几里得算法（ Euclid's algorithm ，又称辗转相除法）（下文中用 gcd(a,b) 来表示整数 a 与 b 的最大公约数，由于 gcd(a,b)=gcd(|a|,|b|) ，我们将参数的范围限定在非负整数的范围内）。

将 a 、 b 分别分解质因数后即可求得二者的最大公约数，然而即使是最好的算法也不能使位操作意义下的时间复杂度达到多项式级（数论算法的时间复杂度分析常常指在位操作的意义下）。高校的欧几里德算法是建立在下面的定理上的：

对于任何的非负整数 a 和正整数 b, 有 gcd(a,b)=gcd(b,a mod b) 。欧几里得算法的流程为：

EUCLID (a,b)

if b=0

then return a

else return EUCLID (b,a mod b)

比如我们要求 31 和 20 的最大公约数 , 则 EUCLID 的运行过程为：

EUCLID(30,21)=EUCLID(21,9)=EUCLID(9,3)=EUCLID(3,0)=3 。

关于该算法的效率分析，有如下定理阐述：

对任何 k>=1 ，若 a>b>=1 且 b<Fk+1(Fx 为 Fibonacci 数列的第 x 项 ) ，那么 EUCLID(a,b) 的递归调用必小于 k 次。

连续的 Fibonacci 数是 EUCLID 输入的最坏情况，因为对任何 k>2 有 Fk+1 mod Fk=Fk-1 ，所以 gcd(Fk+1,Fk)=gcd(Fk,(Fk+1 mod Fk))=gcd(Fk,Fk-1) ，而 CLID(F3,F2)=EUCLID(2,1) 的计算是需要进行一次递归调用的。

现在我们来改写欧几里得算法以计算附加的有用信息，使得我们可以在算法结束时获得两个整数（可能为 0 或负数） x 、 y 使等式 d=gcd(a,b)=ax+by 成立。如下所示的 EXTENDED-EUCLID 过程接收一对非负整数，返回满足上述等式的三元组 (d,x,y) ：

EXTENDED-EUCLID(a,b)

if b=0

then return (a,1,0)

(d',x',y') ← EXTEND-EUCLID(b,a mod b)

(d,x,y)=(d',y',x'-floor(a/b)*y')

return (d,x,y)

EXTENDED-EUCLID 过程与 EUCLID 过程的时间复杂度在渐近意义下相同，差别仅在于常数因子，而由此得出的 x 、 y 却对解决一些问题十分有用。

三、求解同余线性方程

在一些重要的应用中常会涉及到求解如下同余线性方程的问题： a*x 模 n 余 b ，求 x 。这个问题有解的条件是 d=gcd(a,n)|b ，在已知有解 x0 的情况下，可以推出其在模 n 意义下的全部解 , 即 xi=x0+i*(n/d)(i=0,1...d-1) ，具体过程如下：

MODULAR-LINEAR-EQUATION-SOLVER(a,b,n)

(d,x',y') ← EXTENDED-EUCLID(a,n)

if d|b

then x0 ← x'(b/d) mod n

for i ← 0 to d-1

do print (x0+i(n/d)) mod n

else print "no solutions"

线性同余方程


在数论中，线性同余方程是最基本的同余方程，“线性”表示方程的未知数次数是一次，即形如：


 



的方程。此方程有解当且仅当 b 能够被 a 与 n 的最大公约数整除（记作 gcd(a,n) | b）。这时，如果 x0 是方程的一个解，那么所有的解可以表示为：


 



其中 d 是a 与 n 的最大公约数。在模 n 的完全剩余系 {0,1,…,n-1} 中，恰有 d 个解。





目录 


1 例子 


2 求特殊解 


3 线性同余方程组 


4 参见 





例子


在方程 


3x ≡ 2 (mod 6) 


中， d = gcd(3,6) = 3 ，3 不整除 2，因此方程无解。





在方程 


5x ≡ 2 (mod 6) 


中， d = gcd(5,6) = 1，1 整除 2，因此方程在{0,1,2,3,4,5} 中恰有一个解: x=4。





在方程 


4x ≡ 2 (mod 6) 


中， d = gcd(4,6) = 2，2 整除 2，因此方程在{0,1,2,3,4,5} 中恰有两个解: x=2 and x=5。





求特殊解


对于线性同余方程





ax ≡ b (mod n)      
(1) 


若 d = gcd(a, n 整除 b ，那么为整数。由裴蜀定理，存在整数对 (r,s) （可用辗转相除法求得）使得 ar+sn=d，因此 是方程 (1) 的一个解。其他的解都关于与 x 同余。





举例来说，方程





12x ≡ 20 (mod 28) 


中 d = gcd(12,28) = 4 。注意到 ，因此 是一个解。对模 28 来说，所有的解就是 {4,11,18,25} 。





线性同余方程组


线性同余方程组的求解可以分解为求若干个线性同余方程。比如，对于线性同余方程组：





2x ≡ 2 (mod 6) 


3x ≡ 2 (mod 7) 


2x ≡ 4 (mod 8) 


首先求解第一个方程，得到x ≡ 1 (mod 3)，于是令x = 3k + 1，第二个方程就变为：





9k ≡ 
−
1 (mod 7) 


解得k ≡ 3 (mod 7)。于是，再令k = 7l + 3，第三个方程就可以化为：





42l ≡ 
−
16 (mod 8) 


解出：l ≡ 0 (mod 4)，即 l = 4m。代入原来的表达式就有 x = 21(4m) + 10 = 84m + 10，即解为：





x ≡ 10 (mod 84) 


对于一般情况下是否有解，以及解得情况，则需用到数论中的中国剩余定理。





参见


二次剩余 


中国剩余定理





谈谈解线性同余方程





因为ACM/ICPC中有些题目是关于数论的，特别是解线性同余方程，所以有必要准备下这方面的知识。关于这部分知识，我先后翻看过很多资料，包括陈景润的《初等数论》、程序设计竞赛例题解、“黑书”和很多网上资料，个人认为讲的最好最透彻的是《算法导论》中的有关章节，看了之后恍然大悟。经过几天的自学，自己觉得基本掌握了其中的“奥妙”。拿出来写成文章。





那么什么是线性同余方程？对于方程：ax≡b(mod   
m)，a，b，m都是整数，求解x 的值。





解题例程：pku1061 青蛙的约会 解题报告





符号说明：





                  
mod表示：取模运算





                  
ax≡b(mod   
m)表示：(ax - b) mod m = 0，即同余





                  
gcd(a,b)表示：a和b的最大公约数





求解ax≡b(mod n)的原理：





对于方程ax≡b(mod n)，存在ax + by = gcd(a,b)，x，y是整数。而ax≡b(mod n)的解可以由x，y来堆砌。具体做法，见下面的MLES算法。





第一个问题：求解gcd(a,b)





定理一：gcd(a,b) = gcd(b,a mod b)





实现：古老的欧几里德算法





int Euclid(int a,int b)


{


if(b == 0)


      
return a;


else


      
return Euclid(b,mod(a,b));


}





附：取模运算





int mod(int a,int b)


{


if(a >= 0)


      
return a % b;


else


      
return a % b + b;


}





第二个问题：求解ax + by = gcd(a,b)





定理二：gcd(b,a mod b) = b * x' + (a mod b) * y'


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


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


              
= a * x + b * y





则：x = y'


y = x' - a / b * y'





实现：





triple Extended_Euclid(int a,int b)


{


triple result;


if(b == 0)


{


   
   
result.d = a;


      
result.x = 1;


      
result.y = 0;


}


else


{


      
triple ee = Extended_Euclid(b,mod(a,b));


      
result.d = ee.d;


      
result.x = ee.y;


      
result.y = ee.x - (a/b)*ee.y;


}


return result;


}





附：三元组triple的定义





struct triple


{


int d,x,y;


};





第三个问题：求解ax≡b(mod n)





实现：由x，y堆砌方程的解





int MLES(int a,int b,int n)


{


triple ee = Extended_Euclid(a,n);


if(mod(b,ee.d) == 0)


      
return mod((ee.x * (b / ee.d)),n / ee.d);


else


      
return -1;


}//返回-1为无解，否则返回的是方程的最小解





说明：ax≡b(mod n)解的个数：





如果ee.d 整除 b 则有ee.d个解；





如果ee.d 不能整除 b 则无解。

四、中国剩余定理

中国有一本数学古书《孙子算经》有这样的记载：

「今有物，不知其数，三三数之，剩二，五五数之，剩三，七七数之，剩二，问物几何？」

答曰：「二十三」

术曰：「三三数之剩二，置一百四十，五五数之剩三，置六十三，七七数之剩二，置三十，并之，得二百三十三，以二百一十减之，即得。凡三三数之剩一，则置七十，五五数之剩一，则置二十一，七七数之剩一，则置十五，即得。」

孙子算经的作者及确实着作年代均不可考，不过根据考证，着作年代不会在晋朝之后，以这个考证来说上面这种问题的解法，中国人发现得比西方早，所以这个问题的推广及其解法，被称为中国剩余定理。中国剩余定理（ Chinese Remainder Theorem ）在近代抽象代数学中占有一席非常重要的地位。它揭示了这样两个系统的一致性：一是模两两互质的一组数的同余方程组，二是模它们的乘积的方程。

中国剩余定理的内容如下：

令 n=n1n2...nk ，其中 ni 是两两互质的数，则对 0<=a<n 与 0<=ai<ni 且 ai=a mod ni ， a 与 (a1,a2...,ak) 之间有一种一一对应的关系，一切对 a 的操作均可被等价的转换为对对应 k 元组中的每一元进行同样的操作。因此我们可以将一种表达经过简单的转换后得出另一种表达，其中从 a 到 (a1,a2...,ak) 的转换十分容易，而从 (a1,a2...,ak) 推得对应的 a 则要稍微复杂一些。

首先定义 mi=n/ni(i=1,2...k) ，则 mi 是除了 ni 以外的所有 nj 的乘积，接下来令 ci=mi 与模 n 意义下 mi 的逆元（即 (m_i *m_i ^-1 )mod ni=1 ）的积，则 a 为 (a1c1+a2c2+...+akck) (mod n) 。

例如我们已知 a 模 5 余 2 且模 13 余 3 ，那么 a1=2,n1=m2=5,a2=3,n2=m1=13 ，则有 c1=13*(2 mod 5)=26 ， c2=5*(8 mod 13)=40 ，所以 a=(2*26+3*40) （ mod 65)=42 。

中国剩余定理 ( 同余方程组 ) 小结

问题简单来说就是 a = ai (mod ni) 求未知数 a,

以下小结略去证明 , 只是对定理作了必要的解释 , 要了解相关定理 , 可查阅数论资料 .

中国余数定理 :

设 n=n1*n2...nk, 其中因子两两互质 . 有 : a-----(a1,a2,...,ak), 其中 ai = a mod ni, 则 a 和 (a1,a2,...,ak) 关系是一一对应的 . 就是说可以由 a 求出 (a1,a2,...,ak), 也可以由 (a1,a2,...,ak) 求出 a

推论 1:

对于 a=ai (mod ni) 的同余方程 , 有唯一解

下面说说由 (a1, a2, ..., ak) 求 a 的方法 :

定义 mi = n1*n2*...nk / ni; ci = mi(mf mod ni); 其中 mi*mf mod ni = 1;

则 a = (a1*c1+a2*c2+...+ak*ck) (mod n) ( 注 : 由此等式可求 a%n, 当 n 很大时 )

中国剩余定理关键是 mf 的求法 , 如果理解了扩展欧几里得 ax+by=d, 就可以想到 :

mi*mf mod ni = 1 => mi*mf+ni*y=1;

代码如下 :

#include <iostream>

#include <cmath>

using namespace std;

const int MAXN = 100;

int nn, a[MAXN], n[MAXN];

int egcd(int a, int b, int &x, int &y) {

int d;

if (b == 0) {

x = 1; y = 0;

return a;

} else {

d = egcd(b, a%b, y, x);

y -= a/b*x;

return d;

}

int lmes() {

int i, tm=1, mf, y, ret=0, m;

for (i=0; i<nn; i++) tm *= n[i];

for (i=0; i<nn; i++) {

m = tm/n[i];

egcd(m, n[i], mf, y);

ret += (a[i]*m*(mf%n[i]))%tm;

}

return (ret+tm)%tm;

}

int main() {

a[0] = 4; a[1] = 5;

n[0] = 5; n[1] = 11;

nn = 2;

printf("%d/n", lmes());

return 0;

}

五、模取幂运算

m^e mod n 叫做模取幂运算，事实上， m^e mod n 可以直接计算，没有必要先算 m^e ，根据简单的数论知识，很容易设计一个分治算法。具体如下 :

设 <b[k], b[k-1],...,b[1],b[0]> 是整数 b 的二进制表示（即 b 的二进制有 k+1 位， b[k] 是最高位），下列过程随着 c 的值从 0 到 b 成倍增加，最终计算出 a^c mod n

Modular-Exponentiation(a, b, n)

　 c ← 0

d ← 1

　设 <b[k],b[k-1],..b[0]> 是 b 的二进制表示

　 for i ← k downto 0

　 do c ← 2c

d ← (d*d) mod n

if b[i] = 1

then c ← c + 1

d ← (d*a) mod n

　 return d

上述伪代码依次计算出的每个幂或者是前一个幂的两倍，或者是前一个幂的两倍再乘上底数 a 。过程依次从右到左逐个读入 b 的二进制表示已控制执行哪一种操作。循环中的每次迭代都用到了下面的两个恒等式中的一个：

a^(2c) mod n = (a^c mod n)^2

a^(2c+1) mod n = a * (a^c mod n)^2

用哪一个恒等式取决于 b[i]=0 还是 1 。由于平方在每次迭代中起着关键作用，所以这种方法叫做“反复平方法 (repeated squaring) ”。在读入 b[i] 位并进行相应处理后， c 的值与 b 的二进制表示 <b[k],b[k-1],..b[0]> 的前缀的值相同。事实上，算法中并不真正需要变量 c ，只是为了说明算法才设置了变量 c ：当 c 成倍增加时，算法保持条件 d = a^c mod n 不变，直至 c=b 。