您的位置:首页 > 其它

高斯消元

2014-08-24 22:13 218 查看
#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
using namespace std;

int a[40][40],x[40];//系数矩阵,解数组
int equ,var;//方程数,变元数

void Debug()
{
for(int i=0;i<equ;i++)
{
for(int j=0;j<=var;j++)
printf("%d ",a[i][j]);
puts("");
}
puts("");
}

int gcd(int a,int b)  //最大公约数
{
if(b==0) return a;
if(a<0) return gcd(-a,b);
if(b<0) return gcd(a,-b);
return gcd(b,a%b);
}

int lcm(int a,int b)  //最小公倍数
{
return a/gcd(a,b)*b;
}

int Gauss()
{
int k,col,max_r;//k为当前处理的方程的下标,col为x的下标,max_r为最大元素行
//求出max_r,并与当前行交换是为了求解时减小误差

for(k=0,col=0;k<equ&&col<var;k++,col++)      //处理所有的方程
{
max_r=k;
for(int i=k+1;i<equ;i++)                  //求出当前列元素最大的那一行
if(a[i][col]>a[max_r][col]) max_r=i;

if(max_r!=k)                              //交换两行的元素
for(int i=0;i<=var;i++)
{
swap(a[max_r][i],a[k][i]);
}

if(a[k][col]==0)                          //若为0,则这列所有元素都为0,无需处理
{
k--;
continue;
}

for(int i=k+1;i<equ;i++)                  //不为0时,将所有行当前列的元素置0,并处理所有其他元素(这里是二进制)
{
if(a[i][col])
for(int j=col;j<=var;j++)
{
a[i][j]=a[i][j]^a[k][j];
}
}

/*                                        //通用性处理各行元素
for(int i=k+1;i<equ;i++)
{
if(a[i][col])
{
int LCM=lcm(a[k][col],a[i][col]);
int ta=LCM/a[i][col],tb=LCM/a[k][col];

if(ta*tb<0) tb=-tb;
for(int j=col;j<=var;j++)
{
a[i][j]=(a[i][j]*ta%2-a[k][j]*tb%2+2)%2;
}
}
}
*/
}

Debug();

for(int i=k;i<equ;i++)                      //无解的情况
if(a[i][var]) return -1;

for(int i=0;i<equ;i++)                     //处理系数矩阵对角线元素均不为0
if(!a[i][i])
{
int j;
for(j=i+1;j<var;j++)
if(a[i][j]) break;

if(j>=var) break;

for(int r=0;r<equ;r++)
swap(a[r][i],a[r][j]);
}

if(k<var)                                  //多解情况,下面分两种情况讨论(二进制枚举,全部设置为0)
{
//二进制枚举
int lim=(1<<(var-k)),ans=10000;
for(int i=0;i<lim;i++)
{
int temp=i,sum=0;
for(int j=var-1;j>=k;j--)
{
x[j]=temp%2;temp/=2;
}
for(int j=k-1;j>=0;j--)
{
temp=a[j][var]%2;
for(int r=j+1;r<var;r++)
temp=(temp-a[j][r]*x[r]%2+2)%2;

x[j]=(temp/a[j][j])%2;
}
for(int j=0;j<var;j++)
if(x[j]!=0)
{
sum++;
}
ans=min(ans,sum);
}
return ans;

/*  全部设置为0
for(int j=var-1;j>=var-k;j--)
{
x[j]=0;
}
for(int j=var-k-1;j>=0;j--)
{
temp=a[j][var];
for(int r=j+1;r<var;r++)
temp-=a[j][r]*x[r];

x[j]=temp/a[j][j];
}
return var-k;
*/
}

int ans=0;                                //唯一解的情况
for(int j=var-1;j>=0;j--)
{
int temp=a[j][var]%2;
for(int r=j+1;r<var;r++)
temp=(temp-a[j][r]*x[r]%2+2)%2;

x[j]=(temp/a[j][j])%2;
}
for(int j=0;j<var;j++)
if(x[j]!=0)
{
ans++;
}
return ans;
}

void init()               //系数矩阵初始化
{
memset(a,0,sizeof(a));
memset(x,0,sizeof(x));

}

int main()
{

return 0;
}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: