function [ A, b ]  = splitting_monohedral_matrix ( r_shape, p_shape, C )

%*****************************************************************************80
%
%% splitting_monohedral_matrix() constructs a monohedral tiling linear system.
%
%  Discusssion:
%
%    This function constructs the linear system Ax = b for a pariomino 
%    splitting problem associated with a monohedral polyomino tiling problem.
%
%    In a monohedral tiling problem, the tiling is to be carried out using
%    copies of a single polyomino shape.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    16 February 2023
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer r_shape(rowsR,colsR): a rectangular matrix whose entries are 0 or 1.
%    The 1 values define the shape of the tiling region.
%
%    integer p_shape(rows,cols): a rectangular matrix whose entries are 0 or 1.
%    The 1 values define the shape of the polyomino used for tiling.
%
%    integer C: the number of colors.
%
%  Output:
%  
%    integer A(C,C): the matrix of constraint coefficients.
%
%    integer b(C): the right hand sides of the constraint equations.
%

%
%  Allocate A and b.
%
  n = C;
  m = C; 
  b = zeros ( m, 1 );
  A = zeros ( m, n ); 
%
%  Construct the RHS vector b first. 
%
%  If r_shape does NOT represent a rectangle, we use the first of 
%  the C variants associated with the r_shape.
%  If r_shape DOES represent a rectangle, we use the variant of r_shape 
%  that has a '1' in the (1,1) position to be consistent
%  with a lemma in the paper.
%
  [ rowsR, colsR ] = size ( r_shape );
  
  if ( any ( r_shape(:)==0 ) )
    for j = 1 : colsR 
      for i = 1 : rowsR
        if ( r_shape(i,j) == 1 )
          r_shape_pari(i,j) = mod ( i + j, C ) + 1;
        end
      end
    end
  else
    for page = 1 : C
      variantsR(:,:,page) = zeros ( rowsR, colsR);
    end
    for j = 1:colsR
      for i = 1:rowsR
        if r_shape(i,j) == 1
          for page = 1 : C
            variantsR(i,j,page) = mod(i+j+page-1,C) + 1;
          end
        end
      end
    end
%
%  Find which variant of r_shape has a '1' in the (1,1) position.\
%
    for page = 1 : C
      ind = variantsR(1,1,page);
      if ( ind == 1 )
        keep = page;
        break
      end
    end
    r_shape_pari = variantsR(:,:,keep);
  end
  
  for colour = 1 : C
    b(colour) = length ( find ( r_shape_pari == colour ) );
  end   
%
%  Construct the coefficient matrix A.  Initialize the array that stores 
%  all n variants of the pariominoes.  Each variant is packed into a
%  matrix with the same size as r_shape.
%
  for page = 1 : C
    variants(:,:,page) = zeros ( rowsR, colsR );
  end

  polyomino = p_shape;
  [ rows, cols ] = size ( polyomino );
  for j = 1 : cols
    for i = 1 : rows
      if ( polyomino(i,j) == 1 )
        for colour = 1 : C
          variants(i,j,colour) = mod ( i+j+colour-1, C ) + 1;
        end
      end
    end
  end
%
%  Count how many 1's, 2's, ..., C's there are in each variant and
%  construct the matrix A.
%
  for colour = 1 : C
    for page = 1 : C
      polyomino = variants(:,:,page);
      A(colour,page) = length ( find ( polyomino == colour ) );
    end
  end
    
  return
end
          
