//Stuart Anderson stuart.errol.anderson@gmail.com
// program to count unique pairs of NxN matrices A, B that satisfy AxB = A+B; 
// Nov 15 2024 

#include <iostream>
#include <vector>
using namespace std;
unsigned int pow(unsigned int b, unsigned int e) {
    if (e == 0) return 1;
    if (e == 1) return b;
    if (e % 2 == 0) {
        return pow(b, e/2) * pow(b, e/2);
    }
    return b * pow(b, e - 1);
}
void displayMatrix(vector<vector<int>> &A) {
    for (const auto &row : A) {
        for (const auto value : row) {
            cout << value << ' ';
        }
        cout << '\n';
    }
    cout << '\n';  
}
bool equalMatrices(vector<vector<int>> &A, vector<vector<int>> &B, int N) {
  bool equal = true;
  for(int i=0; i<N; i++) {
    for(int j=0; j<N; j++) {
      if(A[i][j] != B[i][j]) {
        equal = false;
        break;
      }
    }
  }
  return equal;
}
void multiplyMatrices( int N,
    vector<vector<int>> &A,
    vector<vector<int>> &B, 
    vector<vector<int>> &C){
    for(int i = 0; i < N; i++){
        for(int j = 0; j < N; j++){
            C[i][j] = 0;
            for(int k = 0; k < N; k++){
                C[i][j] += A[i][k] * B[k][j];
            }
        }
    }
}
void addMatrices( int N,
 vector<vector<int>> &A,
 vector<vector<int>> &B,
 vector<vector<int>> &D) {
  for(int i=0; i<N; i++) {
    for(int j=0; j<N; j++) {
      D[i][j] = A[i][j] + B[i][j];
    }
  }
}
void populateMatrix( int N, int c, vector<vector<int>> &A) {
	for (int i = 0; i < N && c != 0; i++) {
	   for (int j = 0; j < N && c != 0; j++) {
		A[i][j] = c % 2;
		c /= 2; 
	    }        
        }  
    }
void zeroMatrix( int N, vector<vector<int>> &B) {
  for(int i=0; i<N; i++) {
    for(int j=0; j<N; j++) {
      B[i][j] = 0;
    }
  }
}	
void generateMatrices(int N, vector<vector<int>> &A , vector<vector<int>> &B, vector<vector<int>> &C , vector<vector<int>> &D) {
    unsigned int max = pow(2, N*N);
    int z = 0;
    int y = 0;
    for (unsigned int counter = 0; counter < max; ++counter) {
        int c = counter;
        zeroMatrix( N, B) ;
         for (unsigned int dcounter = 0; dcounter <= counter; ++dcounter) {
           int d = dcounter;
           populateMatrix( N, c, A);
           populateMatrix( N, d, B); 
           addMatrices(N, A, B, C);
           multiplyMatrices(N, A, B, D);
           if (equalMatrices(C,D,N)==true) {
               displayMatrix(A);
	       displayMatrix(B);
               cout<<"A "<<counter<<", B "<<dcounter<<endl;
               z++;
	       if (equalMatrices(A,B,N)==true) {
	          y++;
	          }
            } 
        }  
    }
    cout<<z<<" unique solutions for "<<N<<"x"<<N<<" where A+B = A*B"<<endl;
    cout<<y<<" solutions where A == B"<<endl;
}

int main() {
int size;
cout<< "Program to generate and count unique pairs of NxN matrices A, B that satisfy AxB = A+B "<<endl;
cout<<"input matrix dimension N: ";
cin >> size;
vector<vector<int>> A(size, vector<int>(size));
vector<vector<int>> B(size, vector<int>(size));
vector<vector<int>> C(size, vector<int>(size));
vector<vector<int>> D(size, vector<int>(size));
    generateMatrices(size, A, B, C, D);
    return 0;
}