线性代数学习_线性代数什么时候x几等于自由变元-CSDN博客

本文链接：https://blog.csdn.net/g1342522389/article/details/84928327

这周学习了线性代数，感觉很厉害的样子。。。

高斯消元

对于一组多项式方程（增广矩阵中，x[i, n+1]表示式子的值；x[i,j]表示第i个方程第j项的系数，在这里，增广矩阵可能不一定是n个，可能多可能少；opt表示运算规则）：
(x[1,1]*a[1]) opt (x[1,2]*a[2]) opt … opt (x[1,n]*a[n])=x[1, n+1]
(x[2,1]*a[1]) opt (x[2,2]*a[2]) opt … opt (x[2,n]*a[n])=x[2, n+1]
…
(x[n,1]*a[1]) opt (x[n,2]*a[2]) opt … opt (x[n,n]*a[n])=x[n, n+1]
常见的opt有+和异或
消元：我们对于每个方程，假设为i，且x[i][now]不为0（now是当前被消元的变元，更一般的，只要能支持消掉后边方程的系数即可），然后对所有的j>i的方程，我们将所有x[j][now]不为0（或者需要消元的）用一种方式消掉。

回代：最后我们得到的是一个x[i][i]均不为0的方程（除非无解或此方程无意义）的倒三角矩阵。然后根据性质3将当前方程的其它j>i的a[j]和系数根据性质3转移到右式消去即可（此时a[j]一定是求出来的）

无解的情况：当存在一个矩阵使得左式=0而右式!=0那么无解。
自由变元：自由变元就是当这些未知量一旦确定，整个方程就确定了。但是这些量是未知的。（例如x+y=5，自由变元就是1，因为无论是x还是y确定，另一个就能唯一确定），而答案要求的是方案，那么显然因为自由变元是可以随便赋值的，而如果这些值只有2个，开和不开，那么方案数就是2^自由变元。
仅当n个不同的方程（就是无论怎么通过其它方程都不会将这两个方程变成一样）才能确定n个解。那么我们如果只确定了x个方程，那么自由变元的数量就是n-x。（这个x可以轻易得到，因为在高斯消元过程中，会将所有不同的方程消元，因为消元会将相同的方程消成这个样子：0=0。所以就能得到x了。
复杂度O（n^3）
代码

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
#define f(i,a,b) for(int i=a;i<=b;i++)
const double eps=0.000000001;
const int N=101;
double a[N][N+1];
int n;
int main(){
  cin>>n;
  f(i,1,n){
      f(j,1,n+1) cin>>a[i][j];
  }
  f(i,1,n){
      int top=i;
      f(j,i,n) if(fabs(a[j][i]-a[top][i])<=eps) top=j;
      f(j,1,n+1) swap(a[i][j],a[top][j]);
      if(fabs(a[i][i])<=eps){
          cout<<"No Solution"<<endl;
          return 0;
      }
      f(j,i+1,n+1) a[i][j]/=a[i][i];
      f(j,1,n)
        if(i!=j)
          f(k,i+1,n+1) a[j][k]-=a[j][i]*a[i][k];
  }
  f(i,1,n) printf("%.2lf\n",a[i][n+1]);
  return 0;
}

矩阵快速幂

同正常快速幂

矩阵加速线性递推

例如求斐波那契数列，可以构造矩阵然后矩阵的幂就是答案

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstring>
using namespace std;

struct Matrix {
    long long a[4][4];
}ans;

int mod = 1e9 + 7;;

Matrix Matrix_mul(Matrix A, Matrix B) {
    Matrix C;
    memset(C.a, 0 , sizeof(C.a));
    for(int i = 1; i <= 3; ++i)
        for(int j = 1; j <= 3; ++j)
            for(int k = 1; k <= 3; ++k)
                C.a[i][j] = (C.a[i][j] + A.a[i][k] * B.a[k][j] % mod) % mod;
    return C;
}

Matrix qpow(Matrix A, int n) {
    if(n <= 1) return A;
    Matrix B = A; --n;
    while(n) {
        if(n & 1) B = Matrix_mul(B, A);
        n >>= 1;
        A = Matrix_mul(A, A);
    }
    return B;
}

