Introduction

數論是什麼？

• 數論是屬於純數學的一個領域，專注於研究「整數」的性質。

• 因為很常為基礎數學訓練的一個項目，所以也被稱為「The Queen of Mathematics」

• 在資訊領域中最常應用到的是在密碼學

• 競賽中數論算是個通常用到的coding技巧不難，但是吃重數學的思考跟小技巧的領域。

• 意思是如果你會coding量少到讓你開心？

我們等等會講啥？

• 基礎知識
• 因數與質數
• 著名的數論定理
• 循環節
• 中國剩餘定理
• 質數專題

基礎知識

• 數論通常會用到啥呢？

+ , - , * , /

• 有人不知道這些操作在幹嘛嗎 OA O?

MOD, 模, 同餘

"%"

從除法開始

• 給你兩個數字\(a, m\)

• 計算\(a\)「除以」\(m\)的結果

  \(a = b \times m + r\)

• 這裡特別有\(0 \le r \le m-1\)

同餘

• 我們會說\(a,b\)在\(m\)下同餘意思是，

  \(a,b\)除以\(m\)的餘數相同

• 數學上我們會用

  \(a \equiv r \bmod{m}\) 表示

同餘系

• 在\(\mod{m}\)之下

  通常我們會用餘數\(r\)來分類整數們

  因為\(0 \le r \le m-1\)

  所以我們可以把整數分成\(m\)類

• 我們可以寫\(\mathbb{Z}_m\)

  來表示\(\mod{m}\)下的同餘系

在程式裡呢？

• 在C++裡面可以用「%」來算出餘數

  a % m = r

• 但要特別特別注意的是，在C++裡面當\(a\)是負數時

  a % m = r所獲得的r

  \(-m+1 \le r \le 0\)

GCD

• Greatest Common Divisor

• 通常數學上會寫為\(\gcd\)

LCM

• Least common multiple

• 對於兩個數字\(a, b\)，有\(ab = \gcd(a,b) \times lcm(a,b)\)

要怎樣算呢？

• 要計算GCD我們會使用輾轉相除法

• 用到性質如果\(a=br +c\)

• 會有\(\gcd(a,b)=\gcd(b,c)\)

  有趣的是\(a \ge b, b > c\)

• 這樣遞迴下去就會得到最大公因數

GCD code

int gcd(int a,int b){
    if(b==0) return a;
    else return gcd(b,a%b);
}

GCD 精簡版

std::__gcd(a,b)

• 耶XD

快速冪

• 我們很常要算\(a^n \mod{m}\)

• 當然可以直接\(O(n)\)一個一個乘

  但是\(n\)很大就糟了

• 因此我們需要快速冪

基本想法

巧妙地利用\(1,a,a^2,a^4,a^8 \ldots\) 等\(\bmod{m}\)的結果，來算出\({a^n}\bmod{m}\)的答案

答案初始值設為1，將\(n\)按照二進位分拆

如果\(n\)在二進位的第\(i\)的數值是\(1\)，那我們就幫答案乘上\(a^{2^{i}}\)

基本想法

• 不難看出是一個\(O(\log{n})\)的算法。

快速冪 Code

int pow(int a, int b, int mod) {
    int ans = 1, tmp = a;
    while(b) {
        if(b&1) ans = ans * tmp % mod;
        tmp = tmp * tmp % mod;
        b >>= 1;
    }
    return a;
}

反元素, 模逆元

我們常會想要在模底下求某個數字的反元素

也就是找到\(a^{-1}\)使得\(aa^{-1} \equiv 1 \mod{m}\)，

如果在\(\mod{m}\)底下\(a\)要有反元素，就必須有\(\gcd(a,m)=1\)，

裴蜀定理

• 對於任意整數\(a, m\)，存在整數\(x,y\)使得

  \(ax+my=\gcd(a,m)\)

• 如果今天\(\gcd(a,m)=1\)，那就會有，\(ax+my = 1\)

• 也就是\(ax \equiv 1 \bmod{m}\)

• 此時的\(x\)其實就是\(a\)的反元素，

• 因此可以改進GCD的計算，讓我們計算出\(a\)的反元素。

Extended Euclidean Algorithm

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

因數 質數

因數

• 一個整數的因數們是整數重要的性質

• 如果想要計算一個數字的所有因數

• 可以利用一個重要的性質

