数论相关知识

ACM 进阶学习第一课----同余相关之欧几里得算法及其扩展（2）

分类： C/C++学习进阶之路 2013-11-03 02:16  747人阅读  评论(0)  收藏  举报

目录(?)[+]

最大公约数算法分析

欧几里德算法

伪代码

while b>0
    do r←a%b
         a←b
         b←r
return a

算法分析：

欧几里德算法是计算最大公约数的传统算法,也是最简单的算法,效率很高
时间复杂度:O(lgn)(最坏情况:斐波那契数列相邻的两项)
空间复杂度:O(1)
但是,对于大整数来说,取模运算非常耗时

Stein算法

伪代码：

Stein算法(假设0<=b<a):
r←0
while b>0
   do if a偶,b偶 then a←a>>1   b←b>>1   r←r+1
        else if a偶,b奇 then a←a>>1
        else if a奇,b偶 then b←b>>1
        else if a奇,b奇 then a←(a-b)>>1
        if a<b then 交换a,b
return a<<r

算法分析

渐近时间,空间复杂度均与欧几里德算法相同
原理:gcd(ka,kb)=k*gcd(a,b)
最大特点:只有移位和加减法计算,避免了大整数的取模运算

代码：

[cpp]  view plain copy print ?

1. unsigned MaxDivisor(unsigned a, unsigned b)   
2. {   
3.     unsigned c = 0;   
4.     while(1)  
5.     {   
6.     // 退出条件   
7.         if(a==0)   
8.             return b << c;  
9.         else if(b == 0)   
10.             return a << c;  
11.     // 为提高速度，采用位的与运算，避免用取模判断奇偶   
12.         if(!(a & 1) && !(b & 1)) //a,b 都是偶数   
13.         {   
14.             a >>= 1; b >>= 1; ++c;   
15.         }   
16.         else if(!(a & 1) && (b & 1)) //a偶 b奇   
17.         {   
18.             a >>= 1;   
19.         }   
20.         else if((a & 1) && !(b & 1)) //a奇 b偶   
21.         {  
22.              b >>= 1;   
23.         }   
24.         else if((a & 1) && (b & 1)) //a,b都是奇数   
25.         {   
26.             unsigned tmp = a>b?b:a; //取较小的一个   
27.             a = a>b?a-b:(b-a); //绝对差值  
28.             b = tmp;   
29.         }   
30.     }  
31. }   

扩展欧几里德算法

基本算法

对于不完全为 0 的非负整数 a，b，gcd（a，b）表示 a，b 的最大公约数，必然存在整数对 x，y ，使得 gcd（a，b）=ax+by。

【证明】

设 a>b

1，显然当 b=0，gcd（a，b）=a。此时 x=1，y=0；
2，ab!=0 时
设 :ax1+by1=gcd(a,b);
显然也有：bx2+(a mod b)y2=gcd(b,a mod b);
根据朴素的欧几里德 原理有 gcd(a,b)=gcd(b,a mod b);
则:ax1+by1=bx2+(a mod b)y2;

即:ax1+by1=bx2+(a-(a/b)*b)y2=ay2+bx2-(a/b)*by2;

根据恒等定理得：x1=y2; y1=x2-(a/b)*y2;

这样我们就得到了求解 x1,y1 的方法：x1，y1 的值基于 x2，y2.

上面的思想是以递归定义的，因为 gcd 不断的递归求解一定会有个时候 b=0，所以递归可以结束。

扩展欧几里得递归代码：

[cpp]  view plain copy print ?

1. int exgcd(int a,int b,int &x,int &y)  
2. {  
3.     if(b==0)  
4.     {  
5.         x=1;  
6.         y=0;  
7.         return a;  
8.     }  
9.     int r=exgcd(b,a%b,x,y);  
10.     int t=x;  
11.     x=y;  
12.     y=t-a/b*y;  
13.     return r;  
14. }  

