快速冪和矩陣快速冪

 2023/04/05 

快速冪 & 矩陣快速冪

快速冪

今天介紹快速冪和矩陣快速冪，首先今天如果我們要計算$3^{15}$這個數字，我們可以用一個簡單的迴圈去實作。

typedef long long ll;
ll num = 1;
for(int i = 0;i < 15;i++)
    num *= 3;
cout << num << '\n';

然而這樣的計算方式的時間複雜度是$O(N)$，太慢了，因此我們可以使用快速冪的方式來計算的更快，至於快速冪是什麼就是我接下來要介紹的。他的做法會先將指數的部分轉成二進位，以$3^{15}$為例，$(15)_{10}$轉成二進位為$(1111)_{2}$，那我們就可以把$3^{15}$寫成$3^8\times 3^4 \times 3^2 \times 3^1$，只需要計算4次。再舉一個例子，$5^{18}$這個數，將$(18)_{10}$轉成二進位為$(10010)_{2}$這麼一來我們就可以把$5^{18}$寫成$5^2\times 5^{16}$，只需要計算5次。經過上述兩個例子應該可以看出當該指數轉成二進位時，如果該位數為1即乘上底數的$2^{位數 - 1}$，而這樣的計算就是我們的快速冪啦，這樣計算的時間複雜度為$O(lgN)$，很明顯的我們從原本迴圈一次一次計算的$O(N)$進步到了$O(lgN)$。接下來我們將這樣計算的code寫出來。遞迴版本

#include <bits/stdc++.h>
#define IOS ios_base::sync_with_stdio(false);cin.tie(0);cout.tie(0);
using namespace std;
typedef long long ll;
ll mod = 1000000007;
ll fast_pow(int base, int exp){
    if(exp == 0) return 1;
    ll res = 1;
    if(exp & 1){
        res = fast_pow(base, exp >> 1);
        return res * res * base % mod;
    }
    res = fast_pow(base, exp >> 1);
    return res * res % mod;
}
int main(){
    IOS
    int base = 3, exp = 15;
    cout << fast_pow(base, exp) << '\n';

    return 0;
}

迴圈版本

#include <bits/stdc++.h>
#define IOS ios_base::sync_with_stdio(false);cin.tie(0);cout.tie(0);
using namespace std;
typedef long long ll;
ll mod = 1000000007;
ll fast_pow(int base, int exp){
    ll res = 1;
    while(exp > 0){
        if(exp & 1) res = res * base % mod;
        base = base * base % mod;
        exp >>= 1;
    }
    return res;
}
int main(){
    IOS
    int base = 3, exp = 15;
    cout << fast_pow(base, exp) << '\n';

    return 0;
}

由於可能會有超大的整數，因此我們通常會$MOD$一個很大$prime$。 $$$$

矩陣快速冪

首先要懂得矩陣乘法，矩陣乘法如下

\[ A= \left[ \begin{matrix} a & b \\ c & d \end{matrix} \right] \tag{1}, B= \left[\begin{matrix} e & f \\ g & h \end{matrix}\right] \]

\[ A\times B = \left[ \begin{matrix} ae+bg & af+bh \\ ce+dg & cf+dh \end{matrix} \right] \] 說白了就是前列後行做內積而已。再來我們要寫一個矩陣乘法的code，矩陣乘法的時間複雜度為$O(N^3)$，這邊我們拿二階方陣當作示範。

struct Matrix{
    ll mat[2][2] = {{0}};
    Matrix operator * (Matrix &inp){
        Matrix tmp;
        for(int i = 0;i < 2;i++){
            for(int j = 0;j < 2;j++){
                for(int k = 0;k < 2;k++){
                    tmp.mat[i][j] = ((tmp.mat[i][j] + (mat[i][k] % mod) * (inp.mat[k][j] % mod)) % mod) % mod;
                }
            }
        }
        return tmp;
    }
};

我把它寫成一個struct，並用operator overloading來重載*。接著我們要來寫快速冪的部分了。 code:

Matrix fast_pow(int n){
    if(n == 1) return base;
    if(n % 2 == 0){
        Matrix res = fast_pow(n >> 1);
        return res * res;
    }
    Matrix res = fast_pow(n >> 1);
    return base * res * res;
}

這其實就是剛剛寫的遞迴版的快速冪，寫法有一點點不一樣而已，迴圈版就不再演示，因為也跟剛剛是差不多的。再來我們將整個矩陣快速冪組裝起來就變成下面的code了。

#include <bits/stdc++.h>
#define IOS ios_base::sync_with_stdio(false);cin.tie(0);cout.tie(0);
#define INF 0x3f3f3f3f
using namespace std;
typedef long long ll;
ll mod = 1000000007;
struct Matrix{
    ll mat[2][2] = {{0}};
    Matrix operator * (Matrix &inp){
        Matrix tmp;
        for(int i = 0;i < 2;i++){
            for(int j = 0;j < 2;j++){
                for(int k = 0;k < 2;k++){
                    tmp.mat[i][j] = ((tmp.mat[i][j] + (mat[i][k] % mod) * (inp.mat[k][j] % mod)) % mod) % mod;
                }
            }
        }
        return tmp;
    }
};
Matrix base;
Matrix fast_pow(int exp){
    if(exp == 1) return base;
    if(exp % 2 == 0){
        Matrix res = fast_pow(exp >> 1);
        return res * res;
    }
    Matrix res = fast_pow(exp >> 1);
    return base * res * res;
}
int main(){
    IOS
    base.mat[0][0] = 1;
    base.mat[0][1] = 4;
    base.mat[1][0] = 2;
    base.mat[1][1] = 3;
    Matrix output = fast_pow(10);
    for(int i = 0;i < 2;i++){
        for(int j = 0;j < 2;j++){
            cout << output.mat[i][j] << ' ';
        }
        cout << '\n';
    }

    return 0;
}