• 如果\(a=xy\)，則會有 \(x\le \sqrt{a}\) 或者 \(y \le \sqrt{a}\)

  所以只需要檢查所有小於等於\(\sqrt{a}\)的數字

  就可以算出\(a\)所有的因數，複雜度為\(O(\sqrt{a})\)

因數 Code

vector<int> factor(int a){
    vector<int> res;
    for(int i=1;i<sqrt(a);i++){
        if(a%i==0){
            res.push_back(i);
            res.push_back(a/i);
        }
    }
    int s=(int)sqrt(a);
    if(a%s==0)res.push_back(s);
}

重要!!!

• 「如果\(a=xy\)，則\(x\le \sqrt{a}\) or \(y \le \sqrt{a}\)」

• 常常可以利用這個性質做出複雜度根號的做法

質數

• 質數在數論中可是很重要的東西

• 從因數的觀點來看是一種整數的單元

• 而在模質數\(p\)底下除了0以外的所有數字都存在模逆元！

• 就是我們昨天講到\(\mathbb{Z}_m\)是個Filed這件事(還記得嗎？)

• 這代表著一件有趣的事，我們可以在模\(p\)底下進行「除法」。

質數判定

• 想要把整數用質數來處理，當然要先知道要怎樣判斷一個質數

• 如果我們想要判定單一個數字是不是質數，可以簡單地利用因數分解，只要一個數的因數只有1和自己就是質數，

Primality Test

bool isPrime(int a) {
    if(a==1) return false;
    else return (factor(a).size() == 2);
}

• 利用剛剛寫過的因數分解，需要\(O(\sqrt{a})\)的時間判斷。

• 最後我們會介紹更快驗質數的方法

篩法

• 如果今天想要重複的問不同數字是不是質數

• 我們可以開個陣列存某個數字是不是質數，利用篩法來算出有哪些質數

篩法中心精神

• 每遇到一個質數\(p\)，我們就把所有\(p\)的倍數\(2p, 3p \ldots\)在陣列中都標記成\(false\)，將它們「篩掉」

• 一樣每個正整數只需要知道自己「根號」以內的數字都不整除自己，該正整數就是個質數

• 因此只需要篩到\(\sqrt{n}\)就可以確定\(n\)以內有哪些是質數

質數個數

給定整數\(n\)，請問有多少質數小於\(n\)？

質數個數 Code

int prime[N], ptop, numOfPrime[N];
bool isPrime[N];
void init(){
    pTop=0;
    for(int i=2;i<N;i++)isPrime[i]=true;
    for(int i=2;i<N;i++){
        if(isPrime[i]){
            prime[pTop++]=i;
            for(int j=2*i;j<N;j+=i)isPrime[j]=false;
            numOfPrime[i]=numOfPrime[i-1]+1;
        }
        else numOfPrime[i]=numOfPrime[i-1];
    }
}

質數個數\(\pi(n)\)

而數學上對於質數個數有這樣個重要的結論

\(\pi(n) \approx \frac n {\ln n}\)

質因數分解

一個整數N我們可以分解成他質因數冪次的乘積

\(N = {p_1}^{e_1} \times {p_2}^{e_2} \times ... \times {p_n}^{e_n}\)

通常我們會用vector of pairs來表示\(N\)質因數分解後的結果

質因數個數約莫是\(O(\log{N})\)

實作?

• 我們可以掃過所有小於\(\sqrt{N}\)的所有質數，然後用while從把N除掉

• 但是這樣最差還是\(O(\sqrt{N})\)

• 有沒有辦法只看質因數個數次？

實作

• 我們可以利用篩法，對於每個數字記錄他其中一個質因數

• 就是primeFactor[i]

修改的質數篩法片段

        if(isPrime[i]){
            prime[pTop++]=i;
            firstPrime[i]=i;
            for(int j=2*i;j<N;j+=i)
                isPrime[j]=false,primeFactor[j]=i;
            numOfPrime[i]=numOfPrime[i-1]+1;
        }

質因數分解 Code

讓我們回到質因數分解

typedef pair<int,int> PII;
typedef vector<PII> VPII;
#define pb(x) push_back(x)
#define mp(x,y) make_pair(x,y)
VPII factorInt(int n){
    VPII res;
    while(n>1){
        int p = primeFactor[n];
        int e =0;
        while(n%p==0)n/=p,e++;
        res.pb(mp(p,e));
    }
    return res;
}

