hdu 5755(高斯消元——模线性方程组模板)

时间:2023-03-09 08:28:44
hdu 5755(高斯消元——模线性方程组模板)

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

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

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

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

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

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <queue>
#include <stack>
#include <set>
#include <map>
#include <string> #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 1010
using 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 == ) 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 = ; i < equ; i++)
{
for (j = ; j < var + ; 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 , y
void extendgcd(ll a,ll b,ll &d,ll &x,ll &y)
{
if(b==){d=a;x=;y=;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=,rey=;
ll td=;
extendgcd(A,MM,td,rex,rey);
return (rex%MM+MM)%MM;
} int Guass(){
int i,j,k,col;
CL(X,); CL(free_x,); for (k = ,col = ; k < equ && col < var; k++, col++){
int max_r = k;
for (i = k + ; 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 + ; ++i) swap(a[k][i],a[max_r][i]);
}
if (a[k][col] == ){
//可以随意定义的变量
X[col] = ;//强制赋值为0
free_x[col] = ;
k--;
//cout<<k<<endl;
continue;
}
for (i = k + ; i < equ; ++i){
if (a[i][col] != ){
int lcm = LCM(a[k][col],a[i][col]);
int ta = lcm/iabs(a[i][col]); int tb = lcm/iabs(a[k][col]);
if (a[i][col]*a[k][col] < ) tb = -tb;
for (j = col; j < var + ; ++j){
a[i][j] = ((ta*a[i][j] - tb*a[k][j])%+)%;
}
}
}
}
//Debug();
// 1. 无解的情况:
for (i = k; i < equ; ++i){
if (a[i][col] != ) return -;
}
// 2. 无穷解的情况:
if (k < var){ int num = ,freeidx=;
for (i = k - ; i >= ; --i){
num = ;
int tmp = a[i][var];
for (j = ; j < var; ++j){
if (a[i][j] != && free_x[j]){
num++;
freeidx = j;
}
}
if (num > ) continue;
tmp = a[i][var];
for (j = ; j < var; ++j){
if (a[i][j] && j != freeidx) tmp -= a[i][j]*X[j];
}
//这里也要 int k2 = (tmp%+)%;
int k1 = (a[i][freeidx]%+)%;
if(k1!=)
{
X[freeidx] = k2*(int)GetNi(k1, );
}
else
{
X[freeidx] = ;
//printf("X[%d]为任意?\n",i);
} //X[freeidx] = tmp/a[i][freeidx];
free_x[freeidx] = ;
}
return var - k;
}
// 3. 唯一解的情况:
for (i = k - ; i >= ; --i){
int tmp = a[i][var];
for (j = i + ; j < var; ++j){
tmp -= a[i][j]*X[j];
}
//强行搞一发
int k2 = (tmp%+)%;
int k1 = (a[i][i]%+)%;
if(k1!=)
{
X[i] = k2*(int)GetNi(k1, );
}
else
{
X[i] = ;
//printf("X[%d]为任意?\n",i);
}
//X[i] = tmp/a[i][i];//不整除?
}
return ;
} int g[][]; int getid(int i,int j,int n,int m)
{
return (i-)*m + (j-);
}
int up[]={,-,,};
int rl[]={,,,-}; 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=;i<=n;i++)
for(j=;j<=m;j++)
scanf("%d",&g[i][j]);
equ = n*m;
var = n*m; memset(a,,sizeof(a)); int id = ;
for(i=;i<=n;i++)
for(j=;j<=m;j++)
{
for(int k=;k<;k++)
{
int ti = i+up[k];
int tj = j+rl[k];
if( ti>=&&ti<=n && tj>=&&tj<=m )
{
int tid = getid(ti, tj, n, m);
a[id][tid] = ;
}
} a[id][id] = ;
a[id][n*m] = (-g[i][j])%;
id++;
}
CL(X,);
CL(free_x,);
//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 == -) printf("无解!\n");
else if (free_num == -) 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 = ;
for (i = ; i < var; ++i){
if(X[i]% != ) ans += (X[i]%+)%;
//printf("X%d %d\n",i + 1,X[i]);
}
printf("%d\n",ans);
for(i=;i<var;i++)
{
int tmp = (X[i]%+)%;
for(j=;j<tmp;j++)
{
int x,y;
x = i/m;
y = i%m;
x++; y++;
printf("%d %d\n",x,y);
}
}
}
//printf("\n");
}
return ;
}
/*
10
3 3
0 0 0
0 0 0
0 0 1
1 2
0 0
2 30
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
*/

2. m个未知数的情况