int main() {
    memset(ans.a, 0, sizeof(ans.a));
    ans.a[1][1] = 1; ans.a[1][2] = 1; ans.a[1][3] = 0;
    ans.a[2][1] = 0; ans.a[2][2] = 0; ans.a[2][3] = 1;
    ans.a[3][1] = 1; ans.a[3][2] = 0; ans.a[3][3] = 0;
    int T, x;
    cin >> T;
    while(T--) {
        cin >> x;
        Matrix A;
        memset(A.a, 0, sizeof(A.a));
        A = qpow(ans, x - 1);
        cout << A.a[1][1] << endl;
    }
    return 0;
}

线性基

xor运算神器

#include<iostream>
#include<algorithm>
#include<cstdio>
#include<vector>
#include<cstring>
#include<queue>
#include<cmath>
#include<map>
#include<set>
using namespace std;
#define REP(i, a, b) for(int i = a; i <= b; i++)
#define PER(i, a, b) for(int i = a; i >= b; i--)
#define LL long long
inline LL read(){
    LL x = 0, flag = 1;char ch = getchar();
    while(!isdigit(ch)) {
        if(ch == '-') flag = - 1;
        ch = getchar();
    }
    while(isdigit(ch)) x = x * 10 + ch - '0', ch = getchar();
    return x * flag;
}

template<typename T> struct linear_base {
    static const int maxlog = 60; 
    T a[61];
    void ins(T x) {
        for(T i = maxlog; i >= 0; --i) {
            if(x & ((T)(1) << i)) {
                if(!a[i]) {
                    a[i] = x;
                    return;
                } else x ^= a[i];
            }
        }
    }

    bool check(T x) {
        for(T i = maxlog; i >= 0; --i) 
            if(x & ((T)(1) << i)) {
                if(!a[i]) return 0;
                else x ^=a[i];
            }
        return 1;
    }

    T qmax() {
        T res = 0;
        for(T i = maxlog; i >= 0; --i)
            if((res ^ a[i]) > res) res ^= a[i];
        return res;
    }

    T qmin() {
        for(T i = 0; i <= maxlog; ++i)
            if(a[i]) return a[i];
    }

    T query(int k) {
        T tmp[61], res = 0, cnt = 0;
        for(T i = 0; i <= maxlog; ++i) {
            for(int j = i - 1; j >= 0; --j)
                if(a[i] & ((T)(1) << j))
                    a[i] ^= a[j];
            if(a[i]) tmp[cnt++] = a[i];
        }
        for(T i = 0; i < cnt; ++i)
            if(k & ((T)(1) << i)) res ^= tmp[i];

        return res;
    }
};

linear_base<LL> s;

int main() {   
    int n = read();
    REP(i, 1, n) s.ins(read());
    printf("%lld\n", s.qmax());
	return 0;
}

经典例题

给出一些运用矩阵运用线性代数的典型例题

UVa10870

Consider recurrent functions of the following form:
f(n) = a1f(n − 1) + a2f(n − 2) + a3f(n − 3) + . . . + adf(n − d), for n > d,
where a1, a2, . . . , ad are arbitrary constants.
A famous example is the Fibonacci sequence, defined as: f(1) = 1, f(2) = 1, f(n) = f(n − 1) +
f(n − 2). Here d = 2, a1 = 1, a2 = 1.
Every such function is completely described by specifying d (which is called the order of recurrence),
values of d coefficients: a1, a2, . . . , ad, and values of f(1), f(2), . . . , f(d). You’ll be given these numbers,
and two integers n and m. Your program’s job is to compute f(n) modulo m.

矩阵加速线性递推，构造相伴矩阵 $F_n = A * F_{n-1}$ A 为对角线为1最后一行为a的矩阵, $F(n)=A^{n-d}*F(d)$

#include<iostream>
#include<string>
#include<cstring>
using namespace std;

const int maxn = 20;
typedef long long Matrix[maxn][maxn];
typedef long long Vector[maxn];

int sz, mod;
void matrix_mul(Matrix A, Matrix B, Matrix res) {
  Matrix C;
  memset(C, 0, sizeof(C));
  for(int i = 0; i < sz; i++)
    for(int j = 0; j < sz; j++)
      for(int k = 0; k < sz; k++) C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % mod;
  memcpy(res, C, sizeof(C));
}