應用

• 質因數分解可以做啥呢？

• 因數個數\(d(N)\)

  \(d(N) = \prod\limits_{i=1}^n {(e_{i}+1)}\)

• 尤拉函數 \(\varphi(N)\)

  \(\varphi(N)= N \times \prod\limits_{i=1}^n {(1- \frac{1}{p_i})}\)

排容原理

我們可以利用質因數來做排容原理。

給定整數\(N,K\)， 請輸出一個答案表示在\([1,N]\)之中與\(K\)互質的整數個數。

Solution

• 想要計算一個數字\(k\)在\([1,N]\)中有多少倍數

  可以簡單地用\([N/k]\)做到

• 從質因數\(S = \{p_1,p_2,...,p_n\}\)中，窮舉所有的非空子集合

• 對於任意一種集合\(s\)，計算\([1,N]\)中有多少是這個\(s\)裡面質數乘積的倍數

• 然後乘上\((-1)^{|s|}\)加起來

Solution

• 所以答案其實是

  \(N - (\sum [N/p_i] - \sum [N/(p_ip_j)] + \sum [N/(p_ip_jp_k)]......)\)

Code

typedef vector< pair<int, int> > VPII;
int solve(int N, int K) {
    VPII factor = factorInt(N);
    int ans = 0;
    for(int i = 1; i<((1<<factor.size())-1) ; i++){
        int tmp = 1, now = i;
        int count = 0, idx = 0;
        while(now){
            if(now&1){
                count++;
                tmp *= factor[idx].first;
            }
            idx++, now>>=1;
        }
        int s = (count%2 == 0) ? 1:-1;
        ans += s*(N/tmp);
    }
    return N - ans;
}

綜合應用

給定\(N\)和\(K\)個小於等於\(N\)的質數 \(p_1 \ldots p_K\) ，\(N\)代表有標號\(1\)到\(N\)的\(N\)盞燈，每個質數\(p_i\)是一個開關，按下去的瞬間會把 \([1,N]\) 之中所有是這個質數倍數的燈會狀態改變，亮的變暗，暗的變亮。如果一開始所有燈都是暗的，你可以自由地使用質數開關，請問最後最多可以有多少燈是亮的呢？

有名的數論定理

費馬小定理

\(a^p \equiv a \bmod{p}\)

特別的是如果\(\gcd(a,p)=1\)，可以有

\(a^{p-1} \equiv 1 \bmod{p}\)

反元素 again

• 如果我們想要計算一個整數\(a\)在\(\mod p\)之下的反元素\(a^{-1}\)

  在\(\gcd(a,p) = 1\)的狀況下我們知道，\(a^{p-1} \equiv 1 \bmod{p}\)

• 有趣的是\(a\times (a^{p-2}) \equiv 1 \bmod{p}\)

  所以\((a^{p-2})\equiv a^{-1} \bmod{p}\)

反元素 Code

int inverse(int a, int p){
    return pow(a, p-2, p); 
}

尤拉函數

尤拉函數是在計算對於一個整數\(N\)， 小於等於\(N\)的整數，且跟\(N\)互質的有多少個，常用\(\varphi(N)\)表示。

\(\varphi(N)= N \times \prod\limits_{i=1}^n {(1- \frac{1}{p_i})}\) \(\prod\limits_{i=1}^n {{p_i}^{e_i-1}*(p_i-1)}\)

令人困擾的\(\gcd\)

給定整數\(N, d\)，請問\([1,N]\)中有多少整數\(i\)， \(\gcd(i,N) = d\)呢？

Solution

因為\(\gcd(i,N)\)整除\(N\)，所以如果\(d\)不是\(N\)的因數就答案一定是0。

如果\(N\)是\(d\)的倍數，那 \(N/d\)就會是個整數。

• 又因為\(\gcd(i,N)=d\)，所以\(i\)一定是\(d\)的倍數，在\([1,N]\)中其實只有\(d, 2d ...,j\times d, ..., (N/d) \times d\)有可能是\(i\)，讓我們重新回到等式並做些化簡

  \(d = \gcd(j\times d, N) = d \times \gcd(j,N/d)\)

Solution

• 從上面等式中我們會發現我們其實要找

  \(j \in [1,N/d]\)且\(\gcd(j,N/d) =1\)

