博客
关于我
强烈建议你试试无所不能的chatGPT,快点击我
hdu 5755(高斯消元——模线性方程组模板)
阅读量:6827 次
发布时间:2019-06-26

本文共 7631 字,大约阅读时间需要 25 分钟。

PS. 看了大神的题解,发现确实可以用m个未知数的高斯消元做。因为确定了第一行的情况,之后所有行的情况都可以根据第一行推。 这样复杂度直接变成O(m*m*m)

 

知道了是高斯消元后,其实只要稍加处理,就可以解决带模的情况。

1 是在进行矩阵行变化的时候,取模。

2 最后的除法用逆元。(因为a[i][i]必定非0 且小于模数) 

然后对于无穷多解的情况,只需要将那些列全为0的未知数定义一个固定值。(这里设的是0)其余操作不变。

#include 
#include
#include
#include
#include
#include
#include
#include
#include
#include
#define CL(a,num) memset((a),(num),sizeof(a))#define iabs(x) ((x) > 0 ? (x) : -(x))#define Min(a,b) (a) > (b)? (b):(a)#define Max(a,b) (a) > (b)? (a):(b)#define ll long long#define inf 0x7f7f7f7f#define MOD 100000007#define lc l,m,rt<<1#define rc m + 1,r,rt<<1|1#define pi acos(-1.0)#define test puts("<------------------->")#define maxn 100007#define M 100007#define N 1010using namespace std;//freopen("din.txt","r",stdin);int a[N][N],X[N];//分别记录增广矩阵和解集int free_x[N];//记录自由变量int equ,var;//分别表示方程组的个数和变量的个数int GCD(int x,int y){ if (y == 0) return x; return GCD(y,x%y);}int LCM(int x,int y){ return x/GCD(x,y)*y;}void Debug(void){ int i, j; for (i = 0; i < equ; i++) { for (j = 0; j < var + 1; j++) { cout << a[i][j] << " "; } cout << endl; } cout << endl;}// 高斯消元法解方程组(Gauss-Jordan elimination).(-2表示有浮点数解,//但无整数解,-1表无整数解,-1表示无解,0表示唯一解,大于0表示无穷解,//并返回自由变元的个数)//ax + by = gcd(a,b)//传入固定值a,b.放回 d=gcd(a,b), x , yvoid extendgcd(ll a,ll b,ll &d,ll &x,ll &y){ if(b==0){d=a;x=1;y=0;return;} extendgcd(b,a%b,d,y,x); y-=x*(a/b);}//Ax=1(mod M),gcd(A,M)==1//输入:10^18>=A,M>=1//输出:返回x的范围是[1,M-1]ll GetNi(ll A,ll MM){ ll rex=0,rey=0; ll td=0; extendgcd(A,MM,td,rex,rey); return (rex%MM+MM)%MM;}int Guass(){ int i,j,k,col; CL(X,0); CL(free_x,1); for (k = 0,col = 0; k < equ && col < var; k++, col++){ int max_r = k; for (i = k + 1; i < equ; ++i){ if (iabs(a[i][col]) > iabs(a[max_r][col])) max_r = i; } if (max_r != k){ for (i = k; i < var + 1; ++i) swap(a[k][i],a[max_r][i]); } if (a[k][col] == 0){ //可以随意定义的变量 X[col] = 0;//强制赋值为0 free_x[col] = 0; k--; //cout<
<
= 0; --i){ num = 0; int tmp = a[i][var]; for (j = 0; j < var; ++j){ if (a[i][j] != 0 && free_x[j]){ num++; freeidx = j; } } if (num > 1) continue; tmp = a[i][var]; for (j = 0; j < var; ++j){ if (a[i][j] && j != freeidx) tmp -= a[i][j]*X[j]; } //这里也要 int k2 = (tmp%3+3)%3; int k1 = (a[i][freeidx]%3+3)%3; if(k1!=0) { X[freeidx] = k2*(int)GetNi(k1, 3); } else { X[freeidx] = 0; //printf("X[%d]为任意?\n",i); } //X[freeidx] = tmp/a[i][freeidx]; free_x[freeidx] = 0; } return var - k; } // 3. 唯一解的情况: for (i = k - 1; i >= 0; --i){ int tmp = a[i][var]; for (j = i + 1; j < var; ++j){ tmp -= a[i][j]*X[j]; } //强行搞一发 int k2 = (tmp%3+3)%3; int k1 = (a[i][i]%3+3)%3; if(k1!=0) { X[i] = k2*(int)GetNi(k1, 3); } else { X[i] = 0; //printf("X[%d]为任意?\n",i); } //X[i] = tmp/a[i][i];//不整除? } return 0;}int g[33][33];int getid(int i,int j,int n,int m){ return (i-1)*m + (j-1);}int up[4]={ 1,-1,0,0};int rl[4]={ 0,0,1,-1};int main(){ //freopen("din.txt","r",stdin); int i,j; int T; cin>>T; while(T--) { int n,m; scanf("%d%d",&n,&m); for(i=1;i<=n;i++) for(j=1;j<=m;j++) scanf("%d",&g[i][j]); equ = n*m; var = n*m; memset(a,0,sizeof(a)); int id = 0; for(i=1;i<=n;i++) for(j=1;j<=m;j++) { for(int k=0;k<4;k++) { int ti = i+up[k]; int tj = j+rl[k]; if( ti>=1&&ti<=n && tj>=1&&tj<=m ) { int tid = getid(ti, tj, n, m); a[id][tid] = 1; } } a[id][id] = 2; a[id][n*m] = (3-g[i][j])%3; id++; } CL(X,0); CL(free_x,0); //for (i = 0; i < equ; ++i) //for (j = 0; j < var + 1; ++j) scanf("%d",&a[i][j]); //Debug(); int free_num = Guass(); if (free_num == -1) printf("无解!\n"); else if (free_num == -2) printf("无整数解\n"); else {// else if (free_num > 0){ //printf("无穷多解! 自由变元个数为%d\n", free_num); //我懂了,我需要确定free_num个数// for (i = 0; i < var; ++i){// if (free_x[i]) printf("X%d 是不确定的\n",i + 1);// else printf("X%d %d\n",i + 1,X[i]);// } // }// else{ //我觉得可以! int ans = 0; for (i = 0 ; i < var; ++i){ if(X[i]%3 != 0) ans += (X[i]%3+3)%3; //printf("X%d %d\n",i + 1,X[i]); } printf("%d\n",ans); for(i=0;i

 

2. m个未知数的情况

////  main.cpp//  hdu5755.1////  Created by New_Life on 16/8/4.//  Copyright © 2016年 chenhuan001. All rights reserved.//#include 
#include
#include
#include
#include
#include
#include
#include
#include
#define CL(a,num) memset((a),(num),sizeof(a))#define iabs(x) ((x) > 0 ? (x) : -(x))#define Min(a,b) (a) > (b)? (b):(a)#define Max(a,b) (a) > (b)? (a):(b)#define ll long long#define inf 0x7f7f7f7f#define MOD 100000007#define lc l,m,rt<<1#define rc m + 1,r,rt<<1|1#define pi acos(-1.0)#define test puts("<------------------->")#define maxn 100007#define M 100007#define N 33using namespace std;//freopen("din.txt","r",stdin);int a[N][N],X[N];//分别记录增广矩阵和解集int free_x[N];//记录自由变量int equ,var;//分别表示方程组的个数和变量的个数int GCD(int x,int y){ if (y == 0) return x; return GCD(y,x%y);}int LCM(int x,int y){ return x/GCD(x,y)*y;}void Debug(void){ int i, j; for (i = 0; i < equ; i++) { for (j = 0; j < var + 1; j++) { cout << a[i][j] << " "; } cout << endl; } cout << endl;}// 高斯消元法解方程组(Gauss-Jordan elimination).(-2表示有浮点数解,//但无整数解,-1表无整数解,-1表示无解,0表示唯一解,大于0表示无穷解,//并返回自由变元的个数)//ax + by = gcd(a,b)//传入固定值a,b.放回 d=gcd(a,b), x , yvoid extendgcd(ll a,ll b,ll &d,ll &x,ll &y){ if(b==0){d=a;x=1;y=0;return;} extendgcd(b,a%b,d,y,x); y-=x*(a/b);}//Ax=1(mod M),gcd(A,M)==1//输入:10^18>=A,M>=1//输出:返回x的范围是[1,M-1]ll GetNi(ll A,ll MM){ ll rex=0,rey=0; ll td=0; extendgcd(A,MM,td,rex,rey); return (rex%MM+MM)%MM;}int Guass(){ int i,j,k,col; CL(X,0); CL(free_x,1); for (k = 0,col = 0; k < equ && col < var; k++, col++){ int max_r = k; for (i = k + 1; i < equ; ++i){ if (iabs(a[i][col]) > iabs(a[max_r][col])) max_r = i; } if (max_r != k){ for (i = k; i < var + 1; ++i) swap(a[k][i],a[max_r][i]); } if (a[k][col] == 0){ //可以随意定义的变量 X[col] = 0;//强制赋值为0 free_x[col] = 0; k--; //cout<
<
= 0; --i){ num = 0; int tmp = a[i][var]; for (j = 0; j < var; ++j){ if (a[i][j] != 0 && free_x[j]){ num++; freeidx = j; } } if (num > 1) continue; tmp = a[i][var]; for (j = 0; j < var; ++j){ if (a[i][j] && j != freeidx) tmp -= a[i][j]*X[j]; } //这里也要 int k2 = (tmp%3+3)%3; int k1 = (a[i][freeidx]%3+3)%3; if(k1!=0) { X[freeidx] = k2*(int)GetNi(k1, 3); } else { X[freeidx] = 0; //printf("X[%d]为任意?\n",i); } //X[freeidx] = tmp/a[i][freeidx]; free_x[freeidx] = 0; } return var - k; } // 3. 唯一解的情况: for (i = k - 1; i >= 0; --i){ int tmp = a[i][var]; for (j = i + 1; j < var; ++j){ tmp -= a[i][j]*X[j]; } //强行搞一发 int k2 = (tmp%3+3)%3; int k1 = (a[i][i]%3+3)%3; if(k1!=0) { X[i] = k2*(int)GetNi(k1, 3); } else { X[i] = 0; //printf("X[%d]为任意?\n",i); } //X[i] = tmp/a[i][i];//不整除? } return 0;}int g[33][33];int save[33][33][33];int getni(int x){ if(x==1) return 2; if(x==2) return 1; return 0;}void func(int x,int y){ g[x][y] = (g[x][y]+2)%3; g[x+1][y] = (g[x+1][y]+1)%3; g[x][y+1] = (g[x][y+1]+1)%3; g[x-1][y] = (g[x-1][y]+1)%3; g[x][y-1] = (g[x][y-1]+1)%3;}int main() { int T; cin>>T; while(T--) { int n,m; scanf("%d%d",&n,&m); for(int i=1;i<=n;i++) for(int j=1;j<=m;j++) scanf("%d",&g[i][j]); memset(save,0,sizeof(save)); for(int i=1;i<=m;i++) { save[1][i][i] = 1; } for(int i=2;i<=n;i++) { for(int j=1;j<=m;j++) { for(int k=0;k<=m;k++) { int tmp = (save[i-1][j][k]*2+save[i-1][j+1][k]+save[i-1][j-1][k]) + save[i-2][j][k]; if(k == 0) tmp += g[i-1][j]; tmp %= 3; save[i][j][k] = getni(tmp); } } } //然后构建矩阵 memset(a,0,sizeof(a)); equ = m; var = m; for(int j=1;j<=m;j++) { for(int k=0;k<=m;k++) { int tmp = save[n][j][k]*2+save[n][j+1][k]+save[n][j-1][k] + save[n-1][j][k]; if(k==0) { tmp += g[n][j]; a[j-1][m] = getni(tmp%3); } else a[j-1][k-1] = tmp%3; } } CL(X,0); CL(free_x,0); int free_num = Guass(); if(free_num == -1 || free_num == -2){ cout<
<<" 错误"<

 

转载地址:http://zcykl.baihongyu.com/

你可能感兴趣的文章
WinForm更新文件
查看>>
setprecision **fixed
查看>>
JVM系列五:JVM监测&工具[整理中]
查看>>
局部自适应自动色阶/对比度算法在图像增强上的应用。
查看>>
CMD命令
查看>>
Spring中@Autowired与@Resource的区别
查看>>
Python 学习笔记 -- 类和实例
查看>>
Android 静默安装/后台安装
查看>>
java 非空判断类
查看>>
【html5】如何让Canvas标签自适应设备
查看>>
SecureCRT最佳配色方法+直接修改默认配置方法 - imsoft.cnblogs
查看>>
通俗地介绍下---数据结构之堆
查看>>
JQuery实现简单的服务器轮询效果
查看>>
2017.6.26 工作记录
查看>>
“Too many open files” 小记
查看>>
《Effective C#》读书笔记——条目4:使用Conditional特性而不是#if条件编译<C#语言习惯>...
查看>>
浅谈异常与恋爱
查看>>
分享:http-watcher更新,改进对动态web程序的支持
查看>>
设计模式---->经典设计模式一览
查看>>
Asp.Net生命周期系列一
查看>>