zoukankan      html  css  js  c++  java
  • 高斯消元模板

      1 #include<stdio.h>
      2 #include<algorithm>
      3 #include<iostream>
      4 #include<string.h>
      5 #include<math.h>
      6 using namespace std;
      7 
      8 const int MAXN=50;
      9 
     10 
     11 
     12 int a[MAXN][MAXN];//增广矩阵
     13 int x[MAXN];//解集
     14 bool free_x[MAXN];//标记是否是不确定的变元
     15 
     16 
     17 
     18 /*
     19 void Debug(void)
     20 {
     21     int i, j;
     22     for (i = 0; i < equ; i++)
     23     {
     24         for (j = 0; j < var + 1; j++)
     25         {
     26             cout << a[i][j] << " ";
     27         }
     28         cout << endl;
     29     }
     30     cout << endl;
     31 }
     32 */
     33 
     34 
     35 inline int gcd(int a,int b)
     36 {
     37     int t;
     38     while(b!=0)
     39     {
     40         t=b;
     41         b=a%b;
     42         a=t;
     43     }
     44     return a;
     45 }
     46 inline int lcm(int a,int b)
     47 {
     48     return a/gcd(a,b)*b;//先除后乘防溢出
     49 }
     50 
     51 // 高斯消元法解方程组(Gauss-Jordan elimination).(-2表示有浮点数解,但无整数解,
     52 //-1表示无解,0表示唯一解,大于0表示无穷解,并返回自由变元的个数)
     53 //有equ个方程,var个变元。增广矩阵行数为equ,分别为0到equ-1,列数为var+1,分别为0到var.
     54 int Gauss(int equ,int var)
     55 {
     56     int i,j,k;
     57     int max_r;// 当前这列绝对值最大的行.
     58     int col;//当前处理的列
     59     int ta,tb;
     60     int LCM;
     61     int temp;
     62     int free_x_num;
     63     int free_index;
     64 
     65     for(int i=0;i<=var;i++)
     66     {
     67         x[i]=0;
     68         free_x[i]=true;
     69     }
     70 
     71     //转换为阶梯阵.
     72     col=0; // 当前处理的列
     73     for(k = 0;k < equ && col < var;k++,col++)
     74     {// 枚举当前处理的行.
     75 // 找到该col列元素绝对值最大的那行与第k行交换.(为了在除法时减小误差)
     76         max_r=k;
     77         for(i=k+1;i<equ;i++)
     78         {
     79             if(abs(a[i][col])>abs(a[max_r][col])) max_r=i;
     80         }
     81         if(max_r!=k)
     82         {// 与第k行交换.
     83             for(j=k;j<var+1;j++) swap(a[k][j],a[max_r][j]);
     84         }
     85         if(a[k][col]==0)
     86         {// 说明该col列第k行以下全是0了,则处理当前行的下一列.
     87             k--;
     88             continue;
     89         }
     90         for(i=k+1;i<equ;i++)
     91         {// 枚举要删去的行.
     92             if(a[i][col]!=0)
     93             {
     94                 LCM = lcm(abs(a[i][col]),abs(a[k][col]));
     95                 ta = LCM/abs(a[i][col]);
     96                 tb = LCM/abs(a[k][col]);
     97                 if(a[i][col]*a[k][col]<0)tb=-tb;//异号的情况是相加
     98                 for(j=col;j<var+1;j++)
     99                 {
    100                     a[i][j] = a[i][j]*ta-a[k][j]*tb;
    101                 }
    102             }
    103         }
    104     }
    105 
    106   //  Debug();
    107 
    108     // 1. 无解的情况: 化简的增广阵中存在(0, 0, ..., a)这样的行(a != 0).
    109     for (i = k; i < equ; i++)
    110     { // 对于无穷解来说,如果要判断哪些是自由变元,那么初等行变换中的交换就会影响,则要记录交换.
    111         if (a[i][col] != 0) return -1;
    112     }
    113     // 2. 无穷解的情况: 在var * (var + 1)的增广阵中出现(0, 0, ..., 0)这样的行,即说明没有形成严格的上三角阵.
    114     // 且出现的行数即为自由变元的个数.
    115     if (k < var)
    116     {
    117         // 首先,自由变元有var - k个,即不确定的变元至少有var - k个.
    118         for (i = k - 1; i >= 0; i--)
    119         {
    120             // 第i行一定不会是(0, 0, ..., 0)的情况,因为这样的行是在第k行到第equ行.
    121             // 同样,第i行一定不会是(0, 0, ..., a), a != 0的情况,这样的无解的.
    122             free_x_num = 0; // 用于判断该行中的不确定的变元的个数,如果超过1个,则无法求解,它们仍然为不确定的变元.
    123             for (j = 0; j < var; j++)
    124             {
    125                 if (a[i][j] != 0 && free_x[j]) free_x_num++, free_index = j;
    126             }
    127             if (free_x_num > 1) continue; // 无法求解出确定的变元.
    128             // 说明就只有一个不确定的变元free_index,那么可以求解出该变元,且该变元是确定的.
    129             temp = a[i][var];
    130             for (j = 0; j < var; j++)
    131             {
    132                 if (a[i][j] != 0 && j != free_index) temp -= a[i][j] * x[j];
    133             }
    134             x[free_index] = temp / a[i][free_index]; // 求出该变元.
    135             free_x[free_index] = 0; // 该变元是确定的.
    136         }
    137         return var - k; // 自由变元有var - k个.
    138     }
    139     // 3. 唯一解的情况: 在var * (var + 1)的增广阵中形成严格的上三角阵.
    140     // 计算出Xn-1, Xn-2 ... X0.
    141     for (i = var - 1; i >= 0; i--)
    142     {
    143         temp = a[i][var];
    144         for (j = i + 1; j < var; j++)
    145         {
    146             if (a[i][j] != 0) temp -= a[i][j] * x[j];
    147         }
    148         if (temp % a[i][i] != 0) return -2; // 说明有浮点数解,但无整数解.
    149         x[i] = temp / a[i][i];
    150     }
    151     return 0;
    152 }
    153 int main(void)
    154 {
    155     freopen("in.txt", "r", stdin);
    156     freopen("out.txt","w",stdout);
    157     int i, j;
    158     int equ,var;
    159     while (scanf("%d %d", &equ, &var) != EOF)
    160     {
    161         memset(a, 0, sizeof(a));
    162         for (i = 0; i < equ; i++)
    163         {
    164             for (j = 0; j < var + 1; j++)
    165             {
    166                 scanf("%d", &a[i][j]);
    167             }
    168         }
    169 //        Debug();
    170         int free_num = Gauss(equ,var);
    171         if (free_num == -1) printf("无解!
    ");
    172    else if (free_num == -2) printf("有浮点数解,无整数解!
    ");
    173         else if (free_num > 0)
    174         {
    175             printf("无穷多解! 自由变元个数为%d
    ", free_num);
    176             for (i = 0; i < var; i++)
    177             {
    178                 if (free_x[i]) printf("x%d 是不确定的
    ", i + 1);
    179                 else printf("x%d: %d
    ", i + 1, x[i]);
    180             }
    181         }
    182         else
    183         {
    184             for (i = 0; i < var; i++)
    185             {
    186                 printf("x%d: %d
    ", i + 1, x[i]);
    187             }
    188         }
    189         printf("
    ");
    190     }
    191     return 0;
    192 }
  • 相关阅读:
    Annual Summary-2019
    (翻译向)kafka--简介篇1
    mysql 查看数据库大小
    nohup保证程序后台运行
    NLP(四):文本向量化word2vec
    NLP(三):关键词提取
    NLP(二):jieba高频词提取
    NLP(一):jieba分词
    爬虫(一)百度翻译
    【论文笔记】社交网络中的信息扩散分析及其应用研究
  • 原文地址:https://www.cnblogs.com/zxz666/p/10726077.html
Copyright © 2011-2022 走看看