非递归代码：

[cpp]  view plain copy print ?

1. int exgcd(int m,int n,int &x,int &y)  
2. {  
3.     int x1,y1,x0,y0;  
4.     x0=1; y0=0;  
5.     x1=0; y1=1;  
6.     x=0; y=1;  
7.     int r=m%n;  
8.     int q=(m-r)/n;  
9.     while(r)  
10.     {  
11.         x=x0-q*x1; y=y0-q*y1;  
12.         x0=x1; y0=y1;  
13.         x1=x; y1=y;  
14.         m=n; n=r; r=m%n;  
15.         q=(m-r)/n;  
16.     }  
17.     return n;  
18. }  

应用

扩展欧几里德算法的应用主要有以下三方面：

（1）求解不定方程；

（2）求解模线性方程（线性同余方程）；

（3）求解模的逆元；

（1）使用扩展欧几里德算法解决不定方程的办法：

  对于不定整数方程pa+qb=c，若 c mod Gcd(p, q)=0,则该方程存在整数解，否则不存在整数解。
  上面已经列出找一个整数解的方法，在找到p * a+q * b = Gcd(p, q)的一组解p0,q0后，p * a+q * b = Gcd(p, q)的其他整数解满足：
  p = p0 + b/Gcd(p, q) * t 
  q = q0 - a/Gcd(p, q) * t(其中t为任意整数)
  至于pa+qb=c的整数解，只需将p * a+q * b = Gcd(p, q)的每个解乘上 c/Gcd(p, q) 即可。

  在找到p * a+q * b = Gcd(a, b)的一组解p0,q0后，应该是得到p * a+q * b = c的一组解p1 = p0*(c/Gcd(a,b)),q1 = q0*(c/Gcd(a,b))，

  p * a+q * b = c的其他整数解满足：

   p = p1 + b/Gcd(a, b) * t

  q = q1 - a/Gcd(a, b) * t(其中t为任意整数)

  p 、q就是p * a+q * b = c的所有整数解。

相关证明可参考： http://www.cnblogs.com/void/archive/2011/04/18/2020357.html

 

用扩展欧几里得算法解不定方程ax+by=c;

代码如下：

[cpp]  view plain copy print ?

1. bool linear_equation(int a,int b,int c,int &x,int &y)  
2. {  
3.     int d=exgcd(a,b,x,y);  
4.     if(c%d)  
5.         return false;  
6.     int k=c/d;  
7.     x*=k; y*=k;    //求得的只是其中一组解  
8.     return true;  
9. }  

