#include <vector>
#include <array>
#include <numeric>
#include <iostream>

constexpr size_t M = 3;
constexpr size_t N = 3;
constexpr size_t MN = M * N;
using Mat = std::array<std::array<bool,MN>,MN>;
using Vec = std::array<bool,MN>;

struct point {
    point():x{0},y{0}{}
    point(size_t x_, size_t y_):x{x_},y{y_}{}

    size_t x,y;

    point& inc() {
        point& ret = *this;

        ++y;
        if(y >= N){
            ++x;
            y -= N;
        }

        return *this;
    }

    size_t stride() const {
        return x + y*N;
    }

    bool equal(const point& rhs) const {
        return (x == rhs.x && y == rhs.y);
    }
};

void print_vec(const Vec& x){
    for(size_t j = 0; j < N; ++j){
        for(size_t i = 0; i < M; ++i){
            point ind{i,j};
            std::cout<<x[ind.stride()]<<" ";
        }
        std::cout<<"\n";
    }
    std::cout<<std::endl;
}

void print_perm(const std::array<size_t,MN>& x){
    for(size_t j = 0; j < MN; ++j){
        std::cout<<x[j]<<"\n";
    }
    std::cout<<std::endl;
}

void print_mat(const Mat& x, const std::array<size_t,MN>& P){
    for(size_t j = 0; j < MN; ++j){
        for(size_t i = 0; i < MN; ++i){
            std::cout<<x[P[j]][i]<<" ";
        }
        std::cout<<"\n";
    }
    std::cout<<std::endl; 
}

bool lu_decompose_gf2(Mat& R, std::array<size_t,MN>& P){
    for(size_t i = 0; i < MN; ++i){
        bool found_pivot = false;
        size_t pivot_row = i;
        for(size_t k = i; k < MN && (not found_pivot); ++k){
            if(R[P[k]][i]){
                found_pivot = true;
                pivot_row = k;
            }
        }
        if(not found_pivot){
            return false;
        }else if ( i != pivot_row ){
            std::swap(P[i],P[pivot_row]);
        }
        for(size_t k = i+1u; k < MN; ++k){

            if(R[P[k]][i]){
                for(size_t j = i+1u; j < MN; ++j){
                    R[P[k]][j] ^= (R[P[i]][j]);
                }
            }
        }
    }

    return true;
}

int main(){
    Mat A;
    
    Vec b{
        1,0,1,
        1,0,1,
        0,0,0
    };

    // Works best if we use odd N's
    // even N might not be solvable in roughly half of the cases
    // Init
    std::cout<<"Init"<<std::endl;
    for(point i{}; i.x < M && i.y < N; i.inc()){
        for(point j{}; j.x < M && j.y < N; j.inc()){
            if(
                (i.equal(j)) ||
                (i.x == j.x && (i.y+1) == j.y) ||
                (i.x == j.x && (j.y+1) == i.y) ||
                (i.y == j.y && (i.x+1) == j.x) ||
                (i.y == j.y && (j.x+1) == i.x)
            ){
                A[i.stride()][j.stride()] = true;
            }else{
                A[i.stride()][j.stride()] = false;
            }
            // std::cout<<A[i.stride()][j.stride()]<<std::endl;
        }
        // Populate how you like it.
        //b[i.stride()] = false;
    }
    std::array<size_t, MN> P_ident;
    std::iota(P_ident.begin(),P_ident.end(),0);
    print_mat(A,P_ident);

    std::array<size_t, MN> P;
    std::iota(P.begin(),P.end(),0);
    // Solve it :)
    std::cout<<"Solve"<<std::endl;
    if(!lu_decompose_gf2(A,P)){
        return -1;
    }
    std::array<size_t,MN> P_inv{};
    for(size_t i = 0; i < MN; ++i){
        P_inv[P[i]] = i;
    }

    std::cout<<"LR"<<std::endl;
    print_mat(A,P);

    std::cout<<"P"<<std::endl;
    print_perm(P);

    std::cout<<"Lights-Off puzzle"<<std::endl;
    print_vec(b);

    Vec y;
    for(size_t i = 0; i < MN; ++i){
        y[i] = b[P[i]];
        for(size_t j = 0; j < i; ++j){
            y[i] ^= A[P[i]][j] & y[j];
        }
    }

    Vec x;
    for(size_t i = MN; i > 0; --i){
        x[i-1] = y[i-1];
        for(size_t j = i; j < MN; ++j){
            x[i-1] ^= (A[P[i-1]][j] and x[j]);
        }
    }

    std::cout<<"Solution"<<std::endl;
    print_vec(x);

    return 0;
}