定理一：（除法定理）對於任何整數a和任何正整數n，存在唯一的整數q和r，滿足0<=r<n且a = qn + r

證明略

定理二：gcd(x,y) = d，當ax + by為正整數時有 min(ax + by) = d;

顯然當gcd(x,y) = 1 時。就是ax + by = 1

證明：                   

設集合S為x和y的線性組合ax+by的集合，其中a,b為整數                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         

由0 < x^2 + y^2可知S中有正整數，由最小自然數原理可知，S中存在最小正正數，設為s0，那麼首先d|s0，所以d<=s0

第二，對於x，和y 的其中任意一個，如果除以s0餘數是r1，r2。那麼可知r1,r2也屬於S.所以r1 = r2 = 0.即s0能整除x和y，所以s0為x,y的最大公約數。

定理三：lcm(a,b) = ab/gcd(a,b) 證明略

定理四：假設gcd(n,m) = d，那麼a ≡ b (mod n)且 a≡b( mod m)等價於 a≡b mod( lcm(n,m));

證明略

定理五： a ≡ b (mod nm) 推出  a ≡ b (mod m) ,a ≡

一、擴充套件歐幾里得演算法

擴充套件歐幾里德演算法是用來在已知a, b求解一組x，y使得

ax+by = Gcd(a, b) =d（兩邊除以d，得到a'x+b'y=1,a'與b'是互質的，根據上面的定理，肯定是有正數解的）。擴充套件歐幾里德常用在求解模線性方程及方程組中。

現在求解 ax+by=1

假設a>b，a=kb+l（l<b）則上式可以寫成（k * b+l）x+by = l * x+b * (k * x+y) = 1,令y' = k * x+y,則方程為l * x + b * y' = 1

