zoukankan      html  css  js  c++  java
  • 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 == 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 , y
    void 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<<k<<endl;
                continue;
            }
            for (i = k + 1; i < equ; ++i){
                if (a[i][col] != 0){
                    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] < 0) tb = -tb;
                    for (j = col; j < var + 1; ++j){
                        a[i][j] = ((ta*a[i][j] - tb*a[k][j])%3+3)%3;
                    }
                }
            }
        }
        //Debug();
        // 1. 无解的情况:
        for (i = k; i < equ; ++i){
            if (a[i][col] != 0) return -1;
        }
        // 2. 无穷解的情况:
        if (k < var){
            
            int num = 0,freeidx=0;
            for (i = k - 1; i >= 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]为任意?
    ",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]为任意?
    ",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("无解!
    ");
            else if (free_num == -2) printf("无整数解
    ");
            else {
    //        else if (free_num > 0){
            
                //printf("无穷多解! 自由变元个数为%d
    ", free_num);
                //我懂了,我需要确定free_num个数
    //            for (i = 0; i < var; ++i){
    //                if (free_x[i]) printf("X%d 是不确定的
    ",i + 1);
    //                else printf("X%d %d
    ",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
    ",i + 1,X[i]);
                }
                printf("%d
    ",ans);
                for(i=0;i<var;i++)
                {
                    int tmp =  (X[i]%3+3)%3;
                    for(j=0;j<tmp;j++)
                    {
                        int x,y;
                        x = i/m;
                        y = i%m;
                        x++; y++;
                        printf("%d %d
    ",x,y);
                    }
                }
            }
            //printf("
    ");
        }
        return 0;
    }
    /*
     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 &#169; 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 == 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 , y
    void 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<<k<<endl;
                continue;
            }
            for (i = k + 1; i < equ; ++i){
                if (a[i][col] != 0){
                    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] < 0) tb = -tb;
                    for (j = col; j < var + 1; ++j){
                        a[i][j] = ((ta*a[i][j] - tb*a[k][j])%3+3)%3;
                    }
                }
            }
        }
        //Debug();
        // 1. 无解的情况:
        for (i = k; i < equ; ++i){
            if (a[i][col] != 0) return -1;
        }
        // 2. 无穷解的情况:
        if (k < var){
            
            int num = 0,freeidx=0;
            for (i = k - 1; i >= 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]为任意?
    ",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]为任意?
    ",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<<free_num<<" 错误"<<endl; continue;}
            int ans = 0;
            int saveans[33][33];
            memset(saveans,0,sizeof(saveans));
            for(int i=1;i<=n;i++)
                for(int j=1;j<=m;j++)
                {
                    int tmp = 0;
                    for(int k=0;k<=m;k++)
                    {
                        if(k==0) tmp += save[i][j][0];
                        else tmp += save[i][j][k]*X[k-1];
                        tmp %= 3;
                    }
                    ans += tmp;
                    saveans[i][j] = tmp;
                }
            
            printf("%d
    ",ans);
            for(int i=1;i<=n;i++)
                for(int j=1;j<=m;j++)
                {
                    for(int k=0;k<saveans[i][j];k++)
                    {
                        printf("%d %d
    ",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("
    ");
    //        }
        }
        return 0;
    }
    /*
     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
     */
  • 相关阅读:
    MySQL基础_常见命令
    数据库概述
    Nginx学习笔记
    华为OSPF基础配置实验
    华为RIPv2实验
    华为三层交换实验
    web-debug-server
    花一天时间试玩vsphere6.7(EXSI)服务器版的vmware
    haproxy+keepalived练习
    mailx加163邮箱发邮件
  • 原文地址:https://www.cnblogs.com/chenhuan001/p/5735174.html
Copyright © 2011-2022 走看看