基础线性代数大记（二）三道高消题

目录

• 笔者的睿智前言
• 矩阵逆
• Broken Robot
• HNOI2011 XOR和路径

---

笔者的睿智前言

[我 ext{ TMD } 不想搬砖 ]

所以……

我要搬砖

这傻吊游戏， 不玩也罢

---

矩阵逆

将矩阵 (A_{nn}) 消元成上三角矩阵后， 如果主元的个数为 n， 那么该矩阵可逆。主元不为 n 的 n 阶向量组是线性相关的。

矩阵的初等变换可以通过矩阵乘法实现， 所以通过初等变换将一个矩阵变成单位矩阵， 等于将这个矩阵乘上了这个矩阵的逆矩阵， 记录下初等变换的所有细节， 对一个单位矩阵做同一操作， 就可以得到原来矩阵的逆矩阵了。

至于矩阵是否可逆， 可以用高斯消元判。

写的模板是洛谷的， 我天生自带大常数死得很惨， 开 O2 优化才能过。

#include<bits/stdc++.h>
typedef long long lo;
using namespace std;

const int N = 405;
const int mod = 1000000007;
int n;

lo poW(lo a, int b) {
    lo res = 1ll;
    for(;b;b>>=1,a=(a*a)%mod)if(b&1)res=(res*a)%mod;
    return res%mod;
}

