1. 程式人生 > >矩陣快速冪(模版) 快速乘 + 快速冪 + 矩陣快速冪

矩陣快速冪(模版) 快速乘 + 快速冪 + 矩陣快速冪

矩陣快速冪:首先前置技能:  快速冪 + 矩陣乘法。

1  快速冪

1.1 快速乘法

1.1.1

引用自2009年國家集訓隊論文,駱可強:《論程式底層優化的一些方法與技巧》 (膜膜膜)

可以根據需要換成uLL   (unsigned long long) 

// O(1)

#include <stdio.h>
typedef long long LL;
inline LL mult_mod(LL a, LL b, LL m){ // a * b % m
    a = (a % m + m) % m;
    b = (b % m + m) % m;
    LL ans = a * b - (LL)((long double) a * b / m + 0.5) * m;
    return ans < 0? ans + m : ans;
}

int main(){
    LL x, y, m;
    while(~scanf("%lld%lld%lld", &x, &y, &m)){
        printf("%lld\n", mult_mod(x, y, m));
    }
    return 0;
}

 1.1.2  O(log b)

原理:

 由於計算機底層設計的原因,做加法往往比乘法快的多,因此將乘法轉換為加法計算將會大大提高(大數,比較小的數也沒必要)乘法運算的速度,除此之外,當我們計算a*b%mod的時候,往往較大的數計算a*b會超出long long int的範圍,這個時候使用快速乘法方法也能解決上述問題.    快速乘法的原理就是利用乘法分配率來將a*b轉化為多個式子相加的形式求解(注意這時使用乘法分配率的時候後面的一個乘數轉化為二進位制的形式計算).舉個栗子    20*14 = 20*(1110)2 = 20*(2^3)*1 + 20*(2^2)*1+20*(2^1)*1+20*(2^0)*0 = 160+80+40=280.    這樣計算了4次 而直接算的話 需要14次計算哦

typedef long long LL;
LL multi(LL a, LL b, LL m){
    // a * b % m
    LL ret = 0;
    while(b > 0){
        if(b & 1){
            ret = (ret + a) % m;
        }
        b >>= 1;
        a = (a << 1) % m;
    }
    return ret;
}

 關於用哪個 emmm 大家自己決定吧, 一般都是用log b那個 

1.2  快速冪

其實快速冪的思想和log那個快速乘的思想差不多 也是將b轉為二進位制  不多贅述了

用到快速乘是因為 當a 和 b 太大的時候 用快速乘優化

#include <stdio.h>
typedef long long LL;
LL multi(LL a, LL b, LL m){
    // a * b % m
    LL ret = 0;
    while(b > 0){
        if(b & 1){
            ret = (ret + a) % m;
        }
        b >>= 1;
        a = (a << 1) % m;
    }
    return ret;
}

LL quick_mod(LL a, LL b, LL m){
    LL ans = 1;
    while(b){
        if(b & 1){
            ans = multi(ans,a,m);
            b--;
        }
        b /= 2;
        a = multi(a, a, m);
    }
    return ans;
}

 2 矩陣乘法

2.1矩陣乘法

矩陣屬於線性代數的東西。 受教於MR.ZZZ, 嚶嚶嚶 趙哥教了我好多東西。qaq

2.1.1 矩陣

由 m × n 個數aij排成的m行n列的數表稱為m行n列的矩陣,簡稱m × n矩陣。記作:

\\ \\ A\quad =\quad \begin{bmatrix} a11 & a12 & a13 & ... & a1n \\ a21 & a21 & a23 & ... & a2n \\ a31 & a32 & a33 & ... & a3n \\ ... & ... & ... & ... & ... \\ am1 & am2 & am3 & ... & amn \end{bmatrix}\\ \\ 

2.1.2 矩陣乘法

前提: 矩陣A 乘 矩陣B  可以想乘的前提是 A的列數 == B的行數

兩個矩陣的乘法僅當第一個矩陣A的列數和另一個矩陣B的行數相等時才能定義

矩陣A (m * a)  * 矩陣B(a * n) = 矩陣C(m * n) 