void matrix_pow(Matrix A, int n, Matrix res) {
  Matrix a, r;
  memcpy(a, A, sizeof(a));
  memset(r, 0, sizeof(r));
  for(int i = 0; i < sz; i++) r[i][i] = 1;
  while(n) {
    if(n&1) matrix_mul(r, a, r);
    n >>= 1;
    matrix_mul(a, a, a);
  }
  memcpy(res, r, sizeof(r));
}

void transform(Vector d, Matrix A, Vector res) {
  Vector r;
  memset(r, 0, sizeof(r));
  for(int i = 0; i < sz; i++)
    for(int j = 0; j < sz; j++)
      r[j] = (r[j] + d[i] * A[i][j]) % mod;
  memcpy(res, r, sizeof(r));
}

int main() {
  int d, n, m;
  while(cin >> d >> n >> m && d) {
    Matrix A;
    Vector a, f;

    for(int i = 0; i < d; i++) { cin >> a[i]; a[i] %= m; }
    for(int i = d-1; i >= 0; i--) { cin >> f[i]; f[i] %= m; }

    memset(A, 0, sizeof(A));
    for(int i = 0; i < d; i++) A[i][0] = a[i];
    for(int i = 1; i < d; i++) A[i-1][i] = 1;

    sz = d;
    mod = m;
    matrix_pow(A, n-d, A);
    transform(f, A, f);

    cout << f[0] << endl;
  }
  return 0;
}

UVa1386

有一个细胞自动机，它是一个由 n 个元素组成的环，每个元素的值都必须是 \pmod {m}(modm) 意义下的。现在你需要对这个环进行 k 次操作，每次操作你需要把这个环内每个元素更新成与它距离不超过 d 的所有元素之和（包括自己）。注：每一个新的元素也必须是 \pmod {m}(modm) 意义下的。
首先构造每次操作的矩阵，发现答案是 $v_k=A^k*v_0$ 复杂度 $O(n^3logk)$ 有些高，注意A比较特殊，从第二行开始每一行都是上一行为循环矩阵
这样矩阵乘法的复杂度是 $O(n^2logk)$

#include<cstdio>
#include<vector>
#include<cassert>
using namespace std;

int MOD;

typedef vector<int> CirculantMatrix;

CirculantMatrix operator * (const CirculantMatrix& a, const CirculantMatrix& b) {
  int n = a.size();
  assert(b.size() == n);
  CirculantMatrix c(n);
  for(int i = 0; i < n; i++) { 
    c[i] = 0;
    for(int j = 0; j < n; j++) // ÀÛ¼ÓA(0,j)*B(j,i)
      c[i] = ((long long)a[j]*b[(i-j+n)%n] + c[i]) % MOD; 
  }
  return c;
}

template<typename T>
T fast_pow(const T& a, int k) {
  assert(k > 0);
  if(k == 1) return a;
  T b = fast_pow(a, k/2);
  b = b * b;
  if(k % 2 == 1) b = b * a;
  return b;
}

int main() {
  int n, m, d, k;
  while(scanf("%d%d%d%d", &n, &m, &d, &k) == 4) {
    MOD = m;
    CirculantMatrix mat(n);
    for(int i = 0; i < n; i++)
      mat[i] = min(i, n - i) <= d ? 1 : 0;

    CirculantMatrix v(n);
    for(int i = 0; i < n; i++) scanf("%d", &v[i]);
    v = v * fast_pow(mat, k);

    printf("%d", v[0]);
    for(int i = 1; i < n; i++) printf(" %d", v[i]);
    printf("\n");
  }
  return 0;
}

UVa10828

给出一个程序控制流图，从每个结点出发到每个后继结点的概率相等。当执行完一个没有后继的结点后，整个程序终止。程序总是从编号为1的结点开始执行。你的任务是对于若干个查询结点，求出每个结点的期望执行次数。
这道题是关于马尔可夫过程的。设结点i的出度为di，期望执行次数是xi。对于一个拥有3个前驱结点a，b，c的结点i，可以列出方程xi=xa/da+xb/db+xc/dc。
由于要达到结点i必须先经过结点i的前驱结点a，而弧aài的期望经过次数是xa/da。根据期望的线性性质，上述方程成立。
注意到结点1比较特殊，可以加一个虚拟结点0，并且结点0以概率1转移到结点1。
下面就可以根据题目给出的转移边写出一个方程组并求解。但是需要注意的是输出的答案可能是无穷大，这是因为流图中出现了环，并且可以无限地循环下去。这样高斯消元后得到的方程中存在矛盾方程。
除此之外，高斯消元后还可能得到一些无用方程，从而导致一些变量的值无法确定，这是因为有一些结点无法从起点通过转移到达，这些结点的执行次数期望其实是0。
在解方程的过程中，用数组inf来存储那些无穷变量。

