矩陣乘法

時間  2020-06-12
標籤 矩陣 乘法 欄目 應用數學 简体版
原文   原文鏈接

目錄

矩陣乘法

1. 算法分析

利用快速冪的方法來優化矩陣的乘法，使得計算矩陣A(N*N)的M次方的時間優化到O(N3logM)c++

經常使用技巧
通常若是可以把式子寫成 Kn = Kn-1+t，那就能使用矩陣快速冪處理，設Fn=[fn, fn+1, kn], Fn+1=[fn+1, fn+2, kn+1],那麼Fn+1=FnA (A爲矩陣)算法

2. 板子

  1. 計算ax
#include <bits/stdc++.h>

using namespace std;

typedef long long LL;
const int N = 110, mod = 1e9 + 7;
int n;

// 定義一個矩陣
struct mat
{
    int m[N][N];
}unit;

// 定義矩陣乘法
mat operator * (mat a, mat b)
{
    mat res;
    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;
            res.m[i][j] = x % mod;
        }
    return res;
}

// 初始化單位陣
void init_unit()
{
    for (int i = 0; i < N; ++i)
        unit.m[i][i] = 1;
    return ;
}

// 矩陣快速冪
mat pow_mat(mat a, LL n)
{
    mat res = unit;  // 初始爲單位陣
    while (n)
    {
        if (n & 1) res = res * a;
        a = a * a;
        n >>= 1;
    }
    return res;
}

int main()
{
    LL x;
    init_unit();
    while (cin >> n >> x)
    {
        mat a;
        for (int i = 0; i < n; ++i)
            for (int j = 0; j < n; ++j)
                cin >> a.m[i][j];
        a = pow_mat(a, x);
        for (int i = 0; i < n; ++i)
            for (int j = 0; j < n; ++j)
                if (j + 1 == n) cout << a.m[i][j] << endl;
                else cout << a.m[i][j] << " ";
    }
    return 0;
}
  1. 計算斐波那契數列
#include <bits/stdc++.h>

using namespace std;

typedef long long LL;
const int N = 4, mod = 10000 ;
int n;

// 定義一個矩陣
struct mat
{
    int m[N][N];
}unit;

// 定義矩陣乘法
mat operator * (mat a, mat b)
{
    mat res;
    LL x;
    for (int i = 0; i < 2; ++i)
        for (int j = 0; j < 2; ++j)
        {
            x = 0;
            for (int k = 0; k < 2; ++k)
                x += (LL)a.m[i][k] * b.m[k][j] % mod;
            res.m[i][j] = x % mod;
        }
    return res;
}

// 初始化單位陣
void init_unit()
{
    for (int i = 0; i < N; ++i)
        unit.m[i][i] = 1;
    return ;
}

// 矩陣快速冪
mat pow_mat(mat a, LL n)
{
    mat res = unit;  // 初始爲單位陣
    while (n)
    {
        if (n & 1) res = res * a;
        a = a * a;
        n >>= 1;
    }
    return res;
}

int main()
{
    init_unit();
    while (cin >> n && n != -1)
    {
        if (!n) {
            cout << 0 << endl;
            continue;
        }
        mat a;
        a.m[0][0] = 0, a.m[0][1] = 1;
        a.m[1][0] = 1, a.m[1][1] = 1;
        a = pow_mat(a, n - 1);
        LL res = (a.m[1][0] + 0ll + a.m[0][0]) % mod;
        cout << res << endl;
    }
    return 0;
}

3. 例題

3.1 直接處理矩陣運算

acwing225矩陣冪求和
給定n×n矩陣A和正整數k，求和S=A+A2+A3+…+Ak。結果輸出時每一個元素都須要mod m
1≤n≤30,1≤k≤1e9,1≤m<1e4優化

/*
分治的思想
Sn = a+a^2+a^3+...a^n
若是n是奇數那麼:
Sn = (a + a^2 + ... a^(n/2))*(a^(n/2) + 1) + a^n
若是n是偶數那麼:
Sn = (a + a^2 + ... a^(n/2))*(a^(n/2) + 1)
基於這個思想，能夠不斷把前項進行分治處理
這樣時間複雜度爲:O(((logk)^2) * n^3 )
*/

#include <bits/stdc++.h>

using namespace std;

typedef long long LL;
const int N = 31;
int n, k, mod;

// 定義一個矩陣
struct mat
{
    int m[N][N];
}unit;
mat a;

// 定義矩陣乘法
mat operator * (mat a, mat b)
{
    mat res;
    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;
            res.m[i][j] = x % mod;
        }
    return res;
}

mat operator + (mat a, mat b) {
    mat res;
    for (int i = 0 ; i < n; ++i) {
        for (int j = 0 ; j < n; ++j) {
            res.m[i][j] = ((LL)a.m[i][j] + b.m[i][j]) % mod;
        }
    }
    return res;
}

