`
yzmduncan
  • 浏览: 326416 次
  • 性别: Icon_minigender_1
  • 来自: 武汉
文章分类
社区版块
存档分类
最新评论

解线性方程组——高斯消元法

阅读更多



 

例:ZOJ3645

题意:高斯消元模板题(浮点型)

 

/**
高斯消元求解线性方程组.
*/
#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;

///高斯消元模板
const double eps = 1e-12;
const int Max_M = 15;       ///m个方程,n个变量
const int Max_N = 15;
int m, n;
double Aug[Max_M][Max_N+1]; ///增广矩阵
bool free_x[Max_N];         ///判断是否是不确定的变元
double x[Max_N];            ///解集

int sign(double x){ return (x>eps) - (x<-eps);}

/**
返回值:
-1 无解
0 有且仅有一个解
>=1 有多个解,根据free_x判断哪些是不确定的解
*/
int Gauss()
{
    int i,j;
    int row,col,max_r;
    for(row=0,col=0; row<m&&col<n; row++,col++)
    {
        max_r = row;
        for(i = row+1; i < m; i++)  ///找到当前列所有行中的最大值(做除法时减小误差)
        {
            if(sign(fabs(Aug[i][col])-fabs(Aug[max_r][col]))>0)
                max_r = i;
        }
        if(max_r != row)            ///将该行与当前行交换
        {
            for(j = row; j < n+1; j++)
                swap(Aug[max_r][j],Aug[row][j]);
        }
        if(sign(Aug[row][col])==0)  ///当前列row行以下全为0(包括row行)
        {
            row--;
            continue;
        }
        for(i = row+1; i < m; i++)
        {
            if(sign(Aug[i][col])==0)
                continue;
            double ta = Aug[i][col]/Aug[row][col];
            for(j = col; j < n+1; j++)
                Aug[i][j] -= Aug[row][j]*ta;
        }
    }
    for(i = row; i < m; i++)    ///col=n存在0...0,a的情况,无解
    {
        if(sign(Aug[i][col]))
            return -1;
    }
    if(row < n)     ///存在0...0,0的情况,有多个解,自由变元个数为n-row个
    {
        for(i = row-1; i >=0; i--)
        {
            int free_num = 0;   ///自由变元的个数
            int free_index;     ///自由变元的序号
            for(j = 0; j < n; j++)
            {
                if(sign(Aug[i][j])!=0 && free_x[j])
                    free_num++,free_index=j;
            }
            if(free_num > 1) continue;  ///该行中的不确定的变元的个数超过1个,无法求解,它们仍然为不确定的变元
            ///只有一个不确定的变元free_index,可以求解出该变元,且该变元是确定的
            double tmp = Aug[i][n];
            for(j = 0; j < n; j++)
            {
                if(sign(Aug[i][j])!=0 && j!=free_index)
                    tmp -= Aug[i][j]*x[j];
            }
            x[free_index] = tmp/Aug[i][free_index];
            free_x[free_index] = false;
        }
        return n-row;
    }
    ///有且仅有一个解,严格的上三角矩阵(n==m)
    for(i = n-1; i >= 0; i--)
    {
        double tmp = Aug[i][n];
        for(j = i+1; j < n; j++)
            if(sign(Aug[i][j])!=0)
                tmp -= Aug[i][j]*x[j];
        x[i] = tmp/Aug[i][i];
    }
    return 0;
}
///模板结束

int main()
{
    int i,j;
    int t;
    double a[12][12];
    scanf("%d",&t);
    while(t--)
    {
        memset(Aug,0.0,sizeof(Aug));
        memset(x,0.0,sizeof(x));
        memset(free_x,1,sizeof(free_x));    ///都是不确定的变元
        for(i = 0; i < 12; i++)
            for(j = 0; j < 12; j++)
                scanf("%lf",&a[i][j]);
        double sum=0;
        for(int i=0;i<11;i++)
            sum+=a[11][i]*a[11][i];
        for(int i=0;i<11;i++)
        {
              for(int j=0;j<11;j++)
              {
                  Aug[i][j]=2*(a[i][j]-a[11][j]);
                  Aug[i][11]+=a[i][j]*a[i][j];
              }
              Aug[i][11]+=-a[i][11]*a[i][11]+a[11][11]*a[11][11]-sum;
        }
        m = n = 11;
        Gauss();
        for(int i = 0; i < n; i++)
        {
            printf("%.2lf",x[i]);
            printf("%c",i==n-1?'\n':' ');
        }
    }
    return 0;
}

 

 

开关问题POJ1830

题意:有n(n<29)个开关,每个开关的状态用0/1表示关/开.
开关之间存在着某些联系,比如打开或关闭某个开关,与它相关的某些开关也会发生变化,即原来为开/关,现在变为关/开.
经过若干次操作后,由初始状态为s变为末端状态为t(s,t均为长度为n的01串).
问有多少种方法可以将s变为t,每个开关最多只能动一次,不考虑开关操作的顺序.
输入格式:
开关个数n
初始状态s和目标状态t(01串)
若干整数对(i,j),表示对开关i的操作,开关j的状态也会发生改变.

解:
利用矩阵思想A*x=b,其中A为n*n的变换矩阵,x,b均为列矩阵;
x 表示对n个开关的操作,1表示动它,0表示不动它,显然最后是求解的个数.
A[i][j]=1表示若操作开关j会引起i的变化.
b 表示若开关i的初始状态和目标状态不同,则b[i]=1,否则b[i]=0.即b表示最终所有开关的变化状态.0为不变,1为变.
b[i] = ^(A[i][j]*x[j]),0<=j<=n.

**********理解*********
考虑最终开关i状态的变化:
j!=i
若x[j]=1, 表示对开关j进行操作
    若A[i][j]=1,表示开关j对i有影响,即A[i][j]*x[j]=1,此时开关i状态改变;
    若A[i][j]=0,表示虽然对j操作,但开关j对i却没有影响,即A[i][j]*x[j]=0,此时开关i状态不变;
若x[j]=0, 表示不对开关j操作,显然对i的状态也没有影响,即A[i][j]*x[j]=0,此时开关i状态不变.

异或方程组的消元类似于普通方程组的消元:

若有ax + by = p,cx + dy = q,则有(a+c)x + (b+d)y = p+q;

同理:若有ax ^ by = p,cx ^ dy = q,则有(a^c)x ^ (b^d)y = p^q。

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

///高斯消元模板
const int Max_M = 35;
const int Max_N = 35;
int Aug[Max_M][Max_N];
int m,n;        ///m个方程,n个未知数

int Gauss()
{
    int i,j;
    int max_r,row,col;
    for(col=0,row=0; row<m && col<n; row++,col++)
    {
        max_r = row;
        for(i = row+1; i < m; i++)
        {
            if(Aug[i][col] > Aug[max_r][col])
                max_r = i;
        }
        if(max_r != row)
        {
            for(j = row; j < n+1; j++)
                swap(Aug[max_r][j],Aug[row][j]);
        }
        if(Aug[row][col]==0)
        {
            row--;
            continue;
        }
        for(i = row+1; i < m; i++)
        {
            if(Aug[i][col]==0)
                continue;
            for(j = col; j < n+1; j++)
                Aug[i][j] ^= Aug[row][j];
        }
    }
    for(i = row; i < m; i++)
    {
        if(Aug[i][col] != 0)
            return -1;
    }
    return 1<<(n-row);          ///n-row个变元,只有0/1两种取值
}

///end

int main()
{
    int i,j;
    int t;
    int nn;
    int start[30],end[30];
    scanf("%d",&t);
    while(t--)
    {
        m = n = 0;
        memset(Aug,0,sizeof(Aug));
        scanf("%d",&nn);
        for(i = 0; i < nn; i++)
            scanf("%d",&start[i]);
        for(i = 0; i < nn; i++)
            scanf("%d",&end[i]);
        while(true)
        {
            scanf("%d %d",&i,&j);
            if(!i&&!j)
                break;
            Aug[j-1][i-1] = 1;
        }
        n = m = nn;
        for(i = 0; i < n; i++)
        {
            Aug[i][i] = 1;
            Aug[i][n] = start[i]^end[i];
        }
        int ret = Gauss();
        if(ret==-1)
            printf("Oh,it's impossible~!!\n");
        else
            printf("%d\n",ret);
    }
    return 0;
}

 

 

     ps: 其实当时看到这题是用双向枚举过的,正向和逆向的复杂度为O(2^(n/2)),对于n比较小的也足够了,大了还是得用高斯消元解异或方程。

#include <iostream>
#include <cstdio>
#include <cstring>
#include <map>
#include <queue>
#include <algorithm>
using namespace std;
const int MaxN = 30;
int nn;                 ///开关数<29
int start,end;
int ch[MaxN];           ///chp[i]表示拨动开关i,会引发的状态改变
bool flag[MaxN][MaxN];  ///防止输入对IJ重复
map<int,int> vis;       ///记录状态及出现的次数
map<int,int>::iterator mit;
int result;             ///次数

///枚举1..mid开关的所有情况,共2^mid种
void Forward(int mid)
{
    int i,j,k;
    int Max = 1<<mid;
    for(i = 0; i < Max; i++)
    {
        j = i;
        k = 0;
        int now = start;
        while(j)
        {
            if(j&1)
                now = now ^ ch[k];
            j >>= 1;
            k++;
        }
        mit = vis.find(now);
        if(mit == vis.end())
            vis.insert(make_pair(now,1));
        else
            mit->second++;
    }
}

void Backward(int mid)
{
    int i,j,k;
    int Max = 1<<(nn-mid);
    for(i = 0; i < Max; i++)
    {
        j = i;
        k = mid;
        int now = end;
        while(j)
        {
            if(j&1)
                now = now ^ ch[k];
            j >>= 1;
            k++;
        }
        mit = vis.find(now);
        if(mit!=vis.end())
            result += mit->second;
    }
}

int main()
{
    int i,j;
    int t;
    scanf("%d",&t);
    while(t--)
    {
        result = 0;
        start = end = 0;
        vis.clear();
        memset(flag,0,sizeof(flag));
        scanf("%d",&nn);
        for(i = 0; i < nn; i++)
            ch[i] = 1<<i,flag[i][i]=true;
        for(i = 0; i < nn; i++)
        {
            scanf("%d",&j);
            if(j) start += 1<<i;
        }
        for(i = 0; i < nn; i++)
        {
            scanf("%d",&j);
            if(j) end += 1<<i;
        }
        while(true)
        {
            scanf("%d %d",&i,&j);
            if(!i&&!j)
                break;
            if(!flag[i-1][j-1])
                ch[i-1] += 1<<(j-1),flag[i-1][j-1]=true;
        }
        Forward(nn/2);
        Backward(nn/2);
        if(result==0)
            printf("Oh,it's impossible~!!\n");
        else
            printf("%d\n",result);
    }
    return 0;
}

 

  • 大小: 21.3 KB
1
0
分享到:
评论

相关推荐

Global site tag (gtag.js) - Google Analytics