演算法學習:矩陣快速冪/矩陣加速

AdviseDY發表於2024-08-11

1.前言

 其實本質上來說,矩陣快速冪或是矩陣加速的題目比較的模板化一些,大體上都是屬於我們要先寫出來一個遞推式子(或者是我們需要遞推的式子),然後由於遞推的次數過大,1e18之類的,會導致複雜度的飈升,所以我們會用到矩陣來幫我們快速處理。
 另外,從題目的型別上大概是分為兩類,一類是線性遞推式,再利用矩陣加速處理;另外一類就是圖上的問題,一般點數較小,可能100不到,然後對邊有限制要求之類的,例如需要走T步,T的範圍是[1,1e18],此時我們也可以用矩陣加速處理。
  本篇的主要內容以矩陣加速為主,矩陣快速冪的內容會比較草率。

2.矩陣快速冪

我們來大概瞭解一下矩陣快速冪,首先要確保我們是會快速冪的。以下貼出來快速冪的程式碼:

int kuai(int a,int b){
	int tem=1;
	while(b){if(b%2) tem=tem*a; a=a*a,b/=2;}
	return tem;
}

具體瞭解快速冪請百度,而矩陣的運算和實數的運算在結合律上是一致的,所以我們完全可以套用快速冪的形式寫出來矩陣快速冪:Mat是我設的矩陣的結構體

static Mat kuai(Mat A, int b) {
        Mat tem(A.n, A.m, 1, A.st);//引數分別代表,行,列,初始化內容,下標從1或者0開始
        while (b) { if (b % 2) tem = tem * A; b = b / 2, A = A * A; }
        return tem;
    }

另外此處貼上我的矩陣的板子,後面的程式碼展示中就省略了板子:

struct Mat {
    int a[130][130];
    int n, m, st;
    Mat(int n_, int m_, int fl = 0, int st_ = 1) {
        n = n_, m = m_, st = st_;
        if (!fl) rep(i, st, n) rep(j, st, m) a[i][j] = 0;
        if (fl == INF) rep(i, st, n) rep(j, st, m) a[i][j] = INF;
        if (fl == 1) rep(i, st, n) rep(j, st, m) a[i][j] = (i == j);
    }
    int* operator [] (int x) { return a[x]; }
    Mat operator * (Mat b) {
        Mat res(n, b.m, 0, b.st);//記得賦值
        rep(i, st, n) rep(j, st, b.m) rep(k, st, m) {
            res[i][j] += a[i][k] * b[k][j] % mod;
            res[i][j] %= mod;
        }
        return res;
    }
    Mat operator + (Mat b) {
        Mat res(n, b.m, 0, b.st);//記得賦值
        rep(i, st, n) rep(j, st, b.m) {
            res[i][j] = a[i][j] + b[i][j];
            res[i][j] %= mod;
        }
        return res;
    }
    static Mat kuai(Mat A, int b) {
        Mat tem(A.n, A.m, 1, A.st);
        while (b) { if (b % 2) tem = tem * A; b = b / 2, A = A * A; }
        return tem;
    }
    //指數從0開始的求和
    static Mat kuai_sum(Mat A, int b) {
        Mat tem(A.n, A.m, 1, A.st); Mat sum = tem;
        while (b) {
            if (b % 2) tem = sum + tem * A;
            sum = sum + sum * A; A = A * A, b /= 2;
        }
        return tem;
    }
    void out() {
        rep(i, st, n) { rep(j, st, m) cout << a[i][j] << " "; cout << endl; }
    }
};

3.矩陣加速

第一類:線性遞推式加速

學習一個東西是什麼,以及怎麼用,最好的辦法就是去看一道經典的例題。
下面來看一道矩陣加速的板題:洛谷P1962

題目描述:

\[F_n = \left\{\begin{aligned} 1 \space (n \le 2) \\ F_{n-1}+F_{n-2} \space (n\ge 3) \end{aligned}\right. \]

請你求出 \(F_n \bmod 10^9 + 7\) 的值。
對於 \(100\%\) 的資料,\(1\le n < 2^{63}\)

題解:

這裡我們已經有了遞推式,\(F_{n-1} + F_{n-2} = F_{n}\),正常我們只需要這樣一直遞推下去就好了,但是n的資料範圍是2的63次方,所以我們不得不考慮別的辦法:我們先把這兩個式子寫成矩陣的形式

\(I_1=\begin{bmatrix}F_{n-1} & F_{n-2} \end{bmatrix}\)

此時我們假設一下,這個矩陣需要乘上一個什麼矩陣可以變成

\(I_2=\begin{bmatrix}F_{n} & F_{n-1} \end{bmatrix}\)

在草稿紙上寫寫畫畫,我們就發現

\(A=\begin{bmatrix}1 & 1 \\ 1 & 0 \end{bmatrix}\)

\(I_1*A\)就會變成\(I_2\),而且,這個矩陣都是常數,且任意的一個\(I_{i}\)矩陣乘上它都會滿足變成\(I_{i+1}\)
所以我們從\(I_{i}\)一直乘A,可以變成我們想要的\(I_{n}\),此時在利用之前學習的矩陣快速冪,我們就可以在\(log\)的時間裡求出來我們最後的\(F_{n}\),另外再提一句,這個A矩陣的構建,是我們要先寫出來\(I_1\)(我們本身有的遞推式),然後將這個遞推式子的下標都+1,這是我們的目標矩陣,根據這兩個來倒推出我們的轉移矩陣,也就是A矩陣。我們會發現當[1,1]和[2,1]的位置都填1,就是我們斐波那契數列本身的遞推式子,然後\(F_{n-1}\)我們直接在[1,2]的位置寫1,再次利用我們上一個矩陣本身就有的即可。

此處貼上程式碼:

int n, m, T;
signed main() {
    read(), kuaidu();
    int st, ed;
    cin >> n;
    if (n == 1 || n == 2) cout << 1 << endl;
    else {
        Mat A(1, 2);
        A[1][1] = 1, A[1][2] = 1;
        Mat base(2, 2);
        base[1][1] = base[1][2] = base[2][1] = 1;
        A = A * Mat::kuai(base, n - 2);
        cout << A[1][1] << endl;
    }
    return 0;
}

第二類:在圖上的應用

我們還是接著來看一道例題:洛谷P2233

題目描述

在長沙城新建的環城公路上一共有 \(8\) 個公交站,分別為 A、B、C、D、E、F、G、H。公共汽車只能夠在相鄰的兩個公交站之間執行,因此你從某一個公交站到另外一個公交站往往要換幾次車,例如從公交站 A 到公交站 D,你就至少需要換 \(3\) 次車。

圖片來自洛谷

Tiger 的方向感極其糟糕,我們知道從公交站 A 到公交 E 只需要換 \(4\) 次車就可以到達,可是 tiger 卻總共換了 \(n\) 次車,注意 tiger 一旦到達公交站 E,他不會愚蠢到再去換車。現在希望你計算一下 tiger 有多少種可能的乘車方案。

輸入格式
僅有一個正整數 \(n\),表示 tiger 從公交車站 A 到公交車站 E 共換了 \(n\) 次車。
\(4\le n\le10^7\)

題解:

首先這個題我們很顯然的能夠知道,轉化成圖上就是有8個點,8個點只有互相相鄰的可以抵達,我們建立一個鄰接矩陣,其中給相鄰的能到達的賦值為1,(除了E這個點,因為到達了就不會再去別的車站了),E車站向其他車站的賦值是0。因為這裡n的數值很大,所以我們還是可以進行矩陣加速。
這裡稍微淺扯一下關於這裡為什麼圖上的鄰接矩陣我們依舊可以進行矩陣加速,先看矩陣乘法的程式碼實現:

Mat operator * (Mat b) {
        Mat res(n, b.m, 0, b.st);//記得賦值
        rep(i, st, n) rep(j, st, b.m) rep(k, st, m) {
            res[i][j] += a[i][k] * b[k][j] % mod;
            res[i][j] %= mod;
        }
        return res;
    }

這裡我們中間的式子,可以抽象的理解成是\(i\)\(k\)的路徑方案數乘以\(k\)\(j\)的路徑方案數,加到\(res[i][j]\)裡面,此刻\(res[i][j]\)代表的就是我們將兩個矩陣相乘後(意義上是步數相加)的方案數,原本這個鄰接矩陣代表著我們走1步的方案數,所以再進行矩陣快速冪即可。

signed main() {
    read(), kuaidu();
    int t;
    cin >> n;
    Mat A(7, 7);
    rep(i, 0, 7) {
        int x = i, y = (i + 1) % 8, z = (i - 1 + 8) % 8;
        A[x][z] = 1;
        A[x][y] = 1;
    }
    rep(i, 0, 7) A[4][i] = 0;
    Mat dp = Mat::kuai(A, n);
    // dp.out();
    // A.out();
    cout << dp[0][4] << endl;
    return 0;
}

到此為止,我們算是將矩陣加速的基礎題過了一遍,本人剩下的訓練內容來自洛谷的矩陣題單,可以自行搜尋。
我在這裡面的題當中抽出來幾道講講:

4.習題

4.1 P4159 [SCOI2009] 迷路

題目描述

該有向圖有 \(n\) 個節點,節點從 \(1\)\(n\) 編號,windy 從節點 \(1\) 出發,他必須恰好在 \(t\) 時刻到達節點 \(n\)
現在給出該有向圖,你能告訴 windy 總共有多少種不同的路徑嗎?
答案對 \(2009\) 取模。
注意:windy 不能在某個節點逗留,且透過某有向邊的時間為該邊長度。

輸入格式

第一行包含兩個整數,分別代表 \(n\)\(t\)
\(2\) 到第 \((n + 1)\) 行,每行一個長度為 \(n\) 的字串,第 \((i + 1)\) 行的第 \(j\) 個字元 \(c_{i, j}\) 是一個數字字元,若為 \(0\),則代表節點 \(i\) 到節點 \(j\) 無邊,否則代表節點 \(i\) 到節點 \(j\) 的邊的長度為 \(c_{i, j}\)

- 對於 \(100\%\) 的資料,保證 \(2 \leq n \leq 10\)\(1 \leq t \leq 10^9\)

題解:

這裡我們知道,如果有向邊的長度是1的話,就直接變成模板題了,但是長度是權值,我們不能這麼做。細心一點,會發現因為輸入的是字串,所以長度最大就是9,邊權的範圍是在1-9之間。我們可以像網路流那樣來進行拆點,每一個點\(i\)拆成9個點,\(i_1,i_2,...,i_9\),而我們每一次從一個點走,都是從\(i_1\)走,舉個例子,我們從\(i\)\(j\),權值是3,我們可以建立\(i_1 -> j_3\),權值是1的邊,然後再從\(j_3->j_2->j_1\),每一條邊權值都是1,這樣用拆點的方法解決了權值的問題,剩下的就和正常的板子題一樣了。

int fan(int x, int y) {
    return (x - 1) * 10 + y;
}
signed main() {
    read(), kuaidu();
    cin >> n >> m;
    Mat A(100, 100);
    rep(i, 1, n) {
        string S; cin >> S;
        int len = S.size();
        rep(j, 0, len - 1) {
            if (S[j] == '0') continue;
            int x = fan(i, 1), y = fan(j + 1, (S[j] - '0'));
            A[x][y] = 1;
        }
    }
    rep(i, 1, n) {
        rep(j, 2, 10) A[fan(i, j)][fan(i, j - 1)] = 1;
    }
    A = Mat::kuai(A, m);
    cout << A[fan(1, 1)][fan(n, 1)] << endl;
    return 0;
}

4.2 P3216 [HNOI2011] 數學作業

題目描述

給定正整數 \(n,m\),要求計算 \(\text{Concatenate}(n) \bmod \ m\) 的值,其中 \(\text{Concatenate}(n)\) 是將 \(1 \sim n\) 所有正整數 順序連線起來得到的數。

例如,\(n = 13\)\(\text{Concatenate}(n) = 12345678910111213\)。小 C 想了大半天終於意識到這是一道不可能手算出來的題目,於是他只好向你求助,希望你能編寫一個程式幫他解決這個問題。

對於 \(100\%\) 的資料,\(1\le n \le 10^{18}\)\(1\le m \le 10^9\)

題解:

這裡我們先用遞推的形式,大概能寫出來,當i是個位數的時候,我們只需要\(F_{i-1}\)乘10,再加上\(i\)就好了,但是我們發現當i是兩位數的時候,我們就是乘以100而不是10,以此類推,當我們是k位數時,我們需要乘\(10^k\),所以這裡我們沒有辦法像之前那樣,找到轉移矩陣一口氣直接乘到最後就好了,我們可以先寫出來當i是個位數時的轉移矩陣,然後需要每當位數+1的時候我們的轉移矩陣也更新,將10變成100,再變成1000...,就可以了,時間複雜度會比之前的log多一個lg出來。
(順便一提,這個題需要開unsigned long long,並且不要使用std的pow函式,返回型別是int,還有判斷位數的時候不要使用log10函式,要手寫!!不要問我怎麼知道的qwq)

int kuai(int a, int b) {
    int tem = 1;
    while (b) {
        if (b % 2) tem = tem * a;
        a = a * a, b /= 2;
    }
    return tem;
}
int fan(int x) {
    int cnt = 0;
    while (x) {
        x /= 10;
        cnt++;
    }
    return cnt;
}
signed main() {
    read(), kuaidu();
    // cout << log10(124551) << endl;
    cin >> n >> mod;
    Mat base(3, 3);
    base[1][1] = 10, base[2][1] = 1, base[2][2] = 1, base[3][2] = 1, base[3][3] = 1;
    Mat dp(1, 3);
    dp[1][1] = 0, dp[1][2] = 1, dp[1][3] = 1;

    int k = fan(n) - 1;
    int las = n - kuai(10, k);
    int q = 1;
    while (1) {
        base[1][1] = q * 10 % mod;
        q *= 10;
        if (q <= n) { dp = dp * Mat::kuai(base, q - q / 10); }
        else break;
    }
    dp = dp * Mat::kuai(base, las + 1);
    cout << dp[1][1] << endl;

    return 0;
}

關於矩陣加速的內容先寫這麼多,之後還會再更一篇矩陣的題,是2024杭電7的迴圈圖。就是因為這個題才讓我開始補矩陣的-_-

相關文章