• 這樣的個數就是

  「小於等於\(N/d\)且與其互質的正整數個數」

  \(\varphi(N/d)\)

尤拉定理

尤拉定理是比較general版本的費馬小定理，如果\(\gcd(a,n)=1\)

\(a^{\varphi (n)} \equiv 1 \bmod{n}\)

特別當\(n\)是質數時是費馬小定理。

又是反元素

如費馬小定理一般，若\(a\)與\(n\)互質， 我們可以利用尤拉定理求出\(a\)在\(n\)下的模逆元， 也就是\(a^{\varphi(n)-1}\)。

來點應用吧？

• RSA加密演算法

加密算法

• 在討論RSA之前，我們應該先來談談什麼是加密

• 通常就是兩個人想要傳訊息，但是有不想讓中間的人看到的訊息

• 最常見的加密方法都會有公鑰與私鑰，我們會將公鑰傳給「會傳東西給我們」的人，私鑰自己保管好

• 如果要傳訊息，那就是擁有公要的人用公鑰把東西加密，我們收到之後用私要解密

RSA製作key

• \((N,e)\)是公鑰，\((N,d)\)是私鑰。

加密

• 如果我們要傳送\(m\)，我們可以先把它分拆成每個單元都小於\(N\)的數字\(n\)

• 對於一個\(n\)我們把東西加密成\(c\)

  \(n^e \equiv c\ (\mathrm{mod}\ N)\)

• 然後把東西送出去

解密

• 收到東西之後我只要計算

  \(c^d \equiv n\ (\mathrm{mod}\ N)\)

  就可以還原出\(n\)

解密原理

• \(c^d \equiv n^{e \cdot d}\ (\mathrm{mod}\ N)\)

  以\(及ed = 1 (\mathrm{mod} p-1)和ed = 1 (\mathrm{mod} q-1)\)

• 由歐拉定理可證明（因為\(p\)和\(q\)是質數）

  \(n^{e \cdot d} \equiv n\ (\mathrm{mod}\ p)\)和\(n^{e \cdot d} \equiv n\ (\mathrm{mod}\ q)\)

• 這說明

  \(n^{e \cdot d} \equiv n\ (\mathrm{mod}\ pq)\)

Wilson's theorem

對於質數\(p\)來說，

\((p-1)! \equiv -1 \bmod{p}\)

C(n,m)

\(C(n,m) = \frac{n!}{m!(n-m)!}\)

預處理 f[i] ,invf[i]。

\(C(n,m) = C(n-1,m) + C(n-1,m-1)\)

循環節

整數在\(\mod{m}\)下是一個元素有限的集合，又寫作\(\mathbb{Z}_m\)。

所以當一個整數\(a\)，把他的所有次方\(1, a, a^2, a^3 \cdots\)列出來， 因為這數列可以無限的下去，鴿籠原理告訴你說一定存在兩個不同的\(i,j\)，其中\(i<j, a^i \equiv a^j \mod{m}\)，對於任意大於等於\(i\)的整數\(k\)

\(a^k \equiv a^{k-i}a^{i} \equiv a^{k-i}a^{j}\)

\(\equiv a^{k +(j-i)} \mod{m}\)

因此對任意這樣的\(k\)來說， 任何正整數\(m\)都會讓下面等式滿足，

\(a^{k} \equiv a^{k+ m\times(j-i)} \mod{m}\)

這時候其實我們就找到了所謂的循環節！

Recall: 費馬小定理 或是 尤拉定理

有趣的結論

• 這...不好說啊

再來是只要從\(1, a, a^2, a^3, a^4...\)找下去，

我們可以用map去記錄看過的數字

在\(O(\log{m} +\varphi(m))\)的時間內一定會找到最短的循環節，

簡單的解釋\(\log{m}\)的存在， 一開始不一定會出現循環節的原因， 其實是因為那些\(a\)與\(m\)之間共有的質因數。

如果那些質因數在\(a^i\)次方夠大時超過了\(m\)裡面同樣質因數的次方數， 那該質因數對於mod就不再會有影響， 所以約莫要花上\(\log{m}\)的時間來做到這件事。

這兩個結論告訴我們， 通常可以很放心的直接一個個跑下去找循環節， 很快就會找到。 很重要的是循環節通常會用兩個數字表示:

一個代表從哪個次方會開始循環。

一個代表循環節長度。

