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

constexpr size_t M = 3;
constexpr size_t N = 3;

template<size_t MM, size_t NN>
using Mat = std::array<std::array<int8_t,MM*NN>,MM*NN>;

template<size_t MM, size_t NN>
using Vec = std::array<int8_t,MM*NN>;

template<size_t MM, size_t NN>
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 >= NN){
            ++x;
            y -= NN;
        }

        return *this;
    }

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

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

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

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

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

template<size_t MM, size_t NN>
bool lu_decompose_gf2(Mat<MM,NN>& R, std::array<size_t,MM*NN>& P){
    for(size_t i = 0; i < MM*NN; ++i){
        bool found_pivot = false;
        size_t pivot_row = i;
        for(size_t k = i; k < MM*NN && (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 < MM*NN; ++k){

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

    return true;
}

template<uint64_t MM, uint64_t NN>
int solve_lights_out(std::array<int8_t,MM*NN>& b_raw){
    Mat<M,N> A;
		constexpr size_t MN = M*N;
    
    Vec<M,N> 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<M,N> i{}; i.x < M && i.y < N; i.inc()){
        for(point<M,N> 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<M,N>(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<M,N>(A,P)){
        return -1;
    }
    std::array<size_t,MN> P_inv{};
    for(size_t i = 0; i < MN; ++i){
        P_inv[P[i]] = i;
    }

    Vec<M,N> 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<M,N> 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]);
        }
    }

    return 0;
}

extern "C"{
int solve_lights_out_3_x_3(std::array<int8_t, 3u*3u>& inp){
	int rc = solve_lights_out<3u,3u>(inp);
	return rc;
}

int fake_main(){
    Mat<M,N> A;
		constexpr size_t MN = M*N;
    
    Vec<M,N> 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<M,N> i{}; i.x < M && i.y < N; i.inc()){
        for(point<M,N> 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<M,N>(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<M,N>(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<M,N>(A,P);

    std::cout<<"P"<<std::endl;
    print_perm<M,N>(P);

    std::cout<<"Lights-Off puzzle"<<std::endl;
    print_vec<M,N>(b);

    Vec<M,N> 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<M,N> 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<M,N>(x);

    return 0;
}
}

int main(){
	return fake_main();
}