(2)用扩展欧几里德算法求解模线性方程的方法：

    同余方程 ax≡b (mod n)对于未知数 x 有解，当且仅当 gcd(a,n) | b。且方程有解时，方程有 gcd(a,n) 个解。

    求解方程 ax≡b (mod n) 相当于求解方程 ax+ ny= b, (x, y为整数)

    设 d= gcd(a,n)，假如整数 x 和 y，满足 d= ax+ ny(用扩展欧几里德得出)。如果 d| b，则方程

    a* x0+ n* y0= d， 方程两边乘以 b/ d，(因为 d|b，所以能够整除)，得到 a* x0* b/ d+ n* y0* b/ d= b。
    所以 x= x0* b/ d，y= y0* b/ d 为 ax+ ny= b 的一个解，所以 x= x0* b/ d 为 ax= b (mod n ) 的解。

    ax≡b (mod n)的一个解为 x0= x* (b/ d ) mod n，且方程的 d 个解分别为 xi= (x0+ i* (n/ d ))mod n {i= 0... d-1}。

    设ans=x*(b/d),s=n/d;

    方程ax≡b (mod n)的最小整数解为：(ans%s+s)%s;

    相关证明：

    证明方程有一解是: x0 = x'(b/d) mod n;
    由 a*x0 = a*x'(b/d) (mod n)
         a*x0 = d (b/d) (mod n)   (由于 ax' = d (mod n))
                 = b (mod n)

    证明方程有d个解: xi = x0 + i*(n/d)  (mod n);
    由 a*xi (mod n) = a * (x0 + i*(n/d)) (mod n)
                             = (a*x0+a*i*(n/d)) (mod n)
                             = a * x0 (mod n)             (由于 d | a)
                             = b

     

首先看一个简单的例子：

=4(mod3)

解得x = 2,5,8,11,14.......

由此可以发现一个规律，就是解的间隔是3.

那么这个解的间隔是怎么决定的呢？

如果可以设法找到第一个解，并且求出解之间的间隔，那么就可以求出模的线性方程的解集了.

我们设解之间的间隔为dx.

那么有

a*x = b(mod n);

a*(x+dx) = b(mod n);

两式相减，得到:

a*dx(mod n)= 0;

也就是说a*dx就是a的倍数，同时也是n的倍数，即a*dx是a 和 n的公倍数.为了求出dx,我们应该求出a 和 n的最小公倍数,此时对应的dx是最小的.

设a 和 n的最大公约数为d,那么a 和 n 的最小公倍数为(a*n)/d.

即a*dx = a*n/d;

所以dx = n/d.

因此解之间的间隔就求出来了.

    代码如下：

[cpp]  view plain copy print ?

1. bool modular_linear_equation(int a,int b,int n)  
2. {  
3.     int x,y,x0,i;  
4.     int d=exgcd(a,n,x,y);  
5.     if(b%d)  
6.         return false;  
7.     x0=x*(b/d)%n;   //特解  
8.     for(i=1;i<d;i++)  
9.         printf("%d\n",(x0+i*(n/d))%n);  
10.     return true;  
11. }  

（3）用欧几里德算法求模的逆元：

       同余方程ax≡b (mod n)，如果 gcd(a,n)== 1，则方程只有唯一解。

      在这种情况下，如果 b== 1，同余方程就是 ax=1 (mod n ),gcd(a,n)= 1。

      这时称求出的 x 为 a 的对模 n 乘法的逆元。

      对于同余方程 ax= 1(mod n )， gcd(a,n)= 1 的求解就是求解方程

      ax+ ny= 1，x, y 为整数。这个可用扩展欧几里德算法求出，原同余方程的唯一解就是用扩展欧几里德算法得出的 x 。

ACM 进阶学习第一课----同余相关之中国剩余定理

分类： 算法总结 2013-11-04 01:05  10983人阅读  评论(30)  收藏  举报

目录(?)[+]

问题引入

"物不知数"问题：

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

答曰:'二十三.'

术曰: 三三数之剩二,置一百四十,五五数之剩三,置六十三,七七数之剩二,置三十,并之,得二百三十三,以二百一十减之,即得.凡三三数之剩一,则置七十,五五数之剩一,则置二十一,七七数之剩一,则置十五,即得." 

 --孙子算经

当代解释

令任意固定整数为M，当M/A余a，M/B余b，M/C余c，M/D余d，…，M/Z余z时，这里的A，B，C，D，…，Z为除数，除数为任意自然数（如果为0，没有任何意义，如果为1，在孙子定理中没有计算和探讨的价值，所以，不包括0和1）时；余数a，b，c，d，z为自然整数时。
1、当命题正确时，在这些除数的最小公倍数内有解，有唯一的解，每一个最小公倍数内都有唯一的解；当命题错误时，在整个自然数范围内都无解。
2、当M在两个或两个以上的除数的最小公倍数内时，这两个或两个以上的除数和余数可以定位M在最小公倍数内的具体位置，也就是M的大小。
3、正确的命题，指没有矛盾的命题：分别除以A，B，C，D，…，Z不同的余数组合个数=A，B，C，D，…，Z的最小公倍数=不同的余数组合的循环周期．

问题分析

一、

" 物不知数"问题的抽象表示:
x≡2(mod 3)
x≡3(mod 5)
x≡2(mod 7)
求满足上述条件的最小正整数x

二、

" 物不知数"解法的数学表示:
任取被3除余2的5和7的倍数:140
任取被5除余3的7和3的倍数:63
任取被7除余2的3和5的倍数:30
140+63+30=233
减去3,5,7的公倍数中不超过233的最大的数210得到答案 23

三、问题归化

一般性问题: 给定两两互质的正整数n1,n2,...,nk,要求找到最小的正整数a,满足a≡ai(mod ni)
将问题分解成若干次求解二元模线性方程组的解

算法步骤

令n=n1n2···nk，mi=n/ni
利用扩展欧几里德算法计算出xi满足mixi ≡ 1(mod ni),由于n1,n2,...,nk两两互质,必有gcd(mi,ni)=1,即可保证一定有解
则a≡a1x1m1 + a2x2m2 + ... + akxkmk (mod n)

典型应用

大整数的表示
选取两两互素的正整数n1,n2,...,nk
已知对每个ni取模的值ri,就可以唯一确定一个1~n1n2...nk的大整数
做大整数加,减,乘法时,只要保证在这个范围内,均可转化为分别对相应的余数进行计算

模板代码

[cpp]  view plain copy print ?

1. <pre name="code" class="cpp">/*****  ACM之中国剩余定理 ********/  
2.   
3. /******** written by C_Shit_Hu ************/  
4.   
5. 扩展欧几里得算法的运用///  
6.   
7. /****************************************************************************/  
8. /*  
9. 由于VC下面无法使用cout或者cin输出64位的整数。 
10. 故改用printf. 
11. */  
12. /****************************************************************************/  
13.   
14. #include<iostream>  
15. using namespace std;  
16. typedef _int64 llong;  
17. llong b[1000],w[1000];  
18.   
19.   
20. // 扩展欧几里得算法  
21. // 递归的形式  
22. llong extended_euclid(llong a, llong b, llong &x, llong  &y)   
23. {  
24.     llong d;  
25.     if(b == 0)   
26.     {x = 1; y = 0; return a;}  
27.       
28.     d = extended_euclid(b, a % b, y, x);  
29.     y -= a / b * x;  
30.     return d;  
31. }  
32.   
33. // 中国剩余定理  
34. llong chinese_remainder(int len)   
35. {   
36.     llong i, d, x, y, m, n, ret;  
37.     ret = 0; n = 1;   
38.     for(i=0; i < len ;i++) n *= w[i];  
39.     for(i=0; i < len ;i++)   
40.     {   
41.         m = n / w[i];  
42.         d = extended_euclid(w[i], m, x, y);  
43.         ret = (ret + y*m*b[i]) % n;  
44.     }  
45.     return (n + ret%n) % n;  
46. }  
47.   
48. int main()  
49. {  
50.     int  n,i;  
51.     llong res;  
52.     // 输入测试的除数和余数的组数  
53.     cout << "输入测试的除数和余数的组数: " ;  
54.     while (scanf("%d",&n)!=EOF)  
55.     {  
56.         // 输入除数和余数  
57.         cout << "输入余数和除数:" << endl;  
58.         for(i=0;i<n;i++)  scanf("%I64d%I64d",&b[i],&w[i]);  
59.         res=chinese_remainder(n);  
60.         cout << "结果为：" ;  
61.         printf("%I64d\n",res);  
62.         cout << "输入测试的除数和余数的组数: " ;  
63.     }  
64.     return 0;  
65. }</pre><br><br>  

运行结果：

【未完待续】。。。。

 

ACM 进阶学习第一课----简单数学问题之同余相关（1）

分类： 算法总结 2013-11-03 01:05  1290人阅读  评论(44)  收藏  举报

目录(?)[+]

前言

    在ACM竞赛中,经常可以看到数学问题的身影，可以是纯数学问题,也可以是需要利用数学上的一些公式,定理,算法来辅助解决的问题。会者不难,而不会的选手在赛场上一般很难推出公式或进行证明，往往想起来费劲,写起来却很轻松。

常见的数学问题：

数论
组合数学
博弈论
线性代数
高等数学
线性规划
概率统计
...

关于数论

简而言之,数论就是研究整数的理论，在ACM竞赛中,经常用到数论的相关知识。纯数论的题目不多,大部分是和其他类型的问题结合起来的。约数,倍数,模线性方程,欧拉定理,素数。

数论的主要内容

第一部分:同余相关
整除的性质->欧几里德算法
->扩展欧几里德算法->中国剩余定理
第二部分:素数相关
算术基本定理->欧拉定理
->素数测试-> Pollard rho方法

基本概念：

1、素数合数

如果大于1的正整数p仅有的正因子是1和p, 则称p为素数(prime)。
大于1又不是素数的正整数称为合数(compound)，如果n是合数, 则n必有一个小于或等于n1/2的素因子。

2、算数基本定理

·····每个正整数都可以惟一地表示成素数的乘积，其中素数因子从小到大依次出现（这里的“乘积”可以有0个、1个或多个素因子）。
·····换句话说, 任意正整数n可以写成n=2a1*3a2*5a3*…,其中a1,a2,a3等为非负整数
·····这个定理也叫做惟一分解定理。它是一个定理而不是公理！虽然在大多人看来，它是“显然成立”的，但它的确是需要证明的定理

3、除法和同余

---令a为整数，d为正整数，那么有惟一的整数q和r，其中0≤r<d，使得a=dq+r
---可以用这个定理来定义除法：d叫除数，a叫被除数，q叫商，r叫余数。如果两个数a,b除以一个数c的余数相等，说a和b关于模c同余，记作a≡b(mod c) 

本次先学习同余相关。

同余相关

同余相关主要内容

整除的性质
欧几里德算法
扩展欧几里德算法
中国剩余定理

整除的性质

• 若a|b, a|c, 则a|(b+c)
• 若a|b, 那么对所有整数c, a|bc
• 若a|b, b|c, 则a|c
• 若a|b,b|a => a=±b
  
• 若a=kb±c => a,b的公因数与b,c的公因数完全相同
• 整除关系具有传递性.
• 整除显然有自反性和反对称性，所以它是一个偏序关系。(partial order), <|,Z>是一个格

FIR：最大公约数·最小公倍数

最大公约数

令a和b是不全为0的两个整数，能使d|a和d|b的最大整数称为a和b的最大公约数 (greatest common divisor) ，用gcd(a,b)表示，或者记为(a,b)。

显然满足以下性质：
·gcd(a,b)=gcd(b,a)=gcd(-a,b)=gcd(|a|,|b|)
·gcd(a,0)=|a|
·gcd(a,ka)=|a|,k为任意整数

·a,b互素,则gcd(a,b)=1

最小公倍数lcm
令a和b是不全为0的两个整数，能使a|d和b|d的最小整数称为a和b的最小公倍数，用lcm(a,b)表示，或者记为[a,b] 

则有定理：

ab=gcd(a,b)*lcm(a,b)

即是：最小公倍数=两数的乘积/最大公约（因）数

证明如图：

明显已经证明。

实例测试

【1】11年腾讯笔试题目

给出一个数N，含数字1、2、3、4，把N的所有数字重新排列一下组成一个新数，使它是7的倍数。

【分析】

把数字1、2、3、4从中抽出，然后把其他数字按照原顺序排列（事实上，怎么排列都无所谓）组成自然数w
w*10,000整除7取余有7种可能，即是为0、1、2、3、4、5、6。这时如果能用数字1、2、3、4排列出7个数，使它们整除7取余的值分别为0、1、2、3、4、5、6，把这个4位数接在w后面即为问题的解。

幸运的是，有这样的7个数，如下：

余数	0（7）	1	2	3	4	5	6
排列	3241	1324	1234	2341	1243	1342	2134

    

选取某一个排列作为后缀时，若w×10000模7余d，则选取（7-d）为余数的那个排列即是W可用表达式：W=d + 7 * m, 而后四位的排列则可用L=（7-d） + 7 * n,即是L+W = 7 + 7*（m+n），即是7的倍数 。

参考代码如下：

[cpp]  view plain copy print ?

1. /*****  简单数论题目 ********/  
2.   
3. /******** written by C_Shit_Hu ************/  
4.   
5. 数论///  
6.   
7. /****************************************************************************/  
8. /*  
9.  把数字1，2，3，4从中各抽出1个，然后把其他数字按原顺序（其实任一顺序都可以）排列，组成自然数w。 
10.  w×10000模7有7种可能，即是0，1，2，3，4，5，6，这时若能用数字1，2，3，4排列出7个数， 
11.  使它们整除7取余的值分别为0，1，2，3，4，5，6。则把这个4位数接在w后面即为原问题的解。 
12.  
13. */  
14. /****************************************************************************/  
15.   
16. #include <stdlib.h>  
17. #include <stdio.h>  
18. #include <iostream.h>  
19.   
20. int main()  
21. {  
22.     int data;  
23.     scanf("%d",&data);//读取输入数据  
24.     /**找出1,2,3,4组合中对7取余数分别为0,6,5,4,3,2,1,情况*/  
25.     int add[7]={3241,2134,2413,1243,2341,1234,1324};  
26.     int temp1=0;  
27.     int temp2=0;  
28.     bool flat[4]={false,false,false,false};//用于标记抽取到数字1,2,3,4的情况，即是验证输入数据是否符合要求  
29.     /*以下循环把数字1,2,3，4从中抽取出来，然后其他数字按照原输入数据的逆序排成自然数*/  
30.     while(data!=0)  
31.     {  
32.         temp1 = data%10;//抽取末位数字  
33.         data = data/10;//去掉末位数字，向右移  
34.           
35.         if((temp1==1 || temp1==2 || temp1==3 || temp1==4)&&(!flat[temp1-1]))  
36.         {  
37.             flat[temp1-1]=true;  
38.         }  
39.         else  
40.         {  
41.             temp2 =temp2*10 + temp1;//如果不是数字1234则按照原输入数据的逆序排成自然数  
42.         }  
43.     }  
44.     if(flat[0]&&flat[1]&&flat[2]&&flat[3])//判断输入数据是否合理  
45.     {  
46.         data=temp2*10000;        //为后面的加运算腾出末四位数  
47.         temp2=data%7;            //求余            
48.         if(temp2==0) //根据求出的余数加上由1,2，3,4组成的四位数  
49.         {  
50.             data = data+3241;  
51.         }  
52.         else if(temp2 ==6)  
53.         {  
54.             data = data+2134;  
55.         }  
56.         else if(temp2 ==5)  
57.         {  
58.             data = data+2413;  
59.         }  
60.         else if(temp2 ==4)  
61.         {  
62.             data = data+1243;  
63.         }  
64.         else if(temp2 ==3)  
65.         {  
66.             data = data+2341;  
67.         }  
68.         else if(temp2 ==2)  
69.         {  
70.             data = data+1234;  
71.         }  
72.         else if(temp2 ==1)  
73.         {  
74.             data = data+1324;  
75.         }  
76.         printf("%d\n",data);//输出符合要求的结果  
77.     }  
78.     else  
79.     {  
80.         printf("the input is illegal\n");  
81.     }  
82.       
83.     return 0;  
84. }  
85.   
86.   
87.   
88.   
89. /********************  心得体会  **********************/  
90. /* 
91. 分析数据的关键部分在拆分数字部分 
92. */  

运行结果：

【未完待续】。。。。