次方次方

給定正整數\(a_1, a_2 ,a_3\)， 再給定一個數字\(m\)，

請問\(a_1^{a_2^{a_3}} \bmod{m}\) 的答案應該要是多少?

Solution

• 如果\(a_1\)與\(m\)互質，我們可以很簡單的用尤拉定理，

• 因為\({a_1}^{\varphi(m)}\equiv 1 \mod{m}\)，

  所以可以讓\(t\equiv{a_2}^{a_3} \mod{\varphi(m)}\)，

• 然後就可以直接用\({a_i}^t\)算出答案。

• 但如果不互質的話，

  我們其實可以去跑\(a_1\)在\(m\)中的起點與循環節，

  假設起點是\(k\)循環節是\(r\)。

要先判斷大小，

如果\({a_2}^{a_3} < k\)我們直接去計算答案就好，

因為\(k\)不會太大。

如果\({a_2}^{a_3} \ge k\)

那我們其實可以去計算\(t \equiv {a_2}^{a_3}-k \mod{r}\)，

那答案就是\(a^k \times a^t \mod{m}\)。

好寫的小技巧

對於 \(a^k \equiv x \mod{m}\)

• 訂個上限 \(len\)
• 如果 \(k <= len\) 就直接算 \(a^k \mod{m}\)