#include<algorithm>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<vector>
using namespace std;

const double eps = 1e-8;
const int maxn = 100 + 10;
typedef double Matrix[maxn][maxn];

void gauss_jordan(Matrix A, int n) {
  int i, j, k, r;
  for(i = 0; i < n; i++) {
    r = i;
    for(j = i+1; j < n; j++)
      if (fabs(A[j][i]) > fabs(A[r][i])) r = j;
    if(fabs(A[r][i]) < eps) continue;
    if(r != i) for(j = 0; j <= n; j++) swap(A[r][j], A[i][j]);
    for(k = 0; k < n; k++) if(k != i)
      for(j = n; j >= i; j--) A[k][j] -= A[k][i]/A[i][i] * A[i][j];
  }
}

Matrix A;
int n, d[maxn];
vector<int> prv[maxn];
int inf[maxn];

int main() {
  int kase = 0;
  while(scanf("%d", &n) == 1 && n) {
    memset(d, 0, sizeof(d));
    for(int i = 0; i < n; i++) prv[i].clear();

    int a, b;
    while(scanf("%d%d", &a, &b) == 2 && a) {
      a--; b--; 
      d[a]++; 
      prv[b].push_back(a);
    }
    memset(A, 0, sizeof(A));
    for(int i = 0; i < n; i++) {
      A[i][i] = 1;
      for(int j = 0; j < prv[i].size(); j++)
        A[i][prv[i][j]] -= 1.0 / d[prv[i][j]];
      if(i == 0) A[i][n] = 1;
    }
    gauss_jordan(A, n);
    memset(inf, 0, sizeof(inf));
    for(int i = n-1; i >= 0; i--) {
      if(fabs(A[i][i])<eps && fabs(A[i][n])>eps) inf[i] = 1; 
      for(int j = i+1; j < n; j++)
        if(fabs(A[i][j])>eps && inf[j]) inf[i] = 1; 
    }

    int q, u;
    scanf("%d", &q);
    printf("Case #%d:\n", ++kase);
    while(q--) {
      scanf("%d", &u); u--;
      if(inf[u]) printf("infinity\n");
      else printf("%.3lf\n", fabs(A[u][u])<eps ? 0.0 : A[u][n]/A[u][u]);
    }
  }
  return 0;
}

UVa11542

给出n个整数，从中选出1个或多个，使得选出整数的乘积是完全平方数。一共有多少种选法？（每个数的素因子小于500）
不大于500的素因子可以求出每个数的唯一分解式用（0，1）表示这个数取或不取，对于每一个因子要求奇偶情况下xor起来=0可以高斯消元答案是自由变量个数的2次方

#include<algorithm>
#include<cmath>
#include<iostream>
#include<cstdio>
#include<cstring>
#include<vector>
using namespace std;

const int maxn = 500 + 10;
const int maxp = 100;

int vis[maxn];
int prime[maxp];

int gen_primes(int n) {
  int m = (int)sqrt(n+0.5);
  memset(vis, 0, sizeof(vis));
  for(int i = 2; i <= m; i++) if(!vis[i])
    for(int j = i*i; j <= n; j+=i) vis[j] = 1;
  int c = 0;
  for(int i = 2; i <= n; i++) if(!vis[i])
    prime[c++] = i;
  return c;
}

typedef int Matrix[maxn][maxn];

int _rank(Matrix A, int m, int n) {
  int i = 0, j = 0, k, r, u;
  while(i < m && j < n) { 
    r = i;
    for(k = i; k < m; k++)
      if(A[k][j]) { r = k; break; }
    if(A[r][j]) {
      if(r != i) for(k = 0; k <= n; k++) swap(A[r][k], A[i][k]);
      for(u = i+1; u < m; u++) if(A[u][j])
        for(k = i; k <= n; k++) A[u][k] ^= A[i][k];
      i++;
    }
    j++;
  }
  return i;
}