到這邊我想應該可以清楚看出矩陣快速冪就是矩陣乘法加上快速冪而已，真的很有趣。

$$$$

快速計算費氏數列

通常我們會使用遞迴或是迴圈的方式計算費氏數列，因為我們有一個$f(i) = f(i - 1) + f(i - 2)$的關係式，我們就可以推得費氏數列的第$N$項為多少。然而我們可以使用矩陣乘法的方式得到費氏數列的第$N$項，並利用快速冪來增加運算速度。首先我們先看這個矩陣 \[ A = \left[ \begin{matrix} 1 & 1 \\ 1 & 0 \end{matrix} \right] \] 這一個矩陣是一個2階方陣，看似平平無奇的矩陣卻是我們今天的主角。先講結論，如果我們想要求費氏數列第$n$項，我們算出$A$矩陣的的$n - 1$次方，而此費氏數列的第$n$項就會是這個矩陣的第一個元素，也就是 \[ fib(n) = A^{n - 1}[0][0] \] \[ {\left[ \begin{matrix} 1 & 1\\ 1 & 0 \end{matrix} \right]}^{n-1}= \left[ \begin{matrix} fib(n) & ...\\ ... & ... \end{matrix} \right] \] 也就是左上方那個元素就是我們費氏數列第$n$項

$$$$

Proof by Mathematical induction

到現在我們有費氏數列的兩種計算方式，一個氏藉由遞迴關係式遞推而來，一個是藉由矩陣乘法來推得費氏數列第$N$項。 \[ A= \left[ \begin{matrix} 1 & 1\\ 1 & 0 \end{matrix} \right] \tag{1} \] \[ f(n)=\begin{cases} 1 &\mbox{if } n\mbox{=1 or 2} \\ fib(n - 1) + fib(n - 2) &\mbox{n } \geq \mbox{3} \end{cases} \tag{2} \]

#include <bits/stdc++.h>
#define IOS ios_base::sync_with_stdio(false);cin.tie(0);cout.tie(0);
#define INF 0x3f3f3f3f
using namespace std;
typedef long long ll;
ll mod = 1000000007;
struct Matrix{
    ll mat[2][2] = {{0}};
    Matrix operator * (Matrix &inp){
        Matrix tmp;
        for(int i = 0;i < 2;i++){
            for(int j = 0;j < 2;j++){
                for(int k = 0;k < 2;k++){
                    tmp.mat[i][j] = ((tmp.mat[i][j] + (mat[i][k] % mod) * (inp.mat[k][j] % mod)) % mod) % mod;
                }
            }
        }
        return tmp;
    }
};
Matrix base;
Matrix fast_pow(int n){
    if(n == 1) return base;
    if(n % 2 == 0){
        Matrix res = fast_pow(n >> 1);
        return res * res;
    }
    Matrix res = fast_pow(n >> 1);
    return base * res * res;
}
int main(){
    IOS
    base.mat[0][0] = 1;
    base.mat[0][1] = 1;
    base.mat[1][0] = 1;
    base.mat[1][1] = 0;
    Matrix output = fast_pow(20);
    cout << output.mat[0][0] << '\n';

    return 0;
}

$$$$

Python version

import time
mod = 1000000007
def fib(n: int) -> int:
    dp = {}
    dp[0], dp[1] = 0, 1
    for i in range(2, n + 1):
        dp[i] = (dp[i - 1] % mod + dp[i - 2] % mod) % mod
    return dp[n]

if __name__ == '__main__':
    output = []
    for i in range(10000):
        fib(i)

import time
mod = 1000000007
class Matrix:
    def __init__(self):
        self.mat = [[0 for i in range(2)] for j in range(2)]
    def mat_mul(self, inp: list):
        tmp = Matrix()
        for i in range(2):
            for j in range(2):
                for k in range(2):
                    tmp.mat[i][j] = ((tmp.mat[i][j] + (self.mat[i][k] % mod) * (inp[k][j] % mod)) % mod) % mod
        return tmp

def fast_pow(exp: int, base: Matrix) -> Matrix:
    res = Matrix()
    res.mat[0][0] = 1
    res.mat[0][1] = 0
    res.mat[1][0] = 0
    res.mat[1][1] = 1
    while exp > 0:
        if exp & 1:
            res = res.mat_mul(base.mat)
        base = base.mat_mul(base.mat)
        exp >>= 1
    return res


if __name__ == '__main__':
    output = []
    for i in range(10000):
        n = i
        base = Matrix()
        base.mat[0][0] = 1
        base.mat[0][1] = 1
        base.mat[1][0] = 1
        base.mat[1][1] = 0
        fast_pow(n - 1, base)

Reference

Next Post

GIT
Previous Post

在Linux安裝CUDA，並在Pytorch進行CUDA訓練

CATALOG

1. 快速冪 & 矩陣快速冪