\begin{bmatrix} a11 & a12 \\ a21 & a22 \end{bmatrix}\begin{bmatrix} b11 & b12 \\ b21 & b22 \end{bmatrix}\quad =\quad \begin{bmatrix} a11*b11\quad +\quad a12*b21 & a11*b12\quad +\quad a12*b22 \\ a21*b11\quad +\quad a22*b21 & a21*b12\quad +\quad a22*b22 \end{bmatrix}\\ \\ C(i,\quad j)\quad =\quad \sum _{ r\quad =\quad 1 }^{ n }{ A(i,r)\quad *\quad B(r,\quad j) } \\

矩陣乘法如何用程式碼實現呢?

#include <stdio.h>
const int maxn = 1000 + 10;
int a[maxn][maxn], b[maxn][maxn];
long long c[maxn][maxn];
int main(){
    int n;
    scanf("%d", &n);
    for(int i = 1; i <= n; ++i){
        for(int j = 1; j <= n; ++j){
            scanf("%d", &a[i][j]);
        }
    }
    for(int i = 1; i <= n; ++i){
        for(int j = 1; j <= n; ++j){
            scanf("%d", &b[i][j]);
        }
    }
    for(int i = 1; i <= n; i++){
        for(int j = 1; j <= n; j++){
            for(int k = 1; k <= n; k++){
                c[i][j] += a[i][k] * b[k][j];
            }
            if(j == n){
                printf("%lld\n", c[i][j]);
            }
            else{
                printf("%lld ", c[i][j]);
            }
        }
    }
    return 0;
}

2.2 單位矩陣

矩陣的乘法中,有一種矩陣起著特殊的作用,如同數的乘法中的1,這種矩陣被稱為單位矩陣。它是個方陣,從左上角到右下角的對角線(稱為主對角線)上的元素均為1。除此以外全都為0。

單位矩陣的乘法: 矩陣A * 單位矩陣 = 矩陣A

3 矩陣快速冪

3.1 有了前兩個前置技能後, 那麼矩陣快速冪就很好理解了。 我們可以把矩陣看作一個特殊的數字,然後去用快速冪的思想去寫,最開始的ans為單位矩陣(單位矩陣的特性)。我們重寫一下乘法就好了,具體看程式碼吧。

//
//  1113 - 矩陣快速冪.cpp
//  數論
//
//  Created by Terry on 2018/9/27.
//  Copyright © 2018年 Terry. All rights reserved.
//
//  1: 單位矩陣 * 矩陣A = 矩陣A
//  2: 矩陣快速冪 可以對應 快速冪去理解 就是重寫一下乘法  a
//      把快速冪裡面的數字乘法 重寫成 矩陣乘法
#include <stdio.h>
int read(){
    int x = 0, f = 1;
    char ch=getchar();
    while(ch < '0' || ch > '9'){if(ch=='-') f = -1; ch = getchar();}
    while(ch >= '0' && ch <= '9'){x = x * 10 + ch - '0'; ch = getchar();}
    return x * f;
}
typedef long long LL;
const int mod = 1e9 + 7;
const int maxn = 100 + 10;
struct Matrix{
    int m[maxn][maxn];
}unit;
int n;
Matrix operator * (Matrix a, Matrix b){
    Matrix ret;
    LL x;
    for(int i = 0; i < n; i++){
        for(int j = 0; j < n; j++){
            x = 0;
            for(int k = 0; k < n; ++k){
                x += ((LL)a.m[i][k] * b.m[k][j]) % mod;
            }
            ret.m[i][j] = x % mod;
        }
    }
    return ret;
}
void init_unit(){
    // 單位矩陣
    for(int i  = 0; i < maxn; i++){
        unit.m[i][i] = 1;
    }
}
Matrix pow_mat(Matrix a, LL n){
    Matrix ans = unit;
    while (n) {
        if(n & 1){
            ans = ans * a;
        }
        a = a * a;
        n >>= 1;
    }
    return ans;
}
int main(){
    init_unit();
    n = read();
    LL x = read();
    Matrix a;
    for(int i = 0; i < n; i++){
        for(int j = 0; j < n; j++){
            a.m[i][j] = read();
        }
    }
    a = pow_mat(a, x);
    for(int i = 0; i < n; i++){
        for(int j = 0; j < n; j++){
            if(j == n - 1){
                printf("%d\n", a.m[i][j]);
            }
            else{
                printf("%d ", a.m[i][j]);
            }
        }
    }
    return 0;
}