抱歉,您的浏览器无法访问本站
本页面需要浏览器支持(启用)JavaScript
了解详情 >

本文正在施工中,有缺失的内容请谅解。

线性代数是代数学的一个分支,主要处理线性关系问题。OI 中用到的相关知识并不多,主要是矩阵乘法与高斯消元等,本文将简单介绍这些内容。

矩阵乘法

一个 n×mn\times m 的矩阵可以看作是一个 n×mn\times m 的二维数组。矩阵加法和减法就是把相对应的数加减,而矩阵乘法则很有意思,拥有几何意义,这里不做展开。

矩阵乘法

AAn×mn\times m 的矩阵,BBm×pm\times p 的矩阵,那么 C=A×BC=A\times Bn×pn\times p 的矩阵,并且:

Ci,j=k=1Ai,kBk,jC_{i,j}=\sum_{k=1}A_{i,k}B_{k,j}

模板,代码如下:

for (int i = 0; i < n; ++i)
    for (int k = 0; k < n; ++k) {
        int r = a[i][k];
        for (int j = 0; j < n; ++j)
            c[i][j] += r * b[k][j];
    }

注意循环顺序!这样的顺序对答案不会造成影响,但是缓存访问是连续的,因此跑得较快。

矩阵快速幂

矩阵乘法不满足交换律,但是满足结合律(也满足分配律),因此快速幂的计算方式可以直接移用。注意矩阵乘法有一个“单位矩阵”的概念,类似于数的乘法中“1”。

模板,代码入下:

查看代码
struct Matrix {
    int a[105][105];
    Matrix() { memset(a, 0, sizeof(a)); };
    void reset(void) { // 构建单位矩阵
        for (int i = 1; i <= n; ++i)
            a[i][i] = 1;
    }
};

Matrix operator * (const Matrix &x, const Matrix &y) {
    Matrix z;
    for (int i = 1; i <= n; ++i)
        for (int k = 1; k <= n; ++k) {
            int r = x.a[i][k];
            for (int j = 1; j <= n; ++j)
                z.a[i][j] = (z.a[i][j] + 1ll * r * y.a[k][j]) % MOD;
        }
    return z;
}

Matrix poww(Matrix a, i64 b) {
    Matrix res; res.reset();
    while (b) {
        if (b & 1) res = res * a;
        a = a * a;
        b >>= 1;
    }
    return res;
}

矩阵与递推

矩阵乘法的一大作用就是加速递推。

概述

递推都有一个递推式子,我们可以把这个式子中的变量写进矩阵,然后再乘上一个矩阵进行变换。由于递推式必变,因此每次乘上的矩阵都是相同的,可以使用矩阵快速幂来进行计算。

[Luogu P1962] 斐波那契数列

Portal.

要将 [Fn1Fn2]\begin{bmatrix}F_{n-1}\\F_{n-2}\end{bmatrix} 变换为 [FnFn1]\begin{bmatrix}F_{n}\\F_{n-1}\end{bmatrix},利用待定系数法可以求得需要乘上 A=[1110]\mathbf{A}=\begin{bmatrix}1 & 1\\1 & 0\end{bmatrix},即 [FnFn1]=[1110][Fn1Fn2]\begin{bmatrix}F_{n}\\F_{n-1}\end{bmatrix}=\begin{bmatrix}1 & 1\\1 & 0\end{bmatrix}\begin{bmatrix}F_{n-1}\\F_{n-2}\end{bmatrix}。要多次变换的话,由于矩阵乘法有结合律,利用矩阵快速幂就可以搞定,最后在乘上初始 FF 矩阵。

查看代码
#include <iostream>
#include <cstdio>
#include <cstring>

using namespace std;
typedef long long i64;
const int MOD = 1000000007;

struct Matrix {
    int a[3][3];
    Matrix() { memset(a, 0, sizeof(a)); };
    void reset(void) {
        for (int i = 1; i <= 2; ++i)
            a[i][i] = 1;
    }
    Matrix friend operator * (const Matrix &a, const Matrix &b) {
        Matrix c;
        for (int i = 1; i <= 2; ++i)
            for (int k = 1; k <= 2; ++k) {
                int r = a.a[i][k];
                for (int j = 1; j <= 2; ++j)
                    c.a[i][j] = (c.a[i][j] + 1ll * r * b .a[k][j]) % MOD;
            }
        return c;
    }
};

