矩陣快速幂 算法原理與模闆

一.前導知識(線性代數)

1.矩陣乘法

※運算公式

※代碼表示

const int N = 100;
int c[N][N];
void matrix_multi(int a[][N], int b[][N], int n){
    memset(c, 0, sizeof(c));
    for(register int i = 1; i <= n; i++)
        for(register int j = 1; j <= n; j++)
            for(register int k = 1; k <= n; k++) 
                c[i][j] += a[i][k] * b[k][j];
}      

2.快速幂

● 此部分内容引自OI-WIKI

※實作原理

計算 的 次方表示将 個 乘在一起：。然而當 太大的時侯，這種方法就不太适用了。不過我們知道：。二進制取幂的想法是，我們将取幂的任務按照指數的 二進制表示 來分割成更小的任務。

首先我們将

因為 有 個二進制位，是以當我們知道了 後，我們隻用計算 次乘法就可以計算出 。

于是我們隻需要知道一個快速的方法來計算上述 3 的

是以為了計算 ，我們隻需要将對應二進制位為 1 的整系數幂乘起來就行了：

将上述過程說得形式化一些，如果把 寫作二進制為 ，那麼有：

其中 。那麼就有

根據上式我們發現，原問題被我們轉化成了形式相同的子問題的乘積，并且我們可以在常數時間内從 項推出

這個算法的複雜度是 的，我們計算了 個 次幂的數，然後花費

※代碼表示

​

​遞歸式快速幂​

​

long long binpow(long long a, long long b) {
  if (b == 0) return 1;
  long long res = binpow(a, b / 2);
  if (b % 2) return res * res * a;
  else return res * res;
}      

​

​循環式快速幂​

​

long long binpow(long long a, long long b) {
  long long res = 1;
  while (b > 0) {
    if (b & 1) res = res * a;
    a = a * a;
    b >>= 1;
  }
  return res;
}      

​

​模意義下取幂​

​

long long binpow(long long a, long long b, long long m) {
  a %= m;
  long long res = 1;
  while (b > 0) {
    if (b & 1) res = res * a % m;
    a = a * a % m;
    b >>= 1;
  }
  return res;
}      

注意：根據費馬小定理，如果 是一個質數，我們可以計算

二、矩陣快速幂

1.執行個體導入

給定整數，計算斐波那契數列的第項。

首先，觀察斐波那契數列的特征：遞推公式：

我們不難發現：與和與之間的關系可以用矩陣的乘法關系進行描述：

不妨設，改寫遞推公式，可得：

于是我們可以用矩陣乘法在

實際上，對于斐波那契數列，我們可以根據上述原理繼續優化：即​

​快速倍增法​

​。這裡不深入讨論。

根據上述題目，我們不難看出，對于具有類似特征的問題，我們可以根據題目構造相應的矩陣達到優化計算的目的。

2.矩陣快速幂模闆

const int N = 10;
int tmp[N][N];

void multi(int a[][N], int b[][N], int n){
    memset(tmp, 0, sizeof tmp);
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            for (int k = 0; k < n; k++)
                tmp[i][j] += a[i][k] * b[k][j];
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            a[i][j] = tmp[i][j];
}

int res[N][N];

void Pow(int a[][N], int n){
    memset(res, 0, sizeof res); //n是幂，N是矩陣大小
    for (int i = 0; i < N; i++) res[i][i] = 1;
    while (n){
        if (n & 1) multi(res, a, N); //res=res*a;複制直接在multi裡面實作了；
        multi(a, a, N); //a=a*a
        n >>= 1;
    }
}      

3.應用