則此方程的解x為原方程的解，y = (y' - k*x)

則上述可以證明擴充套件歐幾里得演算法的正確性

給出程式碼：

void exgcd(int a, int b, int& d, int& x, int& y)  
{
   if(!b) { d = a; x = 1; y = 0; }  
   else  { exgcd(b, a % b, d, y, x); y -= x * (a / b); }  
}  

這樣就可以求出一組解x,y

假設另外一組解為x1,y1

那麼a(x-x1)=b(y-y1)

兩邊除以gcd(a,b)

則，a'(x-x1) = b'(y-y1),其中 a'=a/gcd(a.b)，b'=b/gcd(a.b)

由於a',b'互素，x - x1 是b'的倍數，設為kb‘得（y - y1）=ka'可以 求出所有的x和y

關於ax+by=c 當c不是gcd(a,b)時無解

由上面當c為為k*gcd(a,b)時，一組解為（kx,ky）(x,y)為 ax+by=gcd(a,b)的一組解

當c不是gcd(a,b)的倍數時，顯然無解，因為左邊是最小公約數的 倍數，右邊不是

擴充套件：歐幾里得演算法可以求多於2元的模線性方程，顯然不斷迭代就可以解出來了。

擴充套件歐幾里得演算法通常解決如下三種問題

（1）求解不定方程；

（2）求解模線性方程（線性同餘方程）；

（3）求解模的逆元；

另外兩個問題可以轉化成我們熟悉的方式（1）求解

求解逆元：

 模P乘法逆元

 對於整數a、p，如果存在整數b，滿足a*b mod p =1，則說，b是a的模p乘法逆元。

 定理：a存在模p的乘法逆元的充要條件是gcd(a,p) = 1

 證明：

 充分性：

 如果gcd(a,p) = 1，根據尤拉定理，a^φ(p) ≡ 1 mod p，因此

 顯然aφ(p) - 1 mod p是a的模p乘法逆元。

 必要性：

 假設存在 a mod p的乘法逆元為b

 ab ≡ 1 mod p

 則ab = kp +1 ，所以1 = ab - kp

 因為gcd(a,p) = d

 所以d | 1

 所以d只能為1 

可以由1 = ab - kp求出b和k

計算a模p的乘法逆元

下面程式碼算出a,和b的最大公約數和a 模b的逆元ar,以及b模a的逆元br（用非遞迴的計算）

int gcd(int a,int b,int &ar,int &br){
    int x1,x2,x3;
    int y1,y2,y3;
    int t1,t2,t3;
    if(0 == a){//有一個數為0,就不存在乘法逆元
        ar = 0;br = 0 ;
        return b;
    }
    if(0 == b){
        ar = 0;
        br = 0 ;
        return a;
    }
    x1 = 1; x2 = 0; x3 = a; y1 = 0; y2 = 1; y3 = b; int k;
    for( t3 = x3 % y3 ; t3 != 0 ; t3 = x3 % y3){
        k = x3 / y3;
        t2 = x2 - k * y2;t1 = x1 - k * y1;
        x1 = y1;
        x2 = y2;
        x3 = y3;
        y1 = t1;
        y2 = t2;
        y3 = t3;
    }
    if( y3 == 1){//有乘法逆元
        ar = y2; br = x1;
        return 1;
    }
    else{//公約數不為1,無乘法逆元
        ar = 0;br = 0;return y3;
    }
}
 

遞迴計算：

bigint inv(bigint a, bigint n) {
    bigint d, x, y;
    Gcd(a, n, d, x, y);
    if (d == 1) return ( x % n + n) % n;
    else return -1;
}
 

二、中國剩餘定理

有一個數x，除以n個模數a[0],a[1]……a[n-1]（所有模數 兩兩互質）後，得到n個餘數r[0],r[1]……r[n-1]，求x。

先看三個方程的情況

X mod 3 = 2

X mod 5 = 3

X mod 7 = 2

構造法求解：

我們可以求k1,k2,k3

使得： k1 mod 3 = 2, k2 mod 3 = k3 mod 3 = 0

                  k2 mod 5 = 3, k1 mod 5 = k1 mod 5 = 0

    K3 mod 7 = 2,k1 mod 7 = k2 mod 7 = 0

每個k都由獨立的三個方程確定，就可以解出 k1,k2,k3了。

解k1,k1 = k*35。我們可以先解出k。

由 k * 35 mod 3 = 2可以得到不定方程

K*35 - 2 = 3*y，即：35k - 3y = 2，這個可以由擴充套件歐幾里得演算法求出。

K1,k2也可以按照同樣的方法算出來。最後k1+k2+k3，就是一個答案。  

一般的求解過程：

1）令p = a[0]*a[1]*……*a[n-1]，k[i]=p/a[i]

2） 我們要找到這樣的數  d[i]%k[i]==0且

    d[i]%a[i]==r[i]（i 從0到n-1）(這裡可以用

    到擴充套件歐幾里德演算法來求得d[i])

3）當求出d[0],d[1],d[2]……d[n-1]後，相加得  

     到w，w為其中一個解，且通解x = w + t*p  (t  

     為任意自然數)

中國剩餘定理一個條件就是模數要兩兩互素（條件非常苛刻），怎麼求解一般的線性模方程,也就是模數不滿足兩兩互素的情況。數論裡的一般解法是分解因子，但是程式設計難度比較大，時間複雜度也比較高，這裡採用的解法是將所有方程合併成一個方程求解。

假設有 x mod m1 = a

 x mod m2 = b

 Gcd(m1,m2) = d

 m1 = m1’d

 m2 = m2’d

那麼，x = k1m1 + a ①

x = k2m2 + b 

由k1m1 + a = k2m2 + b，得 k1m1’ - k2m2’ = (b - a)/d ② .所以有解的情況肯定有 d|(b - a)  （這是充要條件，是判斷有沒有解的依據）。② 式可以用擴充套件歐幾里得演算法求出k1的一個解，我們先設為k那麼k1的通解就是k + tm2’.將其代入①,得到x = t*lcm(m1,m2) + km1 + a,即 x ≡ km1 + a (mod  lcm(m1,m2))這樣就完成了兩個式子的合併。按照此方法將所有方程合併成一個求解即可。複雜度是O(nlog(max{mi})

練習：poj 2891 

三、埃拉託斯特尼篩法

// 1：這是最原始的篩法，還有待優化 

#define Max 1000000
bool prime[Max];
void IsPrime(){
int t=(int)sqrt(Max*1.0);
     prime[0]=prime[1]=0;prime[2]=1;
     for(int i=3;i<ax;i++)
        prime[i] = i&1;
     for(int i=3;i<=t;i++)
       if(prime[i])
         for(int j= (i << 1);j<Max;j+=i)
            prime[j]=0;
}

優化1：

#include<iostream>
using namespace std;    
const long N = 200000;   
long prime[N] = {0},num_prime = 0;    
int isNotPrime[N] = {1, 1};   
int main()    
{     
     	for(long i = 2 ; i < N ; i ++)       
       	{            
if(! isNotPrime[i])               
 	prime[num_prime ++]=i;  
//關鍵處1        
for(long j = 0 ; j < num_prime && i * prime[j] <  N ; j ++)
    	{               
      	isNotPrime[i * prime[j]] = 1;  
  	if( !(i % prime[j] ) )  //關鍵處2                  
break;           
}        
}        
return 0;   
}  

首先，先明確一個條件，任何合數都能表示成一系列素數的積。

不管 i 是否是素數，都會執行到“關鍵處1”，

①如果 i 都是是素數的話，那簡單，一個大的素數 i 乘以不大於 i 的素數，這樣篩除的數跟之前的是不會重複的。篩出的數都是 N=p1*p2的形式, p1，p2之間不相等

②如果 i 是合數，此時 i 可以表示成遞增素數相乘 i=p1*p2*...*pn, pi都是素數（2<=i<=n），  pi<=pj  ( i<=j )

p1是最小的係數。

根據“關鍵處2”的定義，當p1==prime[j] 的時候，篩除就終止了，也就是說，只能篩出不大於p1的質數*i。

證明略

優化2：

 易於理解。 只算奇數部分，時空效率都還不錯！

half=SIZE/2; 
int sn = (int) sqrt(SIZE); 
for (i = 0; i < half; i++) 
   p[i] = true;// 初始化全部奇數為素數。p[0]對應3，即p[i]對應2*i+3 
for (i = 0; i < sn; i++) {    
if(p[i])//如果 i+i+3 是素數
{     
    for(k=i+i+3, j=k*i+k+i; j < half; j+=k) 
    // 篩法起點是 p[i]所對應素數的平方 k^2                                        
    // k^2在 p 中的位置是 k*i+k+i
    //    下標 i         k*i+k+i
    //對應數值 k=i+i+3   k^2         
       p[j]=false; 
} 
} 
//素數都存放在 p 陣列中，p[i]=true代表 i+i+2 是素數。
//舉例，3是素數，按3*3,3*5,3*7...的次序篩選，因為只儲存奇數，所以不用刪3*4，3*6....

原始的埃拉託斯特尼篩法的漸進時間複雜度是nlnlnn，可以當做O(n)，通常情況下，它的效率已經足夠。改進後的時間複雜度是均攤是O(n)。

練習：hdu 1431 poj 2689

四、尤拉函式 和尤拉定理

在數論，對正整數n，尤拉函式是少於或等於n的數中與n互質的數的數目，它又稱為Euler's totient function、φ函式。

積性函式：如果gcd(n,m) = 1,有f(nm) = f(n)f(m)，則稱函式f是積性的。

φ是積性的。

證明：如果(m,m’) = 1,那麼當a和a’分別取遍模m和模m’的完全剩餘系時，a’m + am’取遍mm’的一個完全剩餘系。又因為

(a’m + am’,mm’) = 1 <------> (a’m +am’,m) = 1,(a’m+am’,m’) = 1

<------>(am’,m) = 1,(a’m,m’) = 1 <-----> (a,m) = 1,(a’,m’) = 1

從而這φ個小於mm’且與mm’互素的數是這φ（m）φ（m'）個數a’m + am’的最小正剩餘，其中a與m互素，而a’與m’互素，從而有φ（mm’）=φ(m)φ(m’)

從而可以推出，其中p是素數，進而推出尤拉公式：。

尤拉函式的性質：

若(N%a==0 && (N/a)%a==0) 則有:φ(N)=φ(N/a)*a;

若(N%a==0 && (N/a)%a!=0) 則有:φ(N)=φ(N/a)*(a-1);

當n為奇數時有

 a^b ≡a^(b%phi(c) + phi(c)) (mod c) （b >= phi(c)） 

求單個尤拉函式值：

Int  euler(int x)
{
int i, res=x;
int m = (int)sqrt(x*1.0) + 1;
for (i = 2; i < m; i++)
if(x%i == 0) {
res = res / i * (i - 1);	//先進行除法是為了防止中間資料的溢位  
while (x % i == 0) x /= i; // 保證i一定是素數
}
if (x > 1) res = res / x * (x - 1);
return res;
}

遞推求尤拉函式值：

 //篩選法打尤拉函式表   
 #define Max 1000001  
 int euler[Max];  
 void Init(){   
      euler[1]=1;  
      for(int i=2;i<Max;i++)  
        euler[i]=i;  
      for(int i=2;i<Max;i++)  
        if(euler[i]==i)   //i是素數
            for(int j=i;j<Max;j+=i)  
               euler[j]=euler[j]/i*(i-1);  
 }  

尤拉定理：

n與a互質，那麼

當n為素數時，a^(n - 1) ≡ 1 （mod n）,這個稱為費馬小定理

引入一個定理：

定理：如果a1,a2,....aphi(m)是個與m互素的完全剩餘系，且(k,m) = 1，那麼(ka1,ka2,.....kaphi(m))仍然是一個與m互素的完全剩餘系。

證明：  首先kai與m顯然互素，然後用反證法假設kai ≡  kaj (mod m) 和（k，m） = 1，可以推出ai = aj 這即證明了定理

先證明費馬小定理

證：

構造素數p

的既約剩餘系

因為，由引理3可得

也是p的一個既約剩餘系。由既約剩餘系的性質，

即

易知，同餘式兩邊可約去

，得到

這樣就證明了費馬小定理。

模仿上面證明，我們來證明尤拉定理，同樣利用構造既約剩餘系來證明。

證明尤拉定理

證：

我們構造模m的一個既約剩餘系xi，那麼

(ax) ≡x （mod m）

即a^phi(m) x  ≡x （mod m）

由於x與m既約，直接消去，得a^phi(m) ≡ 1 （ mod m）

其中a和m既約，即(a,m) = 1，

.QED