[学习笔记]高斯消元求解两种特殊问题(带状矩阵/主元法)

本文章是[学习笔记]概率与期望进阶的一部分
由于时间问题我写的比较简略,所以我把大佬的总结链接贴上来了(应该没什么吧qwq)。

1 概述

最常见的当然是随机游走问题了…

f u = ∑ p u , v ∗ ( f v + w u , v ) f_u=\sum p_{u,v}* (f_{v}+w_{u,v}) fu=pu,v(fv+wu,v)
• 计算期望在这个节点上,停留多少步: f u = ∑ p v , u ∗ f v + [ u = S ] f_u=\sum p_{v,u}* f_v + [u=S] fu=pv,ufv+[u=S]
然后是两道水题:
HNOI 2013 游走
SDOI 2012 走迷宫

2 Band Matrix

LSK大佬的总结博客
IDnumber4大佬的总结博客
Lucky_Glass大佬的总结博客(好难找…)
cqbzlytina大佬的例题题解

带宽:对于 ∀ ( i , j ) ∈ G , a i , j ≠ 0 存 在 一 个 x , ( x , x ) ∈ G 且 ∣ i − x ∣ + ∣ j − x ∣ ≤ d \forall (i,j)\in G,a_{i,j}≠0 存在一个x,(x,x)\in G且|i-x|+|j-x|\leq d (i,j)G,ai,j=0x,(x,x)Gix+jxd,当 d d d取到最小时就为Band Matrix的带宽。

在做高斯消元时,利用带宽的性质,每次消元只需要处理下 d d d行,下 d d d列即可。
复杂度就成了 O ( n d 2 ) O(nd^2) O(nd2)
注意由于Band Matrix 性质,我们这里就没必要去交换行了。

Berlekamp-Massey算法(并没有什么关系)

CF 963 E Circles of Waiting

蒟蒻的部分实现:

int Gauss(int m,int n){//m行n列
    int r=1,c=1;
    for(;r<=m&&c<=n;r++,c++){
        if(!fun[r][c])return -1;
        int v=inv(fun[r][c]);
        for(int i=c;i<=min(c+D,n);i++)
        fun[r][i]=mul(fun[r][i],v);
        fun[r][n+1]=mul(fun[r][n+1],v);
        for(int i=r+1;i<=min(r+D,m);i++)
        if(fun[i][c]){
            v=fun[i][c];
            for(int j=c;j<=min(c+D,n);j++)
            fun[i][j]=dec(fun[i][j],mul(fun[r][j],v));
            fun[i][n+1]=dec(fun[i][n+1],mul(fun[r][n+1],v));
        }
    }
    for(int i=r;i<=m;i++)
    if(fun[i][n+1])return -1;
    if(c<n)return n-c;
    for(int i=n;i>=1;i--){
        int v=fun[i][n+1];
        for(int j=i+1;j<=min(n,i+D);j++)
        v=dec(v,mul(ans[j],fun[i][j]));
        ans[i]=v;
    }
    return 0;
}

3 主元法(代入消元法)

ixRic大佬的例题题解

主要用于网格图,我们只需要保留最后一行和最后一列的元素,会发现其他的元素都可以被最后一行最后一列的元素表示出来,这样我们就只需要用最后一行,最后一列的元素来做高斯消元了。

TopCoder 12984 TorusSailing

struct node{
    DB r[MAXN+5],v;//该点表示的方程
}d[MAXN+5][MAXN+5];

DB fun[MAXN+5][MAXN+5],ans[MAXN+5];
DB pabs(DB x){return (x<0)?-x:x;}
int Gauss(int m,int n){//m行n列
    int col=0,k=0;
    for(;k<m&&col<n;k++,col++){
        int r=k;
        for(int i=k+1;i<m;i++)
        if(pabs(fun[r][col])<pabs(fun[i][col]))r=i;
        if(pabs(fun[r][col])<eps){k--;continue;}
        if(r!=k){
            for(int i=col;i<=n;i++)
            swap(fun[k][i],fun[r][i]);
        }
        for(int i=k+1;i<m;i++)
        if(pabs(fun[i][col])>eps){
            DB Div=fun[i][col]/fun[k][col];
            for(int j=col;j<=n;j++)
            fun[i][j]-=fun[k][j]*Div;
            fun[i][col]=0;
        }
    }
    for(int i=k;i<m;i++)
    if(pabs(fun[i][n])>eps)return -1;
    if(k<n)return n-k;
    for(int i=n-1;i>=0;i--){
        DB tmp=fun[i][n];
        for(int j=i+1;j<n;j++)
        tmp-=ans[j]*fun[i][j];
        ans[i]=tmp/fun[i][i];
    }
    return 0;
}

class TorusSailing{
	public:
    DB expectedTime(int n,int m,int sx,int sy){
        for(int i=0;i<n-1;i++)d[i][m-1].r[i]=1;
        for(int i=0;i<m-1;i++)d[n-1][i].r[i+n-1]=1;
        for(int i=n-2;i>=0;i--)
            for(int j=m-2;j>=0;j--){
                d[i][j].v=0.5*(d[i+1][j].v+d[i][j+1].v)+1;
                for(int k=0;k<n+m-2;k++)
                d[i][j].r[k]=0.5*(d[i+1][j].r[k]+d[i][j+1].r[k]);
            }
        for(int i=0;i<n-1;i++){
            fun[i][i]=1,fun[i][n+m-2]=0.5*d[i][0].v+1;
            for(int k=0;k<n+m-2;k++)
            fun[i][k]-=0.5*(d[i][0].r[k]+d[i+1][m-1].r[k]);
        }
        for(int i=0;i<m-1;i++){
            int j=i+n-1;
            fun[j][j]=1,fun[j][n+m-2]=0.5*d[0][i].v+1;
            for(int k=0;k<n+m-2;k++)
            fun[j][k]-=0.5*(d[0][i].r[k]+d[n-1][i+1].r[k]);
        }
        Gauss(n+m-2,n+m-2);
        sx=n-1-sx,sy=m-1-sy;
        DB res=d[sx][sy].v;
        for(int k=0;k<n+m-2;k++)
        res+=ans[k]*d[sx][sy].r[k];
        return res;
    }
};
  • 0
    点赞
  • 6
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值