struct mat{
    lo a[N][N];
    mat(){memset(a,0,sizeof a);}
    void Swap(int i,int j) {
        for(int k=1;k<=n;++k) swap(a[i][k], a[j][k]);
    }
    void Mul(int x, lo k) {
        for(int i=1;i<=n;++i) {
            a[x][i] *= k;
            a[x][i] %= mod;
        }
    }
    void Ins(int x,int y,lo k) {
        for(int i=1;i<=n;++i) {
            a[x][i] += a[y][i]*k%mod;
            a[x][i] %= mod;
        }
    }
    void print() {
        for(int i=1;i<=n;++i,putchar('
'))
        for(int j=1;j<=n;++j) cout << (a[i][j]+mod)%mod<< ' ';
        putchar('
');
    }
} A,B;



signed main() {
    scanf("%d", &n);
    for(int i=1;i<=n;++i) {
        for(int j=1;j<=n;++j) {
            scanf("%lld", &A.a[i][j]);
            A.a[i][j] %= mod;
        }
    }
    for(int i=1;i<=n;++i) B.a[i][i]=1;
    
    for(int i=1;i<=n;++i) {
        if(A.a[i][i]%mod==0)
            for(int j=i;j<=n;++j) if(A.a[j][i]%mod) {
                B.Swap(i,j);
                A.Swap(i,j);
                break;
            }
        if(A.a[i][i]%mod==0) {
            puts("No Solution");return 0;
        }
        A.a[i][i] = (A.a[i][i]%mod+mod)%mod;
        B.Mul(i,poW(A.a[i][i],mod-2));
        A.Mul(i,poW(A.a[i][i],mod-2));
        for(int j=i+1;j<=n;++j) {
            B.Ins(j, i, -A.a[j][i]);
            A.Ins(j, i, -A.a[j][i]);
        }
    }
    for(int i=n;i>=1;--i) {
        for(int j=i-1;j>=1;--j) {
            B.Ins(j,i,-A.a[j][i]);
            A.Ins(j,i,-A.a[j][i]);
        }
    }
    B.print();
    return 0;
}

---

Broken Robot

列方程， 高消即可。

注意特判 m=1 的情况， 这种数据不能被列方程解决。

#include<bits/stdc++.h>
using namespace std;
const int N = 1003;

int n,m,x,y;
double a[N][N];

void Gauss()
{
    for(int i=1;i<=m;++i) {
        //i,i i,i+1
        double tmp = a[i][i];
        a[i][i+1]/=tmp, a[i][i]/=tmp;
        if(i==m) break;
        a[i][m+1]/=tmp;
        tmp = a[i+1][i];
        a[i+1][i]-=tmp*a[i][i], a[i+1][i+1]-=tmp*a[i][i+1], a[i+1][m+1]-=tmp*a[i][m+1];
    }
    for(int i=m;i>1;--i) {
        a[i-1][m+1] -= a[i-1][i]*a[i][m+1];
    }
}

int main()
{
    scanf("%d%d%d%d", &n,&m,&x,&y);
    if(m==1) {
        printf("%.4lf", (double)2 * (n - 1));
        return 0;
    }
    n -= (x-1);
    // 2/3 f[i,1] - 1/3 f[i,2] =  1/3 f[i+1,1] + 1
    //-1/4 f[i,j-1] + 3/4 f[i,j] - 1/4 f[i,j+1] = 1/4 f[i+1,j] + 1
    // f[i,m] = 1/3 f[i,m] + 1/3 f[i,m-1] + 1/3 f[i+1,m] + 1
    for(int T=1;T<n;++T) {
        a[1][m+1]=a[1][m+1]/3.0+1, a[m][m+1]=a[m][m+1]/3.0+1;
        for(int j=2;j<m;++j) a[j][m+1] = a[j][m+1]/4.0 + 1;
        a[1][1]=2.0/3.0, a[1][2]=-1.0/3.0;
        a[m][m]=2.0/3.0, a[m][m-1]=-1.0/3.0;
        for(int j=2;j<m;++j) {
            a[j][j] = 3.0/4.0;
            a[j][j-1] = a[j][j+1] = -1.0/4.0;
        }
        Gauss();
    }
    printf("%.4lf", a[y][m+1]);
    return 0;
}

但上面那份代码依然是有漏洞的， 把标程放在这里：

//Author:XuHt
#include <cstdio>
#include <iostream>
using namespace std;
const int N = 1006;
int n, m, x, y;
double f[N][N], d[N][N];

void work() {
	for (int i = 1; i <= m; i++) {
		double w = 1.0 / d[i][i];
		d[i][i] *= w;
		d[i][m+1] *= w;
		if (i == m) break;
		d[i][i+1] *= w;
		w = d[i+1][i] / d[i][i];
		d[i+1][i] -= w * d[i][i];
		d[i+1][i+1] -= w * d[i][i+1];
		d[i+1][m+1] -= w * d[i][m+1];
	}
	for (int i = m - 1; i; i--)
		d[i][m+1] -= d[i+1][m+1] * d[i][i+1];
}

int main() {
	cin >> n >> m >> x >> y;
	for (int i = n - 1; i >= x; i--) {
		d[1][1] = d[m][m] = -2 / 3.0;
		d[1][2] = d[m][m-1] = 1 / 3.0;
		for (int j = 2; j < m; j++) {
			d[j][m+1] = -f[i+1][j] / 4.0 - 1;
			d[j][j] = -3 / 4.0;
			d[j][j-1] = d[j][j+1] = 1 / 4.0;
		}
		if (m == 1) d[1][1] = -1 / 2.0;
		d[1][m+1] = -f[i+1][1] / 3.0 - 1;
		d[m][m+1] = -f[i+1][m] / 3.0 - 1;
		if (m == 1) d[m][m+1] = -f[i+1][m] / 2.0 - 1;
		work();
		for (int j = 1; j <= m; j++) f[i][j] = d[j][m+1];
	}
	printf("%.10f
", f[x][y]);
	return 0;
}

---

HNOI2011 XOR和路径

根据期望的线性性， 求出每一位的期望然后相加。

(f(i,x)) 表示第 (i) 位， x 到 n 路径 Xor 值不为 0 的概率。

[f(i,x) = frac{1}{k(x)}sumlimits_{smid w(x,s){ m xor}2^i=1}(1-f(i,s)) + frac{1}{k(x)}sumlimits_{smid w(x,s){ m xor}2^i=0}f(i,s) ]

#include<bits/stdc++.h>
using namespace std;
const int N = 103;
const int M = 10003;

int n,m;
int ct, hd[N], nt[M<<1], vr[M<<1], w[M<<1];

void ad(int x,int y,int z) {
    vr[++ct]=y, w[ct]=z, nt[ct]=hd[x], hd[x]=ct;
}

struct mat{
    double a[N][N];
    mat(){memset(a,0,sizeof a);}
    void Swap(int x,int y) {
        for(int i=1;i<=n+1;++i) swap(a[x][i],a[y][i]);
    }
    void Mul(int x,double k) {
        for(int i=1;i<=n+1;++i) a[x][i] *= k;
    }
    void Ins(int x,int y,double k) {
        for(int i=1;i<=n+1;++i) a[x][i] += a[y][i]*k;
    }
};

double sol(int t) {
    mat A;
    for(int x=1;x<=n;++x) {
        A.a[x][x] = 1;
        if(x==n) continue;
        int sum=0;
        for(int i=hd[x];i;i=nt[i]) {
            ++sum;
        }
        for(int i=hd[x];i;i=nt[i]) {
            int y = vr[i];
            if((w[i]>>t)&1) {
                A.a[x][n+1] += 1.0/sum;
                A.a[x][y] += 1.0/sum;
            } else {
                A.a[x][y] += -1.0/sum;
            }
        }
    }
    
    for(int i=1;i<=n;++i) {
        int mx=i; for(int j=i+1;j<=n;++j)if(fabs(A.a[j][i])>fabs(A.a[mx][i]))mx=j;
        if(mx^i) A.Swap(i,mx);
        A.Mul(i,1/A.a[i][i]);
        for(int j=i+1;j<=n;++j) {
            A.Ins(j,i,-A.a[j][i]);
        }
    }
    
    for(int i=n;i>=1;--i) {
        double tmp = A.a[i][i];
        A.a[i][i]/=tmp, A.a[i][n+1]/=tmp;
        for(int j=i-1;j>=1;--j) A.Ins(j,i,-A.a[j][i]);
    }
    return A.a[1][n+1];
}

int main()
{
    cin>>n>>m;
    for(int i=0;i<m;++i) {
        int u,v,w; scanf("%d%d%d",&u,&v,&w);
        ad(u,v,w); if(u^v)ad(v,u,w);
    }
    double Ans = 0;
    for(int T=0;T<=30;++T) Ans += sol(T)*(1<<T);
    printf("%.3lf", Ans);
    return 0;
}