// 初始化單位陣
void init_unit()
{
    for (int i = 0; i < N; ++i)
        unit.m[i][i] = 1;
    return ;
}

// 矩陣快速冪
mat pow_mat(mat a, LL n)
{
    mat res = unit;  // 初始爲單位陣
    while (n)
    {
        if (n & 1) res = res * a;
        a = a * a;
        n >>= 1;
    }
    return res;
}

mat dfs(int u) {
    if (u == 1) return a;  // 分治到1
    if (u & 1) return dfs(u / 2) * (unit + pow_mat(a, u / 2)) + pow_mat(a, u);  // 奇數
    else return dfs(u / 2) * (unit + pow_mat(a, u / 2));  // 偶數
}

int main()
{
    init_unit();
    cin >> n >> k >> mod;
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j)
            cin >> a.m[i][j];
    mat res = dfs(k);
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j)
            cout << res.m[i][j] << " ";
        cout << endl;
    }
    return 0;
}

acwing226 233矩陣
假設咱們有一個名爲233矩陣的矩陣。
在第一行，它將包含233,2333,23333 …（這意味着a0,1=233，a0,2=2333，a0,3=23333…）。
此外，在233矩陣中，知足ai,j=ai−1,j+ai,j−1（i,j≠0）。
如今給定a1,0,a2,0,…,an,0，請求出在233矩陣中an,m的值。
1≤n≤10,1≤m≤1e9,ai,0~int
Screenshot_20200531_100040_com.huawei.notepad.jpgspa

/*
233矩陣有特殊的特色
記x0=23
則y0=233=10x0+3
y1=10x0+x1+3
y2=10x0+x1+x2+3
y3=10x0+x1+x2+x3+3
...
yn=10x0+x1+x2+...+xn+3
那麼能夠矩陣快速冪求出每一個y
*/
#include <bits/stdc++.h>

using namespace std;

typedef long long LL;
const int N = 20, mod = 10000007;
int n, m;
int x[N];

// 定義一個矩陣
struct mat
{
    int m[N][N];
}unit;

// 定義矩陣乘法
mat operator * (mat a, mat b)
{
    mat res;
    LL x;
    for (int i = 0; i < n + 2; ++i)
        for (int j = 0; j < n + 2; ++j)
        {
            x = 0;
            for (int k = 0; k < n + 2; ++k)
                x += (LL)a.m[i][k] * b.m[k][j] % mod;
            res.m[i][j] = x % mod;
        }
    return res;
}

// 初始化單位陣
void init_unit()
{
    for (int i = 0; i < N; ++i)
        unit.m[i][i] = 1;
    return ;
}

// 矩陣快速冪
mat pow_mat(mat a, LL n)
{
    mat res = unit;  // 初始爲單位陣
    while (n)
    {
        if (n & 1) res = res * a;
        a = a * a;
        n >>= 1;
    }
    return res;
}

int main()
{
    init_unit();
    while (scanf("%d%d", &n, &m) != EOF)
    {
        for (int i = 1; i <= n; ++i) scanf("%d", &x[i]);
        mat a;
        
        // 構造矩陣
        for (int i = 0; i < n + 2; ++i){
            if (i == 0) a.m[0][i] = 10;
            else if (i == 1) a.m[0][i] = 3;
            else a.m[0][i] = 0;
        }
        for (int i = 0; i < n + 2; ++i) {
            if (i == 1) a.m[1][i] = 1;
            else a.m[1][i] = 0;
        }
        for (int i = 2; i < n + 2; ++i) {
            for (int j = 0; j < n + 2; ++j) {
                if (j == 0) a.m[i][j] = 10;
                else if (j == 1) a.m[i][j] = 3;
                else if (j <= i) a.m[i][j] = 1;
                else a.m[i][j] = 0;
            }
        }
        
        // 計算矩陣的m次冪
        a = pow_mat(a, m);
        
        // 計算anm
        LL res = (a.m[n + 1][0] * 23 + 0ll + a.m[n + 1][1] * 1) % mod;
        for (int i = 1; i <= n; ++i) res = (res + 0ll + a.m[n + 1][i + 1] * 1ll * x[i] % mod) % mod;
        cout << res << endl; 
    }
    return 0;
}

3.2 斐波那契數列+矩陣快速冪

acwing1303斐波那契前n項和
求斐波那契數列fn的前n項和Sn mod m
n~2e9, m~1e9 + 10code

/*
斐波那契數列能夠使用矩陣快速冪來處理 
設Fn={fn, fn+1, Sn}, Fn+1={fn+1, fn+2, Sn+1}
Fn+1=Fn * {0 1 0  = Fn * A = F1*A^n = {1, 1, 1} * A^n
           1 1 1
           0 0 1}
則Sn=Fn[3] = {1, 1, 1} * A ^(n - 1)
*/
#include <bits/stdc++.h>

using namespace std;

typedef long long LL;
const int N = 110;
int n, m;

