hello-barcelona-17/day3/g/g.orig.cpp

#include <bits/stdc++.h>

#define FOR(i, n) for(lli i = 0; i < (lli)(n); ++i)

#define X(A) get<0>(A)
#define Y(A) get<1>(A)
#define Z(A) get<2>(A)
#define W(A) get<3>(A)

#define mt make_tuple
#define mp make_pair
#define pb push_back

#define fst first
#define snd second


using namespace std;
using lli = long long int;
using ll = long long int;

using pii   = tuple<lli, lli>;
using piii  = tuple<lli, lli, lli>;
using vi    = vector<lli>;
using vii   = vector<pii>;
using viii  = vector<piii>;
using vvi   = vector<vi>;
using vvii  = vector<vii>;
using vviii = vector<viii>;
using vb    = vector<bool>;
using vvb   = vector<vb>;


template<int N> // N = dimension
struct gauss {
  bitset<N> A[N];     // upper diagonal matrix (free family)
  int V[N]; int nv=0; // identify the columns
  bitset<N> B[N];     // how elements of A are generated ?

  gauss() { reset(); }
  void reset() {
    // memset(this,0,sizeof(gauss<N>)); should work as well
    FOR(i,N) A[i].reset();
    nv=0;
    FOR(i,N) B[i].reset();
  }

  // add a matrix column (a vector in the free family)
  void add(int id, bitset<N> x) {
    if(nv == N) return;
    bitset<N> y; y[nv]=1;
    FOR(i,N) if(x[i]) {
      if(A[i][i]) { x ^= A[i]; y ^= B[i]; }
      else { A[i] = x; V[nv++] = id; B[i] = y; return; }
    }
  }

  // reach a vector x
  // returns true when possible, false otherwise
  // y contains the vectors of V used
  bool solve(bitset<N> x, bitset<N> &y){
    y.reset();
    FOR(i,N) if(x[i]&&A[i][i]) {
      x ^= A[i];
      y ^= B[i];
    }
    if(x.any()) return false;
    return true;
  }
};

// example usage : solve a linear system in Z/2Z
// M is a matrix (list of columns)
// X is a vector
// solves M Y = X
// returns in Y the indices of columns of M with xor X
const int DIM=2500;
bool solve(vector<bitset<DIM> > M, bitset<DIM> X, vi &Y) {
  gauss<DIM> G;
  FOR(i,M.size()) {
    G.add(i, M[i]);
  }
  bitset<DIM> y;
  if (!G.solve(X, y))
      return false;
  for (int i = 0; i < DIM; ++i)
      if (y[i])
        Y.pb(G.V[i]);
  return true;
}


// Inspire de https://github.com/stjepang/snippets/blob/master/gauss.cpp
// Elimination de Gauss
// Resout un systeme d'equations lineaires
// Complexite : O(nb_lins * nb_cols^2)
// Si le systeme a au moins une solution, value contiendra une solution possible
const int MAX_NB_COLS = 2500;
const double eps = 1e-8;

struct Gauss {
  // posnz[i] = -1 si la i-eme composante est libre
  int posnz[MAX_NB_COLS];
  // value[i] = la valeur de X(i) verifiant l'equation ci-dessous 
  int value[MAX_NB_COLS];
  // vrai ssi. le systeme a >= 1 solution
  bool has_solution;

  // mat[0 .. nb_lins-1][0 .. nb_cols-1] * X = mat[0 .. nb_lins-1][nb_cols]
  Gauss(int mat[][MAX_NB_COLS + 1], int nb_lins, int nb_cols) {
    fill(posnz, posnz + nb_cols, -1);


    for (int i = 0; i < nb_lins; ++i, puts(""))
        for (int j = 0; j < nb_cols; ++j)
        printf("%d ", mat[i][j]);
    int posnz_cur = 0; 
    for (int col = 0; col < nb_cols; ++col) {
      int max_lin = posnz_cur;

      for (int lin = max_lin + 1; lin < nb_lins; ++lin)
        if (fabs(mat[lin][col]) > fabs(mat[max_lin][col]))
          max_lin = lin;

      // La colonne est nulle
      // Condition de la forme == 0 si on est dans les entiers
      if (mat[max_lin][col] == 0) continue;

      for (int i = 0; i <= nb_cols; ++i)
        swap(mat[max_lin][i], mat[posnz_cur][i]);

      for (int lin = 0; lin < nb_lins; ++lin) {
        if (lin == posnz_cur) continue;
        // Pour Gauss modulaire : remplacer par l'inverse de mat[posnz_cur][col]
        int factor = mat[lin][col] * (1 ^ mat[posnz_cur][col]);
        for (int i = 0; i <= nb_cols; ++i) {
          mat[lin][i] -= factor * mat[posnz_cur][i];
          mat[lin][i] %= 2;
          mat[lin][i] += 2;
          mat[lin][i] %= 2;
        }
          // Gauss mod : rajouter le modulo
      }

      posnz[col] = posnz_cur++;
    }

    // Genere une solution valide
    for (int col = 0; col < nb_cols; ++col) {
      if (posnz[col] != -1)
        value[col] = mat[posnz[col]][nb_cols] * (1 ^ mat[posnz[col]][col]);
      // Gauss mod
      else 
        value[col] = 0;
    }

    for (int i = 0; i < nb_lins; ++i, puts(""))
        for (int j = 0; j < nb_cols; ++j)
        printf("%d ", mat[i][j]);

    // Verifie que la solution generee est valide
    has_solution = true;
    for (int lin = 0; lin < nb_lins; ++lin) {
      int sum = 0;
      for (int col = 0; col < nb_cols; ++col) {
        sum += mat[lin][col] * value[col]; // Gauss mod
        sum %= 2;
      }
      if (sum - mat[lin][nb_cols] != 0) // Gauss mod
        has_solution = false;
    }
  }
};


const int MAXN = 50, MAXM = 50;
int grid[MAXN][MAXM];
int mat[MAXN * MAXM][MAXN * MAXM + 1];
int N, M;

int linof(int i) { return i / M; }
int colof(int i) { return i % M; }

int main() {
    scanf("%d%d", &N, &M); 
    for (int i = 0; i < N; ++i)
        for (int j = 0; j < M; ++j) {
           char c;
           scanf(" %c", &c);
           grid[i][j] = c == 'W';
        }
    vector<bitset<DIM>> mat(DIM);
    bitset<DIM> X;
    for (int i = 0; i < N * M; ++i)
        for (int j = 0; j < N * M; ++j)
            mat[i][j] = linof(i) == linof(j) || colof(i) == colof(j);
    for (int i = 0; i < N * M; ++i)
        X[i] = grid[linof(i)][colof(i)] ^ 1;
    /*
    for (int i = 0; i < N * M; ++i)
        mat[i][N * M] = grid[linof(i)][colof(i)] ^ 1;
    */
    vi Y;
    if (solve(mat, X, Y)) {
        puts("Yes");
        printf("%d\n", (int)Y.size());
        for (int x: Y)
            printf("%d %d\n", linof(x) + 1, colof(x) + 1);
    } else
        puts("No");
    return 0;
}