• 如果 \(k > len\)

  \(k \equiv k' \mod{phi(m)}\)

  就算 \(a^{k' + phi(m)}\)

中國剩餘定理

中國剩餘定理

這個定理想必大家耳熟能詳?

孫子算經中曾有段話:

「今有物不知其數，三三數之賸二，五五數之賸三，七七數之賸二，問物幾何？」

答曰：「二十三」

或許這就是中國剩餘定理最早的來由吧。

韓信點兵

劉邦想要知道韓信到底有多少能力可以帶兵打仗， 他想要考考韓信估計兵力的能力， 於是問個問題說：

「如果每個帳篷住十個人會剩下兩個，

如果每餐都有十一個人坐在火爐邊吃飯那就會剩下三人，

如果每個隊伍有十三個人就會有一小隊只有五人，

那這個軍隊中應該要有多少人呢？」

雖然韓信很聰明， 但想說如果每天劉邦都亂出一題新的來問就麻煩了， 所以想請你寫個程式來幫他算算答案。

劉邦出題的方式很固定，每次都是\(n\)個問句， 每個問句都是說如果\(m_i\)個人一團會剩下\(a_i\)個人， 這裡特別有\(m_i\)裡面兩兩互質， 因為不互質劉邦自己也不太會。

想請你算出最小可能的軍隊人數\(x\)是？

\(\left\{ \begin{matrix} x \equiv a_1 \pmod {m_1} \\ x \equiv a_2 \pmod {m_2} \\ \vdots \qquad\qquad\qquad \\ x \equiv a_n \pmod {m_n} \end{matrix} \right.\)

Solution

因為\(m_i\)之間兩兩互質，所以一定有解，

我們直接構造答案來證明這件事，這也順便給出了一種算出答案的做法。

• 首先設\(M = \prod_{i=1}^n m_i\) ，並設\(M_i = M/m_i\)。

• 然後我們可以計算\(t_i\)為\(M_i\)在\(\mod{m_i}\)下的模逆元，

  也就是\(t_iM_i \equiv 1 \mod{m_i}\)。

• 因為\(M_i\)為除了\(m_i\)之外的\(m_j\)的乘積，所以與\(m_i\)互質，

  這裡可以使用尤拉定理或是Extended Euclidean Algorithm來算出模逆元\(t_i\)。

• 最後答案\(\sum\limits_{i=1}^{n} a_it_iM_i \mod{M}\)就會是的結果。

  為甚麼答案是對的呢？

• 因為對於任意的\(m_j\)，

  \(\sum\limits_{i=1}^{n} a_it_iM_i \equiv a_jt_jM_j \equiv a_j \mod{m_j}\)

  然後\(M\)被任意的\(m_j\)整除。

  這樣可以獲得一個\(O(n \log{m})\)的中國剩餘定理做法。

進位制

進位制

• 所謂的\(B\)進制，就是一個數每位的數碼是\(1\)到\(B-1\)，

  數碼逢\(B\)進位。

• 所以如果一個數字在\(B\)進位下是\(\overline{a_na_{n-1}...a_0}_{B}\)，

• 那這個數字在十進位之下就是\(a_n \times B^n + a_{n-1} \times B^{n-1}+......+a_0\)。

• 進位制是在不只是在數論中會探討，

  在程式語言中進位制也是個需要知道的知識，

• 我們平常在使用的「bit」就是把整數在電腦中做二進位轉換，

• 而我們也會常常看到所謂的HEX，就是數字在電腦中存成16進位的字串。

• 進位制可說是程式中的一種常識。

十進制 To \(B\)進制

平常我們最常用的進位制就是十進制，

在這章節中我們用\(a_{(10)}\)來表示\(a\)的十進制數碼，

如果今天我們要把\(a\)從10進制轉成轉成\(B\)進制，

10 to \(B\) Code

其實我們不久前才做過??

int BtoTen(int a_B[],int B,int top){
    int res = 0;
    int base = 1;
    for(int i=0;i<top;i++){
        res += base*a_B[i];
        base *= B;
    }
    return res;
}

B進制 To 十進制

• 給定一個數字在\(B\)進位下的數碼，

  要還原成原本的十進位數字，

• 就直接照著進位制的定義可以算出結果。

  \(a_n \times B^n + a_{n-1} \times B^{n-1}+......+a_0\)

\(B\) to 10 Code

int BtoTen(int a_B[],int B,int top){
    int res = 0;
    int base = 1;
    for(int i=0;i<top;i++){
        res += base*a_B[i];
        base *= B;
    }
    return res;
}

進位制的應用

• 來談談一種進位制的應用，進位制很多時候可以用來簡化的表示一個定理的結果。

Lucas's Theorem

我們要談的是一種跟組合數計算有關的定理，

當我們要計算\(\binom{n}{m} \bmod{p}\)時，

如果\(p\)的數量級沒有那麼大，

我們可以利用這個定理快速的算出這個組合數的結果。

Lucas's Theorem

\(\binom{n}{m} \bmod{p} \equiv \prod\limits_{i=0}^{k} \binom{m_i}{n_i} \bmod{p}\)

\(m = \sum\limits_{i=0}^{k} m_i \times p^{i}, n = \sum\limits_{i=0}^{k} n_i \times p^{i}\)

特別的是如果\(m>n\)，那我們讓\(\binom{n}{m} =0\)

Nim Game

• 這種遊戲中文常叫做「拈」。

  是個經典的遊戲問題，

  其結論也跟進位制有關係。

Nim Game

今天如果我們有\(n\)堆石頭，第\(i\)堆石頭有\(a_i\)個石頭， 兩個人輪流玩遊戲，不能拿的人就輸了， 請問先手還是後手會贏？

Nim Game的結論是今天我們把所以\(a_i\)轉成二進位，

對於每一位數統計1的個數，

如果每一位1個數都是偶數，那就是後手會贏，

其他是先手會贏。

• 輸的盤面是所有石頭堆都變為零，

  每位的1個數都是偶數，

• 舉個例如果石頭有3堆，分別是\(3,4,5\)，

  他們二進制分別是\(011_{(2)}, 100_{(2)}, 101_{(2)}\)，

  第二位1的個數不是偶數，因此先手會贏。

• 要判定這件事，

  其實就是看\(t = \bigoplus\limits_{i=1}^n a_i\)是否\(t=0\)，

• 有趣的如果先手會贏，他的第一步就是尋找\(i\)使得\(a_i \oplus t < a_i\)的那堆\(i\)，

• 然後把\(a_i\)拿成\(a_i \oplus t\)，

• 這樣新的\(t' = \bigoplus\limits_{i=1}^n {a'}_i = t \oplus t = 0\)，

• 因此我們讓後手拿到XOR為零的盤面，

  也可以保證我們一定可以做到這件事。

Joseph Problem

\(N\)個人坐成一個圓，瞬時針編號為\(1, 2,\cdots, N\)， 從1號開始開始順時針報數\(1,2,1,2......\)， 只要數到2的人就會從圈圈中被消除， 請問最後活下來的是幾號呢?

如果用函數\(f(N)\)來表示N個人玩這個遊戲活下來的號碼， 那我們會得到下列遞迴式

• \(2f(m)-1\) 如果 \(N = 2m\) 是偶數

• \(2f(m)+1\) 如果 \(N = 2m+1\) 是奇數

• 這樣的\(f(N)\)的結論可以用進位置表示，

• 今天如果\(N = {b_nb_{n-1}...b_0}_{(2)}\)，最高位\(b_n =1\)，

• 那有趣的是答案會是\(f(N) = {b_{n-1}b_{n-2}...b_0b_n}_{(2)}\)。

質數專題

線性算法 linear sieve

希望每個數字都只被篩法篩到一次(原本篩法會被所有質因數遍歷)

這演算法設計上希望所有的合數\(m\)都只被最小的質因數篩到

因此我們改變篩法的方向，

對於每個數字\(m\)我們用那些小於它的那些質數\(p_i\)乘上\(m\)， 去更新\(m \times p_i\)，讓其被篩掉，

在\(p_i\)尚未整除\(m\)之前，\(p_i\)一定比\(m\)中最小的質因數還要小。

線性算法 Code

int p[N], np=0;
int a[N];     
for (int m = 2; m < N; m++) {
  if (a[m] == 0) p[np++] = m;
  for (i = 0; i < np && m*p[i] < N; i++) {
    a[m*p[i]] = 1;
    if (m%p[i] == 0) break;
  }
}

• 只要一發現某個質數\(p_i\)整除\(m\)，表示該質數是\(m\)最小的質數，

• 為了讓之後的數字都只被自己最小的質數篩掉

  就該從迴圈中跳出繼續窮舉。

費馬質數判定

• 根據費馬小定理我們知道\(a^{p-1}\equiv 1 \mod{p}\)，

• 對於一個未知的\(n\)，

  我們在區間\([1,n-1]\)中選擇很多個\(a\)來進行檢驗，

  如果結果都是1那在費馬質數判定中我們會說他是質數。

錯誤?!

費馬質數判定的缺點是

存在著叫做偽質數 (Carmichael Numbers)數字\(n\)，所有\([1,n-1]\)中的\(a\)都會滿足費馬小定理的等式， 但是\(n\)卻不是質數。

Miller–Rabin 質數判定法

• 要判定\(N\)是否為質數。

• 首先將\(N-1\)分解為\(2^sd\)。

• 在每次測試先隨機選一個介於\([1, N-1]\)的整數\(a\)，對所有的\(r \in [0, s-1]\)， 若\(a^d \mod N \neq 1且a^{2^{r}d} \mod N \neq -1\)， 則\(N\)是合數。 不然，\(N\)有\(3/4\)的「機率」為質數。

Code

bool miller_rabin(LL n, LL a) {
    if (__gcd(a, n) == n) return PRIME;
    if (__gcd(a, n) != 1) return COMPOSITE;
    LL d = n-1, r = 0, res;
    while(d%2==0) { ++r; d/=2; }
    res = bigmod(a, d, n); //計算 a 的 d 次方除以 n 的餘數
    if (res == 1 || res == n-1) return PRIME;
    while(r--) {
        res = mul(res, res, n); //小心乘法溢位
        if (res == n-1) return PRIME;
    }
    return COMPOSITE;
}

Code

bool isprime(LL n) {
    LL as[7]={2, 325, 9375, 28178, 450775, 9780504, 1795265022};
    for(int i=0;i<7;i++)
        if(miller_rabin(n, as[i]) == COMPOSITE)
            return COMPOSITE;
    return PRIME;
}

注意到上面程式碼當中的as陣列。

LL as[7]={2, 325, 9375, 28178, 450775, 9780504, 1795265022};

其實在 \(2^{64}\) 以內，只要用這幾個 \(a\) 檢查過之後， 就可以完全確定範圍內的數字是否為質數了！ 順帶一提，檢查 int 範圍內的質數只需要以 2, 7, 61 這三個數字為基底檢查就可以了。

把這幾個數字藏進 codebook 裡是個不錯的選擇唷！

Miller–Rabin

因為是一定「機率」是質數，所以要做夠多次就幾乎確定該數字是否為質數。

但是沒有費馬質數判定法的問題，次數夠多一定可以把不是合數的數字篩掉。

有趣質數 3*1家族

31, 331, 3331, 33331, 333331, 3333331, 33333331都是質數。

但33333331就不是了。

梅森質數

形式為\(2^k -1\)的質數

前四個是\(k =2,3,5,7\)，目前最大的有22,338,618位數。

\(k = 74207281\)。

The End