//
// main.cpp
// hdu5755.1
//
// Created by New_Life on 16/8/4.
// Copyright © 2016年 chenhuan001. All rights reserved.
// #include <iostream>
#include <string.h>
#include <stdio.h>
#include <algorithm>
#include <queue>
#include <stack>
#include <set>
#include <map>
#include <string> #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 33
using 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 == ) 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 = ; i < equ; i++)
{
for (j = ; j < var + ; 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 , y
void extendgcd(ll a,ll b,ll &d,ll &x,ll &y)
{
if(b==){d=a;x=;y=;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=,rey=;
ll td=;
extendgcd(A,MM,td,rex,rey);
return (rex%MM+MM)%MM;
} int Guass(){
int i,j,k,col;
CL(X,); CL(free_x,); for (k = ,col = ; k < equ && col < var; k++, col++){
int max_r = k;
for (i = k + ; 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 + ; ++i) swap(a[k][i],a[max_r][i]);
}
if (a[k][col] == ){
//可以随意定义的变量
X[col] = ;//强制赋值为0
free_x[col] = ;
k--;
//cout<<k<<endl;
continue;
}
for (i = k + ; i < equ; ++i){
if (a[i][col] != ){
int lcm = LCM(a[k][col],a[i][col]);
int ta = lcm/iabs(a[i][col]); int tb = lcm/iabs(a[k][col]);
if (a[i][col]*a[k][col] < ) tb = -tb;
for (j = col; j < var + ; ++j){
a[i][j] = ((ta*a[i][j] - tb*a[k][j])%+)%;
}
}
}
}
//Debug();
// 1. 无解的情况:
for (i = k; i < equ; ++i){
if (a[i][col] != ) return -;
}
// 2. 无穷解的情况:
if (k < var){ int num = ,freeidx=;
for (i = k - ; i >= ; --i){
num = ;
int tmp = a[i][var];
for (j = ; j < var; ++j){
if (a[i][j] != && free_x[j]){
num++;
freeidx = j;
}
}
if (num > ) continue;
tmp = a[i][var];
for (j = ; j < var; ++j){
if (a[i][j] && j != freeidx) tmp -= a[i][j]*X[j];
}
//这里也要 int k2 = (tmp%+)%;
int k1 = (a[i][freeidx]%+)%;
if(k1!=)
{
X[freeidx] = k2*(int)GetNi(k1, );
}
else
{
X[freeidx] = ;
//printf("X[%d]为任意?\n",i);
} //X[freeidx] = tmp/a[i][freeidx];
free_x[freeidx] = ;
}
return var - k;
}
// 3. 唯一解的情况:
for (i = k - ; i >= ; --i){
int tmp = a[i][var];
for (j = i + ; j < var; ++j){
tmp -= a[i][j]*X[j];
}
//强行搞一发
int k2 = (tmp%+)%;
int k1 = (a[i][i]%+)%;
if(k1!=)
{
X[i] = k2*(int)GetNi(k1, );
}
else
{
X[i] = ;
//printf("X[%d]为任意?\n",i);
}
//X[i] = tmp/a[i][i];//不整除?
}
return ;
} int g[][];
int save[][][]; int getni(int x)
{
if(x==) return ;
if(x==) return ;
return ;
} void func(int x,int y)
{
g[x][y] = (g[x][y]+)%;
g[x+][y] = (g[x+][y]+)%;
g[x][y+] = (g[x][y+]+)%;
g[x-][y] = (g[x-][y]+)%;
g[x][y-] = (g[x][y-]+)%;
} int main() {
int T;
cin>>T;
while(T--)
{
int n,m;
scanf("%d%d",&n,&m);
for(int i=;i<=n;i++)
for(int j=;j<=m;j++)
scanf("%d",&g[i][j]);
memset(save,,sizeof(save)); for(int i=;i<=m;i++)
{
save[][i][i] = ;
}
for(int i=;i<=n;i++)
{
for(int j=;j<=m;j++)
{
for(int k=;k<=m;k++)
{
int tmp = (save[i-][j][k]*+save[i-][j+][k]+save[i-][j-][k]) + save[i-][j][k];
if(k == ) tmp += g[i-][j];
tmp %= ;
save[i][j][k] = getni(tmp);
}
}
} //然后构建矩阵
memset(a,,sizeof(a));
equ = m;
var = m;
for(int j=;j<=m;j++)
{
for(int k=;k<=m;k++)
{
int tmp = save[n][j][k]*+save[n][j+][k]+save[n][j-][k] + save[n-][j][k];
if(k==)
{
tmp += g[n][j];
a[j-][m] = getni(tmp%);
}
else a[j-][k-] = tmp%;
}
}
CL(X,);
CL(free_x,);
int free_num = Guass();
if(free_num == - || free_num == -){ cout<<free_num<<" 错误"<<endl; continue;}
int ans = ;
int saveans[][];
memset(saveans,,sizeof(saveans));
for(int i=;i<=n;i++)
for(int j=;j<=m;j++)
{
int tmp = ;
for(int k=;k<=m;k++)
{
if(k==) tmp += save[i][j][];
else tmp += save[i][j][k]*X[k-];
tmp %= ;
}
ans += tmp;
saveans[i][j] = tmp;
} printf("%d\n",ans);
for(int i=;i<=n;i++)
for(int j=;j<=m;j++)
{
for(int k=;k<saveans[i][j];k++)
{
printf("%d %d\n",i,j);
//func(i,j);
}
} // for(int i=1;i<=n;i++)
// {
// for(int j=1;j<=m;j++) printf("%d ",g[i][j]);
// printf("\n");
// }
}
return ;
}
/*
10
1 3
0 0 1
3 3
0 0 0
0 0 0
0 0 1
1 2
0 0
2 30
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
*/