Cod sursa(job #1942498)

Utilizator stoianmihailStoian Mihail stoianmihail Data 28 martie 2017 00:36:32
Problema Algoritmul lui Gauss Scor 0
Compilator cpp Status done
Runda Arhiva educationala Marime 5.03 kb
#include <stdio.h>
#include <math.h>
#include <assert.h>
#include <algorithm>

using std::swap;

#define MAX_N 301
#define NIL -1

const int even[2] = {-NIL, NIL};
const double eps = 1e-7;

int M, N;
int init[MAX_N][MAX_N], expand[MAX_N][MAX_N];
int move[MAX_N][MAX_N], tmp[MAX_N][MAX_N];
int out[MAX_N][MAX_N + 1];
double freq[MAX_N + 1];
double l[MAX_N][MAX_N], u[MAX_N][MAX_N], a[MAX_N][MAX_N], b[MAX_N][MAX_N];

int rank(int mtr[MAX_N][MAX_N], int n, int m) {
  int i, j, r, c, k;

  for (i = 0; i < n; i++) {
    for (j = 0; j < m; j++) {
      a[i][j] = mtr[i][j];
    }
  }

  for (r = c = 0; c < m; c++) {
    j = r;
    for (i = r + 1; i < n; i++) {
      if (fabs(a[i][c]) > fabs(a[j][c])) {
        j = i;
      }
    }
    if (fabs(a[j][c]) < eps) {
      continue;
    }
    for (i = 0; i < m; i++) {
      swap(a[j][i], a[r][i]);
    }

    double s = a[r][c];
    for (j = 0; j < m; j++) {
      a[r][j] /= s;
    }
    for (i = 0; i < n; i++) {
      if (i != r) {
        double t = a[i][c];
        for (j = 0; j < m; j++) {
          a[i][j] -= t * a[r][j];
        }
      }
    }
    r++;
  }
  return r;
}

double GaussJordan(int mtr[MAX_N][MAX_N], int n, int m){
  double det = 1;

  int i, j, k;
  for (i = 0; i < n; i++) {
    for (j = 0; j < n; j++) {
      a[i][j] = mtr[i][j];
    }
  }

  for (k = 0; k < n; k++) {
    j = k;
    for (i = k + 1; i < n; i++) {
      if (fabs(a[i][k]) > fabs(a[j][k])) {
        j = i;
      }
    }

    for (i = 0; i < n; i++) {
      swap(a[j][i],a[k][i]);
    }
    if (j != k) {
      det *= NIL;
    }

    double s = a[k][k];
    for (j = 0; j < n; j++) {
      a[k][j] /= s;
    }

    det *= s;
    for (i = 0; i < n; i++) {
      if (i != k){
        double t = a[i][k];
        for (j = 0; j < n; j++) {
          a[i][j] -= t * a[k][j];
        }
      }
    }
  }
  return det;
}

void build(int n, int col, int task) {
  int i;
  for (i = 0; i < n; i++) {
    move[i][col] = task ? init[i][col] : init[i][0];
  }
}

void process(int n, int line) {
  int i, j;
  for (i = 0; i < n; i++) {
    if (i != line) {
      for (j = 1; j < n; j++) {
        tmp[i - (i > line)][j - 1] = move[i][j];
      }
    }
  }
}

int main(void) {
  int i, j, val;
  FILE *f = fopen("gauss.in", "r");
  freopen("gauss.out", "w", stdout);

  fscanf(f, "%d %d", &M, &N);
  for (i = 0; i < M; i++) {
    for (j = 0; j <= N; j++) {
      fscanf(f, "%d", &val);
      if (j < N) {
        init[i][j] = val;
      }
      expand[i][j] = val;
    }
  }
  fclose(f);

  int rankInit = rank(init, M, N);
  int rankExpand = rank(expand, M, N + 1);

  fprintf(stderr, "%d %d\n", rankInit, rankExpand);

  if (rankInit == rankExpand) {
    if (rankInit > 0) {
      int outofOrder = N - rankInit + 1;
      for (i = 0; i < rankInit; i++) {
        for (j = 0; j < outofOrder; j++) {
          if (j) {
            out[i][j] = -init[i][N - j];
          } else {
            out[i][j] = expand[i][N];
          }
        }
      }

      double down = GaussJordan(init, rankInit, rankInit);

      //fprintf(stderr, "deter = %f\n", down);

      if (down == 0) {
        puts("imposibil");
        //fprintf(stdout, "Se pare ca determinantul sistemului tau este fix 0! Cacealmist!\n");
      } else {
        for (i = 0; i < rankInit; i++) {
          for (j = 0; j < rankInit; j++) {
            move[i][j] = init[i][j];
          }
        }

        int x, sign = 1;
        for (x = 0; x < rankInit; x++) {
          if (x) {
            build(rankInit, x, 0);
            sign = NIL;
          }

          for (i = 0; i < rankInit; i++) {
            process(rankInit, i);
            double plus = (double)sign * even[i & 1] * GaussJordan(tmp, rankInit - 1, rankInit - 1);
            for (j = 0; j < outofOrder; j++) {
              freq[j] += (double)out[i][j] * plus;
            }
          }
          if (x) {
            build(rankInit, x, 1);
          }

          fprintf(stdout, "%.8lf ", freq[0] / down);
          //fprintf(stdout, "x%d  = [", x + 1);
          for (j = 0; j < outofOrder; j++) {
            /*if (j) {
              if (freq[j]) {
                if (freq[j] < 0) {
                  fprintf(stdout, " + %c * (%.0f)", (j - 1) + 'a', freq[j]);
                } else {
                  fprintf(stdout, " + %c * %.0f", (j - 1) + 'a', freq[j]);
                }
              }
            } else {
              fprintf(stdout, "%.0f", freq[j]);
            }*/
            freq[j] = 0;
          }
          //fprintf(stdout, "] / %.0f\n", down);
        }
        for (; x < N; x++) {
          fprintf(stdout, "%.8lf ", 0);
        }
      }
    } else {
      puts("imposibil");
      //fprintf(stdout, "Daca se afiseaza acest mesaj inseamna ca faci cu Onisor...\n");
    }
  } else {
    puts("imposibil");
    //fprintf(stdout, "Ori esti Tavi sau Skior, ori din greseala ai ales un sistem incompatibil!\n");
  }
  fclose(stdout);
  return 0;
}