消元 

 

高斯消元

1 typedef double Matrix[maxn][maxn]; 2 void gauss_elimination(Matrix A, int n){ 3     int i, j, k, r; 4     //消元过程 5     for(i = 0; i < n; ++i){ 6         //选一行r并与i行交换 7         r = i; 8         for(j = i+1; j < n; ++j){ 9             if(fabs(A[j][i]) > fabs(A[r][i])) r = j;10         }11         if(r != i){12             for(j = 0; j <= n; ++j) swap(A[r][j], A[i][j]);13         }14         //与i+1~n行进行消元15         for(k = i+1; k < n; ++k){16             for(j = n; j >= i; --j){ //必须逆序枚举17                 A[k][j] -= A[k][i]/A[i][i]*A[i][j];18             }19         }20     }21     //回代过程22     for(i = n-1; i >= 0; --i){23         for(j = i+1; j < n; ++j){24             A[i][n] -= A[j][n] * A[i][j];25         }26         A[i][n] /= A[i][i];27     }28 }

 

 

 

高斯——约当消元法

　　运算量比高斯消元略大（将系数矩阵化为对角矩阵），但是代码更简单（少了回调过程）

 

1 typedef double Matrix[maxn][maxn]; 2 const double eps = 1e-8; 3 void gauss_jordan(Matrix A, int n){ 4     int i, j, k, r; 5     for(i = 0; i < n; i++){ 6         r = i; 7         for(j = i+1; j < n; j++){ 8             if(fabs(A[j][i]) > fabs(A[r][i])) r = j; 9         }10         if(fabs(A[r][i]) < eps) continue;11         if(r != i){12             for(j = 0; j <= n; j++){13                 swap(A[r][j], A[i][j]);14             }15         }16         for(k = 0; k < n; k++) if(k != i){17             for(j = n; j >= i; j--){18                 A[k][j] -= A[k][i]/A[i][i] * A[i][j];19             }20         }21     }22 }

 

 

 

--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

 

https://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&category=20&problem=1769&mosmsg=Submission+received+with+ID+18927455

随机程序

 

求每个节点的期望执行次数

 

设　i　的出度为　di　期望执行次数为　xi　对于每个有三个前驱点　a　b　c　的节点　i　可以列出方程　xi　=　xa/da　+　xb/db　+　xc/dc　.

所以矩阵就可以构建出来

 

1 #include 
     
       2 #include 
      
        3 #include 
       
         4 #include 
        
          5 #include 
         
           6 using namespace std; 7  8 const double eps = 1e-8; 9 const int maxn = 110;10 typedef double Ma[maxn][maxn];11 12 Ma A;13 int n, d[maxn];14 vector
          
            p[maxn];15 int inf[maxn];16 17 18 void gass(Ma A, int n){19 int i, j, k, r;20 for(i = 0; i < n; ++i){21 r = i;22 for(j = i+1; j < n; ++j){23 if(fabs(A[j][i]) > fabs(A[r][i])) r = j;24 } 25 if(fabs(A[r][i]) < eps) continue;26 if(r != i){27 for(j = 0; j <= n; ++j) swap(A[r][j], A[i][j]);28 }29 for(k = 0; k < n; ++k) if(k != i){30 for(j = n; j >= i; j--) A[k][j] -= A[k][i] / A[i][i] * A[i][j];31 }32 }33 }34 35 36 int main(){37 int kase = 0;38 while(scanf("%d", &n) == 1 && n){39 memset(d, 0, sizeof(d));40 for(int i = 0; i < n; ++i) p[i].clear();41 int a, b;42 while(scanf("%d%d", &a, &b) == 2 && a){43 a--, b--;44 d[a]++;45 p[b].push_back(a);46 }47 memset(A, 0, sizeof(A));48 for(int i = 0; i < n; ++i){49 A[i][i] = 1;50 for(int j = 0; j < p[i].size(); j++){51 int xx = p[i][j];52 A[i][xx] -= 1.0 / d[xx];53 }54 if(i == 0) A[i][n] = 1;55 }56 gass(A, n);57 memset(inf, 0, sizeof(inf));58 for(int i = n-1; i >= 0; --i){59 if(fabs(A[i][i]) < eps && fabs(A[i][n]) > eps) inf[i] = 1;60 for(int j = i+1; j < n; j++){61 if(fabs(A[i][j]) > eps && inf[j]) inf[i] = 1;62 }63 }64 65 int q, u;66 scanf("%d", &q);67 printf("Case #%d:\n", ++kase);68 while(q--){69 scanf("%d", &u);70 u--;71 if(inf[u]) printf("infinity\n");72 else printf("%.3lf\n", fabs(A[u][u]) < eps ? 0.0 : A[u][n]/A[u][u]);73 }74 75 }76 return 0;77 }
          
         
        
       
      
     

 

 

--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

 

kuangbin大神模板

 

1 #include
     
        2 #include
      
         3 #include
       
          4 #include
        
           5 #include
         
            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
          
           abs(a[max_r][col])) max_r=i; 80 } 81 if(max_r!=k) 82 { // 与第k行交换. 83 for(j=k;j
           
            = 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("无解!\n");172 else if (free_num == -2) printf("有浮点数解，无整数解!\n");173 else if (free_num > 0)174 {175 printf("无穷多解! 自由变元个数为%d\n", free_num);176 for (i = 0; i < var; i++)177 {178 if (free_x[i]) printf("x%d 是不确定的\n", i + 1);179 else printf("x%d: %d\n", i + 1, x[i]);180 }181 }182 else183 {184 for (i = 0; i < var; i++)185 {186 printf("x%d: %d\n", i + 1, x[i]);187 }188 }189 printf("\n");190 }191 return 0;192 }
           
          
         
        
       
      
     

 

----------------------------------------------------------------------------------------------------------------

消元

 

1 #include
     
       2 #include
      
        3 #include
       
         4  5 using namespace std; 6 7 typedef __int64 lld; 8 9 lld a[205][205]; 10 11 int sign; 12 lld N,MOD; 13 void solved() 14 { 15 lld ans=1; 16 for(int i=0;i
       
      
     

 

----------------------------------------------------------------------------------------------------------------

 湘潭邀请赛

 

http://202.197.224.59/OnlineJudge2/index.php/Problem/read/id/1260

 

A-1 = A*/|A|

 

1 #include 
     
       2 #include 
      
        3 #include 
       
         4 using namespace std; 5 #define ll long long 6 const int Mod = 1e9+7; 7 int n; 8 int a[210][210], b[210][210]; 9 int mul(int x){10     int xx = Mod - 2, sum = 1;11     while(xx){12         if(xx&1) sum = 1LL * sum * x % Mod;13         x = 1LL * x * x % Mod;14         xx >>= 1;15     }16     return sum;17 }18 void d(){19     int cnt = 1;20     for(int i = 0; i < n; ++i){21         for(int j = 0; j < n; ++j){22             b[i][j] = (i == j);23         }24     }25     for(int i = 0; i < n; ++i){26         int t = i;27         for(int j = i; j < n; ++j){28             if(!a[j][i]){29                 t = j;30             }31         }32         if(t != i) cnt *= -1;33         for(int j = 0; j < n; ++j){34             swap(a[i][j], a[t][j]);35             swap(b[i][j], b[t][j]);36         }37         cnt = 1LL * cnt * a[i][i] % Mod;38         int xx = mul(a[i][i]);  //求逆39         for(int j = 0; j < n; ++j){40             a[i][j] = 1LL * a[i][j]*xx%Mod;41             b[i][j] = 1LL * b[i][j]*xx%Mod;42         }43         for(int k = 0; k < n; ++k){44             if(k == i) continue;45             ll tm = a[k][i];46             for(int j = 0; j < n; ++j){47                 a[k][j] = (a[k][j] - 1LL * tm*a[i][j]%Mod + Mod)%Mod;48                 b[k][j] = (b[k][j] - 1LL * tm*b[i][j]%Mod + Mod)%Mod;49             }50         }51     }52     cnt = (Mod+cnt)%Mod;53     for(int i = 0; i < n; ++i){54         for(int j = 0; j < n; ++j){55             b[i][j] = 1LL * b[i][j]*cnt%Mod;56         }57     }58 }59 int main(){60     while(~scanf("%d", &n)){61         for(int i = 0; i < n; ++i){62             a[0][i] = 1;63         }64         for(int i = 1; i < n; ++i){65             for(int j = 0; j < n; ++j){66                 scanf("%d", &a[i][j]);67             }68         }69         d();70         for(int i = 0; i < n; ++i){71             printf("%d%c", (i&1 ? (Mod-b[i][0])%Mod : b[i][0]), (i == n-1 ? '\n' : ' '));72         }73     }74     return 0;75 }
       
      
     

 

 

 

----------------------------------------------------------------------------------------------------------------

 

https://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&category=27&page=show_problem&problem=2537

 

 乘积是平方数

 

1 #include 
     
       2 #include 
      
        3 using namespace std; 4 #define ll long long 5 const int maxn = 555; 6 int num = 0, n, ma; 7 bool vis[maxn]; 8 int pr[maxn], A[185][105]; 9 void p(){10     for(int i = 2; i <= 500; ++i){11         if(!vis[i]) pr[num++] = i;12         for(int j = 0; j < num && i*pr[j] <= 500; ++j){13             if(i%pr[j]) vis[i*pr[j]] = true;14             else{15                 vis[i*pr[j]] = true;16                 break;17             }18         }19     }20 }21 void solve(){22     int i = 0, j = 0, k, r, u;23     while(i < ma && j < n){24         r = i;25         for(k = i; k < ma; ++k){26             if(A[k][j]){27                 r = k;28                 break;29             }30         }31         if(A[r][j]){32             if(r != i){33                 for(k = 0; k <= n; ++k) swap(A[r][k], A[i][k]);34             }35             for(k = i+1; k < ma; ++k) if(A[k][j]){36                 for(u = i; u <= n; ++u){37                     A[k][u] ^= A[i][u];38                 }39             }40             ++i;41         }42         ++j;43     }44     //cout<<"****************"<
       <
      
     

 

 

 

----------------------------------------------------------------------------------------------------------------

----------------------------------------------------------------------------------------------------------------

----------------------------------------------------------------------------------------------------------------

 

 

只有不断学习才能进步！