hihoCoder1195 1196(高斯消元,异或方程组)

高斯消元     异或方程组

1.高斯消元

我们把每一件商品的价格看作是x[1]..x[n],第i个组合中第j件商品数量记为a[i][j],其价格记作y[i],则可以列出方程式:

a[1][1] * x[1] + a[1][2] * x[2] + ... + a[1][n] * x[n] = y[1]
a[2][1] * x[1] + a[2][2] * x[2] + ... + a[2][n] * x[n] = y[2]
                                  ...
a[m][1] * x[1] + a[m][2] * x[2] + ... + a[m][n] * x[n] = y[m]

我们可以对方程组进行3种操作而不改变方程组的解集:

1. 交换两行。

2. 把第i行乘以一个非0系数k。即对于j = 1..n, 令a[i][j] = k*a[i][j], y[i]=k*y[i]

3. 把第p行乘以一个非0系数k之后加在第i行上。即对于j=1..n, 令a[i][j] = a[i][j]+k*a[p][j],y[i]=y[i]+k*p[i]

以上三个操作叫做初等行变换。

我们可以使用它们,对这个方程组中的a[i][j]进行加减乘除变换,举个例子:

a[1][1] * x[1] + a[1][2] * x[2] + ... + a[1][n] * x[n] = y[1]    式子(1)
a[2][1] * x[1] + a[2][2] * x[2] + ... + a[2][n] * x[n] = y[2]    式子(2)

我们可以通过 式子(1) - 式子(2) * (a[1][1] / a[2][1]),将第1行第1列的a[1][1]变换为0。

对整个方程组进行多次变换之后,可以使得a[i][j]满足:

a[i][j] = 1 (i == j)
a[i][j] = 0 (i != j)

则整个方程组变成了:

x[1] = y'[1]
x[2] = y'[2]
...
x[n] = y'[n]
0 = y'[n + 1]
0 = y'[n + 2]
...
0 = y'[m]

这样的话,y'[1] .. y'[n]就是我们要求的x[1]..x[n]

小Hi:挺不错的嘛,继续?

小Ho:好,关于如何变换,我们可以利用一个叫高斯消元的算法。高斯消元分成了2个步骤:

首先我们要计算出上三角矩阵,也就是将方程组变为:

a[1][1] * x[1] + a[1][2] * x[2] + ... + a[1][n] * x[n] = y'[1]
      0 * x[1] + a[2][2] * x[2] + ... + a[2][n] * x[n] = y'[2]
      0 * x[1] +       0 * x[2] + ... + a[3][n] * x[n] = y'[3]
                                   ...
      0 * x[1] +       0 * x[2] + ... + a[n][n] * x[n] = y'[n]
      0 * x[1] +       0 * x[2] + ... +       0 * x[n] = y'[n + 1]
	                               ...
      0 * x[1] +       0 * x[2] + ... +       0 * x[n] = y'[m]

也就是通过变换,将所有a[i][j](i>j)变换为0。同时要保证对角线上的元素a[i][i]不为0。

方法也很见简单,从第1行开始,我们利用当前行第i列不为0,就可以通过变换将i+1..M行第一列全部变换为0,接着对于第2行,我们用同样的方法将第3..M行第2列也变换为0...不断重复直到第n行为止。

假如计算到第i行时,第i列已经为0,则我们需要在第i+1..M行中找到一行第i列不为0的行k,并交换第i行和第k行,来保证a[i][i] != 0。但这时候还有可能出现一个情况,就是第i..M行中的i列均为0,此时可以判定,该方程组有多解。


当得到上三角矩阵后,就可以从第n行开始逆推,一步一步将a[i][j](i<j)也变换为0.

因为第n行为a[n][n] * x[n] = y'[n],则x[n] = y'[n] / a[n][n]。

第n-1行为a[n-1][n-1] * x[n - 1] + a[n][n] * x[n] = y'[n - 1]。我们将得到的x[n]代入,即可计算出x[n-1]。

同样的依次类推就可以得到所有的x[1]..x[n]。


而对于多解和无解的判定:

当在求出的上三角矩阵中出现了 a[i][1] = a[i][2] = ... = a[i][n] = 0, 但是y'[i] != 0时,产生了矛盾,即出现了无解的情况。

而多解的证明如下:

假设n=3,m=3,而我们计算出了上三角矩阵为:

a * x[1] + b * x[2] + c * x[3] = d
                      e * x[3] = f
                             0 = 0

当我们在第一个式子中消去x[3]后,有a * x[1] + b * x[2] = g,显然x[1]和x[2]有无穷多种可能的取值。

小Hi:既然小Ho你都已经把整个算法讲了,那么我就只能给出伪代码了:

// 处理出上三角矩阵
For i = 1..N
    Flag ← False
    For j = i..M                // 从第i行开始,找到第i列不等于0的行j
        If a[j][i] != 0 Then
            Swap(j, i)          // 交换第i行和第j行
            Flag ← True
            Break
        End If
    End For
    // 若无法找到,则存在多个解
    If (not Flag) Then
        manySolutionsFlag ← True // 出现了秩小于N的情况
        continue;
    End If
    // 消除第i+1行到第M行的第i列
    For j = i+1 .. M
        a[j][] ← a[j][] - a[i][] * (a[j][i] / a[i][i])
        b[j] ← b[j] - b[i] * (a[j][i] / a[i][i])
    End For
End For 

// 检查是否无解,即存在 0 = x 的情况
For i = 1..M
    If (第i行系数均为0 and (b[i] != 0)) Then
        return "No solutions"
    End If
End For

// 判定多解
If (manySolutionsFlag) Then
	return "Many solutions"
End If

// 此时存在唯一解
// 由于每一行都比前一行少一个系数,所以在M行中只有前N行有系数
// 解析来从第N行开始处理每一行的解
For i = N .. 1
    // 利用已经计算出的结果,将第i行中第i+1列至第N列的系数消除
    For j = i + 1 .. N
        b[i] ← b[i] - a[i][j] * value[j]
        a[i][j] ← 0
    End For
    value[i] ← b[i] / a[i][i]
End For

AC代码:

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>

#define fbs 1e-6

using namespace std;

double a[1001][501];
double c[1001];
double x[501];

int N,M;

void Swap(int i,int j)   //交换i,j行
{
    for(int k=1;k<=N;k++)
        swap(a[i][k],a[j][k]);
    swap(c[i],c[j]);
}

int judge(int i)   //判断是否全为0
{
    for(int j=1;j<=N;j++)
        if(a[i][j]!=0)
            return 0;
    return 1;
}

int Gauss()
{
    for(int i=1;i<=N;i++)  //处理出上三角矩阵
    {
        int flag=0;
        for(int j=i;j<=M;j++)
        {
            if(a[j][i]!=0)
            {
                Swap(j,i);
                flag=1;
                break;
            }
        }
        if(!flag)   //若无法找到,则存在多解
        {
            printf("Many solutions\n");
            return 0;
        }

        for(int j=i+1;j<=M;j++)   //消除第i+1行到第M行的第i列
        {
            double mul=a[j][i]/a[i][i];
            for(int k=i;k<=N;k++)
                a[j][k]-=mul*a[i][k];
            c[j]-=c[i]*mul;
        }
    }

    for(int i=1;i<=M;i++)  //判断是否无解
        if(judge(i)&&(c[i]>fbs||c[i]<-fbs))
        {
            printf("No solutions\n");
            return 0;
        }

    //求出唯一解
    for(int i=N;i>0;i--)
    {
        for(int j=i+1;j<=N;j++)
        {
            c[i]-=a[i][j]*x[j];
            a[i][j]=0;
        }
        x[i]=c[i]/a[i][i];
    }
    return 1;
}

int main()
{
    scanf("%d %d",&N,&M);
    for(int i=1;i<=M;i++)
    {
        for(int j=1;j<=N;j++)
            scanf("%lf",&a[i][j]);
        scanf("%lf",&c[i]);
    }
    if(Gauss())
        for(int i=1;i<=N;i++)
            printf("%d\n",int(x[i]+0.5));

    return 0;
}


2.异或方程组

首先对于每一个格子的状态,可能会对它造成影响的是其自身和周围4个格子,这五个格子被按下的总次数也就等于该格子所改变的总次数。

对于任意一个格子,如果这个格子改变了偶数次状态,则等价于没有发生改变。

我们可以将1看作格子亮着,0看作格子暗着,每改变1次就加1,最后格子的状态等于其总数值 MOD 2。

则其运算结果刚好满足异或运算,即每改变一次等于状态值 xor 1。

同样的对于一个格子和它周围的4个格子来说,若格子被按下偶数次,它自身和周围4个格子的状态也等于没有发生改变。所以我们可以知道:任意一个格子至多被按下一次。

假设有数组x[1..30],分别表示这30个格子是否按下1次,若按下则x[i]=1,否则x[i]=0。

则对于1个格子,他最后的状态为:

当前状态 = 初始状态 xor (a[1] * x[1]) xor (a[2] * x[2]) xor ... xor (a[30] * x[30])

其中a[i]表示格子i是否会对当前格子产生影响,若能够则a[i] = 1,否则a[i] = 0

对方程进行变换有:

(a[1] * x[1]) xor (a[2] * x[2]) xor ... xor (a[30] * x[30]) = 当前状态 xor 初始状态

因为我们的目标是要让所有等格子都为亮的状态,故我们需要让 当前状态 = 1,则:

(a[1] * x[1]) xor (a[2] * x[2]) xor ... xor (a[30] * x[30]) = 1 xor 初始状态

不妨设y = 1 xor 初始状态:

(a[1] * x[1]) xor (a[2] * x[2]) xor ... xor (a[30] * x[30]) = y

对于所有的格子,我们可以连立出方程组:

(a[ 1][1] * x[1]) xor (a[ 1][2] * x[2]) xor ... xor (a[ 1][30] * x[30]) = y[ 1]
(a[ 2][1] * x[1]) xor (a[ 2][2] * x[2]) xor ... xor (a[ 2][30] * x[30]) = y[ 2]
                                            ...
(a[30][1] * x[1]) xor (a[30][2] * x[2]) xor ... xor (a[30][30] * x[30]) = y[30]
		

到此,我们的目标就是求出一个x[1..30],使得上面的方程组成立。

小Ho:这个看上去和高斯消元很像啊。

小Hi:没错,这个方程组叫异或方程组,它可以用和高斯消元同样的方法来解决。

其解答过程几乎和高斯消元无异,判定无解和多解的方式也相同。唯一需要注意的是消元过程不再是高斯消元的加减,而是通过xor运算来进行消元。比如消除第j行第i列的1:

a[j][k] = a[j][k] xor a[i][k], y[j] = y[j] xor y[i]

其原理是:

    (a[j][1] * x[1]) xor (a[j][2] * x[2]) xor ... xor (a[j][30] * x[30]) xor (a[i][1] * x[1]) xor (a[i][2] * x[2]) xor ... xor (a[ i][30] * x[30]) = y[j] xor y[i]
<=> ((a[j][1] * x[1]) xor (a[i][1] * x[1])) xor (((a[j][2] * x[2]) xor (a[i][2] * x[2]))) xor ... xor ((a[j][30] * x[30]) xor (a[i][30] * x[30])) = y[j] xor y[i]<=> ((a[j][1] xor a[i][1]) * x[1]) xor ((a[j][2] xor a[i][2]) * x[2]) xor ... ((a[j][30] xor a[i][30]) * x[30]) = y[j] xor y[i]

而且由于给定游戏板是固定的,我们可以知道a[i][j]矩阵一定是固定的,而且通过计算可以知道我们消元得到的上三角矩阵也是固定的,并且在这一次的问题中该上三角矩阵是满秩的,所以其一定存在唯一解。


AC代码:

#include<iostream>
#include<cstdio>
#include<cstring>
#include<vector>
#include<algorithm>

using namespace std;

int c[31];
int a[31][31];   //i受j的影响
int go[4][2]={{1,0},{-1,0},{0,1},{0,-1}};
int x[31];
vector<int>v;

void Init()
{
    int x,y,now,then;
    memset(a,0,sizeof(a));
    v.clear();
    for(int i=1;i<=5;i++)
    {
        for(int j=1;j<=6;j++)
        {
            now=6*(i-1)+j;
            a[now][now]=1;
            for(int k=0;k<4;k++)
            {
                x=i+go[k][0];
                y=j+go[k][1];
                if(x>=1&&x<=5&&y>=1&&y<=6)
                {
                    then=6*(x-1)+y;
                    a[now][then]=1;
                }
            }
        }
    }
}

void Swap(int i,int j)   //交换i,j行
{
    for(int k=1;k<=30;k++)
        swap(a[i][k],a[j][k]);
    swap(c[i],c[j]);
}

void Gauss()
{
    for(int i=1;i<=30;i++)  //处理出上三角矩阵
    {
        for(int j=i;j<=30;j++)
        {
            if(a[j][i])
            {
                Swap(j,i);
                break;
            }
        }

        for(int j=1;j<=30;j++)   //消除第i+1行到第M行的第i列
        {
            if(j!=i&&a[j][i])
            {
                for(int k=i;k<=30;k++)
                    a[j][k]^=a[i][k];
                c[j]^=c[i];
            }
        }
    }
}


int main()
{
    char s[7];
    int k=1;
    for(int i=0;i<5;i++)
    {
        scanf("%s",s);
        for(int j=0;j<6;j++)
            c[k++]=(s[j]-'0')^1;
    }
    Init();
    Gauss();
    int sum=0;
    for(int i=1;i<=30;i++)
    {
        if(c[i])
        {
            v.push_back(i);
            sum++;
        }
    }
    printf("%d\n",sum);
    int x,y;
    sort(v.begin(),v.end());
    for(int i=0;i<v.size();i++)
    {
        x=(v[i]-1)/6+1;
        y=v[i]%6;
        if(y==0)
            y=6;
        printf("%d %d\n",x,y);
    }

    return 0;
}

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值