高斯消元模板

来源：互联网发布：什么赚钱软件最好编辑：程序博客网时间：2024/06/07 01:55
#include <iostream>#include <string>#include <string.h>#include <vector>#include <stdio.h>#include <algorithm>#include <map>#include <math.h>typedef long long LL ;const int N = 1008 ;//高斯消元模板const double eps = 1e-12;double Aug[N][N]; //增广矩阵bool free_x[N];  //判断是否是不确定的变元double x[N]; //解集int sign(double x){ return (x>eps) - (x<-eps);}/**返回值:-1 无解0 有且仅有一个解>=1 有多个解,根据free_x判断哪些是不确定的解*/int Gauss(int equ , int var){    int i,j;    int row,col,max_r;    for(row=0,col=0; row<equ&&col<var; row++,col++)    {        max_r = row;        for(i = row+1; i < equ; i++)  //找到当前列所有行中的最大值(做除法时减小误差)        {            if(sign(fabs(Aug[i][col])-fabs(Aug[max_r][col]))>0)                max_r = i;        }        if(max_r != row)            //将该行与当前行交换        {            for(j = row; j < var+1; j++)                std::swap(Aug[max_r][j],Aug[row][j]);        }        if(sign(Aug[row][col])==0)  //当前列row行以下全为0(包括row行)        {            row--;            continue;        }        for(i = row+1; i < equ; i++)        {            if(sign(Aug[i][col])==0)                continue;            double ta = Aug[i][col]/Aug[row][col];            for(j = col; j < var+1; j++)                Aug[i][j] -= Aug[row][j]*ta;        }    }    for(i = row; i < equ; i++)    //col=n存在0...0,a的情况,无解    {        if(sign(Aug[i][col]))            return -1;    }    if(row < var)     //存在0...0,0的情况,有多个解,自由变元个数为n-row个    {        for(i = row-1; i >=0; i--)        {            int free_num = 0;   //自由变元的个数            int free_index;     //自由变元的序号            for(j = 0; j < var; j++)            {                if(sign(Aug[i][j])!=0 && free_x[j])                    free_num++,free_index=j;            }            if(free_num > 1) continue;  //该行中的不确定的变元的个数超过1个,无法求解,它们仍然为不确定的变元            //只有一个不确定的变元free_index,可以求解出该变元,且该变元是确定的            double tmp = Aug[i][var];            for(j = 0; j < var; j++)            {                if(sign(Aug[i][j])!=0 && j!=free_index)                    tmp -= Aug[i][j]*x[j];            }            x[free_index] = tmp/Aug[i][free_index];            free_x[free_index] = false;        }        return var-row;    }    //有且仅有一个解,严格的上三角矩阵(n==m)    for(i = var-1; i >= 0; i--)    {        double tmp = Aug[i][var];        for(j = i+1; j < var; j++)            if(sign(Aug[i][j])!=0)                tmp -= Aug[i][j]*x[j];        x[i] = tmp/Aug[i][i];    }    return 0;}//模板结束double  Fix(double x){   // 将 -0.000 变成 0.00 ,结果输出Fix(x[])        if(fabs(x) < eps) return 0 ;        return x ;}int main(){    int k , t  , n ;    scanf("%d" , &t) ;    while(t--){          scanf("%d%d" , &n , &k)  ;          memset(Aug , 0 , sizeof(Aug)) ;          memset(x , 0 , sizeof(x)) ;          memset(free_x , 0 , sizeof(free_x)) ;          for(int i = 0 ; i < n ; i++){              if(i == k){                   Aug[i][i] = 1.0 ;                   continue ;              }              Aug[i][i] = 1.0 ;              Aug[i][(i+1)%n ] = -0.5 ;              Aug[i][(i-1+n)%n ] = -0.5 ;              Aug[i][n] = 1.0 ;          }          Gauss(n,n) ;          printf("%.4lf\n" , Fix(x[0]) ) ;    }    return 0 ;}
0 0