Matrix poww(Matrix T, i64 b) {
    Matrix res; res.reset();
    while (b) {
        if (b & 1) res = res * T;
        T = T * T;
        b >>= 1;
    }
    return res;
}

int main(void) {
    i64 n; cin >> n;
    if (n <= 2) return puts("1"), 0;
    Matrix a, f;
    a.a[1][1] = a.a[1][2] = a.a[2][1] = 1;
    f.a[1][1] = 1; // f 记录 F[1], F[0]
    a = poww(a, n - 1); // 从 n = 2 变换,所以 -1
    f = a * f;
    cout << f.a[1][1] << '\n';
    return 0;
}

[Luogu P1939] 【模板】矩阵加速(数列)

Portal.

方法是相同的。

查看代码
#include <iostream>
#include <cstdio>
#include <cstring>

using namespace std;
const int MOD = 1000000007;
const int N = 3;

struct Matrix {
    int a[4][4];
    Matrix(void) { memset(a, 0, sizeof(a)); }
    void build(void) {
        for (int i = 1; i <= N; ++i)
            a[i][i] = 1;
    }
    Matrix friend operator * (const Matrix &a, const Matrix &b) {
        Matrix c;
        for (int i = 1; i <= N; ++i)
            for (int k = 1; k <= N; ++k) {
                int r = a.a[i][k];
                for (int j = 1; j <= N; ++j)
                    c.a[i][j] = (c.a[i][j] + 1ll * r * b.a[k][j]) % MOD;
            }
        return c;
    }
};

Matrix poww(Matrix a, int b) {
    Matrix res; res.build();
    while (b) {
        if (b & 1) res = res * a;
        a = a * a;
        b >>= 1;
    }
    return res;
}

int main(void) {
    int T; cin >> T;
    while (T--) {
        int n; cin >> n;
        if (n <= 3) {
            puts("1");
            continue;
        }
        Matrix a, f;
        a.a[1][1] = a.a[1][3] = a.a[2][1] = a.a[3][2] = 1;
        f.a[1][1] = f.a[2][1] = f.a[3][1] = 1;
        a = poww(a, n - 3);
        f = a * f;
        cout << f.a[1][1] << '\n';
    }
    return 0;
}

[Luogu P1349] 广义斐波那契数列

Portal.

只是我们的变换矩阵不一样了。

查看代码
#include <iostream>
#include <cstdio>
#include <cstring>

using namespace std;
const int N = 2;

int m;
struct Matrix {
    int a[3][3];
    Matrix(void) { memset(a, 0, sizeof(a)); }
    void build(void) {
        for (int i = 1; i <= N; ++i)
            a[i][i] = 1;
    }
    Matrix friend operator * (const Matrix &a, const Matrix &b) {
        Matrix c;
        for (int i = 1; i <= N; ++i)
            for (int k = 1; k <= N; ++k) {
                int r = a.a[i][k];
                for (int j = 1; j <= N; ++j)
                    c.a[i][j] = (c.a[i][j] + 1ll * r * b.a[k][j]) % m;
            }
        return c;
    }
};

Matrix poww(Matrix a, int b) {
    Matrix res; res.build();
    while (b) {
        if (b & 1) res = res * a;
        a = a * a;
        b >>= 1;
    }
    return res;
}

int main(void) {
    int p, q, a1, a2, n;
    cin >> p >> q >> a1 >> a2 >> n >> m;
    if (n <= 2) return (cout << (n == 1 ? a1 : a2) << '\n', 0);
    Matrix a, f;
    a.a[1][1] = p, a.a[1][2] = q, a.a[2][1] = 1;
    f.a[1][1] = a2, f.a[2][1] = a1;
    a = poww(a, n - 2);
    f = a * f;
    cout << f.a[1][1] << '\n';
    return 0;
}

高斯消元

相关知识点

LGV 引理

矩阵树定理

BEST 定理

线性基

Problemset

主要是应用。

