數論專題（二）數論基礎知識

轉載自：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;
    }

      這個函數揭示了一個約定俗成的概念，即任何非零整數和零的最大公約數為它本身。

     【例題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;
        }
    }

      擴展歐幾里德算法和歐幾里德算法的返回值一致，都是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個等式：

圖二-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所示：

圖二-4-1

圖二-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;
                } 
            }
        }
    }

       這里的phi[i]保存了i這個數的歐拉函數，還是利用素數篩選將所有素數篩選出來，然后針對每個素因子計算它的倍數含有該素因子的個數，利用歐拉公式計算該素因子帶來的歐拉函數分量，整個篩選過程可以參考素數篩選。

      c、歐拉定理和費馬小定理

      歐拉定理：若n,a為正整數，且n,a互素，則: 。

      費馬小定理：若p為素數，a為正整數且和p互素，則： 。

      由於當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的三個圓的並面積，則是先將三個圓的面積相加，然后減去兩兩相交的部分，而三個圓相交的部分被多減了一次，予以加回。

圖二-5-1

圖二-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 );
        }
    }

      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個數的組合。