Matrix A;

int main() {
  int m = gen_primes(500);

  int T;
  cin >> T;
  while(T--) {
    int n, maxp = 0;
    long long x;
    cin >> n;
    memset(A, 0, sizeof(A));
    for(int i = 0; i < n; i++) {
      cin >> x;
      for(int j = 0; j < m; j++) 
        while(x % prime[j] == 0) {
          maxp = max(maxp, j); x /= prime[j]; A[j][i] ^= 1;
        }
    }
    int r = _rank(A, maxp+1, n);
    cout << (1LL << (n-r))-1 << endl; 
  }
  return 0;
}

经典题目1

给定n个点，m个操作，构造O(m+n)的算法输出m个操作后各点的位置。操作有平移、缩放、翻转和旋转
这里的操作是对所有点同时进行的。其中翻转是以坐标轴为对称轴进行翻转（两种情况），旋转则以原点为中心。如果对每个点分别进行模拟，那么m个操作总共耗时O(mn)。利用矩阵乘法可以在O(m)的时间里把所有操作合并为一个矩阵，然后每个点与该矩阵相乘即可直接得出最终该点的位置，总共耗时O(m+n)。假设初始时某个点的坐标为x和y，下面5个矩阵可以分别对其进行平移、旋转、翻转和旋转操作。预先把所有m个操作所对应的矩阵全部乘起来，再乘以(x,y,1)，即可一步得出最终点的位置。
在这里插入图片描述

经典题目2

给定矩阵A，请快速计算出A^n（n个A相乘）的结果，输出的每个数都mod p。
由于矩阵乘法具有结合律，因此 $A^4 = A * A * A * A = (A*A) * (A*A) = A^2 * A^2$ 。我们可以得到这样的结论：当n为偶数时， $A^n = A^{(n/2)} * A^{(n/2)}；$ 当n为奇数时， $A^n = A^{(n/2)} * A^{(n/2)} * A$ （其中n/2取整）。这就告诉我们，计算 $A^n$ 也可以使用二分快速求幂的方法。例如，为了算出 $A^25$ 的值，我们只需要递归地计算出 $A^12$ 、 $A^6$ 、 $A^3$ 的值即可。根据这里的一些结果，我们可以在计算过程中不断取模，避免高精度运算。

经典题目3

题目大意：给定矩阵A，求 $A + A^2 + A^3 + … + A^k$ 的结果（两个矩阵相加就是对应位置分别相加）。输出的数据 $mod m。k<=10^9$ 。
这道题两次二分，相当经典。首先我们知道， $A^i$ 可以二分求出。然后我们需要对整个题目的数据规模k进行二分。比如，当k=6时，有：
$A + A^2 + A^3 + A^4 + A^5 + A^6 =(A + A^2 + A^3) + A^3*(A + A^2 + A^3)$
应用这个式子后，规模k减小了一半。我们二分求出 $A^3$ 后再递归地计算 $A + A^2 + A^3$ ，即可得到原问题的答案。

经典题目4

给定一个有向图，问从A点恰好走k步（允许重复经过边）到达B点的方案数mod p的值
把给定的图转为邻接矩阵，即 $A (i, j) = 1$ 当且仅当存在一条边i->j。令 $C = A * A$ ，那么 $C (i, j) = Σ A (i, k) * A (k, j)$ ，实际上就等于从点i到点j恰好经过2条边的路径数（枚举k为中转点）。类似地， $C * A$ 的第i行第j列就表示从i到j经过3条边的路径数。同理，如果要求经过k步的路径数，我们只需要二分求出 $A^k$ 即可。

经典题目5

用1 x 2的多米诺骨牌填满M x N的矩形有多少种方案，M<=5，N<2^31，输出答案mod p的结果
状态压缩先搜出转移再矩乘

练习题

bzoj1013 bzoj4386 bzoj2476 bzoj1770 bzoj2115 bzoj1923 bzoj3105 bzoj2460 bzoj3270 bzoj1951
luogu上高斯消元、线性基、矩阵乘法、矩阵加速通过超过100人的省选及经典题目