// 定義一個矩陣
struct mat
{
    int m[N][N];
}unit;

// 定義矩陣乘法
mat operator * (mat a, mat b)
{
    mat res;
    LL x;
    for (int i = 0; i < 3; ++i)
        for (int j = 0; j < 3; ++j)
        {
            x = 0;
            for (int k = 0; k < 3; ++k)
                x += (LL)a.m[i][k] * b.m[k][j] % m;
            res.m[i][j] = x % m;
        }
    return res;
}

// 初始化單位陣
void init_unit()
{
    for (int i = 0; i < 3; ++i)
        unit.m[i][i] = 1;
    return ;
}

// 矩陣快速冪
mat pow_mat(mat a, LL n)
{
    mat res = unit;  // 初始爲單位陣
    while (n)
    {
        if (n & 1) res = res * a;
        a = a * a;
        n >>= 1;
    }
    return res;
}

int main()
{
    init_unit();
    cin >> n >> m;
    int f[] = {1, 1, 1};  // F1
    mat a;
    int matrix[] = {0, 1, 0, 1, 1, 1, 0, 0, 1};
    for (int i = 0; i < 3; ++i)
        for (int j = 0; j < 3; ++j)
            a.m[i][j] = matrix[i * 3 + j];
    mat res = pow_mat(a, n - 1);  // A^(n - 1)
    
    // 計算Fn[2]
    int ans = 0;
    for (int i = 0; i < 3; ++i)
        ans = (ans + res.m[i][2] * f[i]) % m;
    cout << ans << endl;
    return 0;
}

acwing1304佳佳的斐波那契
用 T(n)=(F1+2F2+3F3+…+nFn)modm 表示 Fibonacci 數列前 n 項變形後的和 modm 的值。
如今佳佳告訴你了一個 n 和 m，請求出 T(n) 的值。
n,m~intci

/*

本題Sn = F1 + F2 + ... + Fn, Tn= F1 + 2F2 + ... + nFn,其中Sn爲n級別，Tn爲n方級別
那麼考慮處理的時候Pn=nSn-Tn=(n-1)F1+(n-2)F2+...+Fn-1+0, Pn-1=(n-1)Sn-1 - Tn-1 = (n-2)F1 + (n-3)F2 + ... + Fn-2+0+0
則Pn=Pn-1+Sn-1，而Sn=Sn-1+Fn, Fn=Fn-1+Fn-2，當出現這種沒有係數的式子時就能夠轉化爲矩陣乘法:
則設Kn={fn, fn+1, Sn, Pn}, Kn-1={fn-1, fn, Sn-1, Pn-1}
Kn = Kn-1{ 0 1 0 0 = Kn-1 * A = K1 * A^(n - 1) = {1 ,1 ,1, 0} * A ^ (n - 1)
           1 1 1 0
           0 0 1 1
           0 0 0 1}
記x = {1 ,1 ,1, 0} * A ^ (n - 1)
那麼要求Tn=nSn - Pn=n*x[2] - x[3]
*/
#include <bits/stdc++.h>

using namespace std;

typedef long long LL;
const int N = 110;
int n, m;

// 定義一個矩陣
struct mat
{
    int m[N][N];
}unit;

// 定義矩陣乘法
mat operator * (mat a, mat b)
{
    mat res;
    LL x;
    for (int i = 0; i < 4; ++i)
        for (int j = 0; j < 4; ++j)
        {
            x = 0;
            for (int k = 0; k < 4; ++k)
                x += (LL)a.m[i][k] * b.m[k][j] % m;
            res.m[i][j] = x % m;
        }
    return res;
}

// 初始化單位陣
void init_unit()
{
    for (int i = 0; i < 4; ++i)
        unit.m[i][i] = 1;
    return ;
}

// 矩陣快速冪
mat pow_mat(mat a, LL n)
{
    mat res = unit;  // 初始爲單位陣
    while (n)
    {
        if (n & 1) res = res * a;
        a = a * a;
        n >>= 1;
    }
    return res;
}

int main()
{
    init_unit();
    cin >> n >> m;
    int f[] = {1, 1, 1, 0};
    mat a;
    int matrix[] = {0, 1, 0, 0, 1, 1, 1, 0, 0 ,0, 1, 1, 0, 0, 0, 1};
    for (int i = 0; i < 4; ++i)
        for (int j = 0; j < 4; ++j)
            a.m[i][j] = matrix[i * 4 + j];
    mat res = pow_mat(a, n - 1);
    int Sn = 0, Pn = 0;
    for (int i = 0; i < 4; ++i)
        Sn = (Sn + res.m[i][2] * f[i]) % m, Pn = (Pn + res.m[i][3] * f[i]) % m;
    // cout << (LL)n * Sn % m << endl << Pn % m << endl;
    cout << (((LL)n * Sn % m - Pn % m + m) % m + m ) % m<< endl;
    return 0;
}
相關文章
相關標籤/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程