线性代数基础

一些简单题。用到了概念和矩阵乘法。

[NOI2012] 随机数生成器

Portal.

这样:

[Fnc]=[a101][Fn1c]\begin{bmatrix}F_n\\c\end{bmatrix}=\begin{bmatrix}a&1\\0&1\end{bmatrix}\begin{bmatrix}F_{n-1}\\c\end{bmatrix}

查看代码
#include <iostream>
#include <cstdio>
#include <cstring>

using namespace std;
typedef long long i64;
typedef __int128_t LL;

i64 m, a, c, x0, n, g;

struct Matrix {
    i64 a[2][2];
    Matrix() { memset(a, 0, sizeof(a)); }
    void build(void) { a[0][0] = a[1][1] = 1; }
    Matrix friend operator * (const Matrix &a, const Matrix &b) {
        Matrix c;
        for (int i = 0; i < 2; ++i)
            for (int k = 0; k < 2; ++k) {
                i64 r = a.a[i][k];
                for (int j = 0; j < 2; ++j)
                    c.a[i][j] = (c.a[i][j] + (LL)r * b.a[k][j]) % m;
            }
        return c;
    }
};

Matrix poww(Matrix a, i64 b) {
    Matrix res; res.build();
    while (b) {
        if (b & 1) res = res * a;
        a = a * a;
        b >>= 1;
    }
    return res;
}

int main(void) {
    cin >> m >> a >> c >> x0 >> n >> g;
    Matrix A, F;
    A.a[0][0] = a; A.a[0][1] = A.a[1][1] = 1;
    F.a[0][0] = x0, F.a[1][0] = c;
    A = poww(A, n);
    F = A * F;
    cout << F.a[0][0] % g << endl;
    return 0;
}

[NOI2013] 向量内积

Portal.

发现模数是 2233,因此考虑乱搞。

维护向量的前缀和,让一个向量乘以它前面的前缀和向量。k=2k=2 时如果不为 (i1)mod2(i-1)\bmod 2 说明有一个乘积为 00,暴力找即可。

k=3k=3 时维护前缀平方和可以达到一样的效果。

随机化做几轮即可。

查看代码
#include <bits/stdc++.h>
using namespace std;

int n, m, k; 
int a[100005][105], id[100005], b[105], c[105][105]; 
mt19937 Rand(time(0)); 

inline bool check(int x, int y) {
    int r = 0; 
    for (int i = 1; i <= m; ++i) r += a[x][i] * a[y][i]; 
    return r % k == 0; 
}

inline int calc(int x) {
    int ans = 0; 
    if (k == 2) {
        for (int i = 1; i <= m; ++i) {
            ans = (ans + b[i] * a[x][i]) % k;
            b[i] = (b[i] + a[x][i]) % k; 
        }
    } else {
        for (int i = 1; i <= m; ++i)
            for (int j = 1; j <= m; ++j) {
                ans = (ans + c[i][j] * a[x][i] * a[x][j]) % k; 
                c[i][j] = (c[i][j] + a[x][i] * a[x][j]) % k; 
            }
    }
    return ans; 
}

int main(void) {
    scanf("%d%d%d", &n, &m, &k); 
    for (int i = 1; i <= n; ++i) for (int j = 1; j <= m; ++j) scanf("%d", &a[i][j]), a[i][j] %= k; 
    for (int i = 1; i <= n; ++i) id[i] = i; 
    for (int T = 1; T <= 7; ++T) {
        shuffle(id + 1, id + n + 1, Rand);
        memset(b, 0, sizeof b); memset(c, 0, sizeof c); 
        for (int i = 1; i <= n; ++i) if (calc(id[i]) != (i - 1) % k) // 有解
            for (int j = 1; j < i; ++j)
                if (check(id[i], id[j])) {
                    if (id[i] > id[j]) swap(id[i], id[j]); 
                    return printf("%d %d\n", id[i], id[j]), 0; 
                }
    }
    puts("-1 -1"); return 0;
}

高斯消元

评论

若无法加载,请尝试刷新,欢迎讨论、交流和提出意见,支持 Markdown 与 LaTeX 语法(公式与文字间必须有空格)!