function reid_tiling ( )

%*****************************************************************************80
%
%% reid_tiling() computes tilings of the Reid region using RREF.
%
%  Discussion:
%
%    A region of 8 squares is to be covered by 4 dominos.
%
%    The requirement that each square be covered becomes 8 equations,
%    that there must be exactly one domino square covering each 
%    region square.
%
%    The right hand side b is of length 8, with each entry being a
%    1, reflecting the number of domino squares covering that region square.
%
%    There are exactly 10 distinct ways that a domino can be placed
%    on the region.  We regard each of these ways as a variable v(j).
%    
%    Our solution x(k) is a 1 if variable v(k) is used in the tiling
%    and 0 otherwise.
%
%    This can be formulated as an 8 x 10 underdetermined linear system
%    of the form A*x=b.  The RREF can be used to find the degrees of
%    freedom in the system.  We are free to set these x variables to 
%    0 or 1.  By considering every possible such combination, we uncover
%    four distinct legal solutions and the corresponding tiling patterns.
%
%  Licensing:
%
%    This code is distributed under the MIT license. 
%
%  Modified:
%
%    25 June 2025
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'reid_tiling():\n' );
  fprintf ( 1, '  reid_tiling() uses the RREF to uncover tilings of the Reid polygon.' );
  fprintf ( 1, '' )
%
%  Step 1: Define Region:
%
%   1 1 0
%   1 1 1
%   1 1 1
%
%  These squares are indexed as
%
%    0 1 
%    2 3 4
%    5 6 7
%
  r_shape = ones ( 3, 3 );
  r_shape(1,3) = 0;
  
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Shape of the Reid polygon:\n' );
  fprintf ( 1, '\n' );
  disp ( r_shape );
%
%  Step 2: Define Tile Shape:
%
  p_num = 1;
  p_shapes = zeros ( 3, 3, p_num );

  p_shapes(1:2,1,1) = 1;
%
%  Step 3: Report Tile Shapes:
%
%  1
%  1
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Tile 1:\n' );
  fprintf ( 1, '\n' );
  disp ( p_shapes(1:2,1,1) );
%
%  Define linear system A*x=b.
%
%  The rows are equations, each associated with a square.
%  Equation I specifies that square I is to be covered exactly 1 time.
%
%  Thus the right hand side is (1,1,1,1,1,1,1).
%
%  Column J indicates which of the squares 0 through 7 are covered
%  by tile J.
%
%  This linear system can be generated (the hard way) by careful
%  consideration of all possible positions and orientations of each tile.
%  Each column of this matrix represents such a position and orientation,
%  and the 1 values indicate the squares covered in that case.  Thus,
%  the first column says that this domino position covers squares 0 and 2.
%
  A = [ ...
    1,0,0,0,0,1,0,0,0,0; ...
    0,0,1,0,0,1,0,0,0,0; ...
    1,1,0,0,0,0,1,0,0,0; ...
    0,0,1,1,0,0,1,1,0,0; ...
    0,0,0,0,1,0,0,1,0,0; ...
    0,1,0,0,0,0,0,0,1,0; ...
    0,0,0,1,0,0,0,0,1,1; ...
    0,0,0,0,1,0,0,0,0,1];

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The matrix A:\n' );
  fprintf ( 1, '\n' );
%
%  We need right hand side to be a column vector, hence double brackets.
%
  b = [ 1; 1; 1; 1; 1; 1; 1; 1 ];
%
%  Append b to A.
%
  Ab = [ A, b ];
%
%  Compute the reduced row echelon form of [ A | b ].
%
  RREF = rref_compute ( Ab );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The RREF of [A|b]\n' );
  fprintf ( 1, '\n' );
  disp ( RREF );
%
%  We need a new matrix C that combines the RREF information with
%  dummy equations for the free variables.
%
  m = size ( RREF, 1 );
  np1 = size ( RREF, 2 );
  n = np1 - 1;

  C = eye ( n );
  b_RREF = zeros ( n, 1 );
  f = 1 : n;
%
%  Active rows of RREF replace corresponding rows of C.
%  f records the "free" variables.
%
  for i = 1 : m
    for j = 1 : n
      if ( RREF(i,j) == 1.0 )
        C(j,:) = RREF(i,1:n);
        b_RREF(j) = RREF(i,end);
        k = find ( f == j );
        if ( k )
          f = [ f(1:k-1),f(k+1:end) ];
        end
        break
      end
    end
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The new linear system C:\n' );
  fprintf ( 1, '\n' );
  disp ( C );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The free variable indices f:\n' );
  fprintf ( 1, '\n' );
  disp ( f );

  nf = length ( f );
%
%  For nz = 0 through 4 (the number of tiles)
%    Set nz free variables to 1
%      Solve linear system
%        If all variables are 0 or 1, accept it.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Legal solutions:\n' );
  fprintf ( 1, '\n' );

  for nz = 0 : nf

    done = true;
    combo = [];

    while ( true )

      [ combo, done ] = comb_next ( nf, nz, combo, done );

      if ( done )
        break;
      end

      bb = b_RREF;
      bb(f(combo)) = 1;

      x = C \ bb;

      if ( all ( ( x == 0 ) | ( x == 1 ) ) )
        disp ( x' );
      end

    end

  end

  return
end
function [ a, done ] = comb_next ( n, k, a, done )

%*****************************************************************************80
%
%% comb_next() computes combinations of K things out of N.
%
%  Discussion:
%
%    The combinations are computed one at a time, in lexicographical order.
%
%    10 April 1009: Thanks to "edA-qa mort-ora-y" for supplying a
%    correction to this code.
%
%  Example:
%
%  Combinations of size 3 from a group of 5.
%
%     1     2     3
%     1     2     4
%     1     2     5
%     1     3     4
%     1     3     5
%     1     4     5
%     2     3     4
%     2     3     5
%     2     4     5
%     3     4     5
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    10 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Charles Mifsud,
%    Combination in Lexicographic Order,
%    ACM algorithm 154,
%    Communications of the ACM,
%    March 1963.
%
%  Input:
%
%    integer N, the total number of things.
%
%    integer K, the number of things in each combination.
%
%    integer A(K), the output value of A on the previous call.
%    This value is not needed on a startup call.
%
%    logical DONE, should be set to true on the first call,
%    and thereafter set to the output value of DONE on the previous call.
%
%  Output:
%
%    integer(K), the next combination.
%
%    logical DONE, is false if the routine can be called
%    again for more combinations, and true if there are no more.
%
  if ( k < 0 )
    a = [];
    done = true;
    return
  end

  if ( k == 0 )
    if ( done )
      a = [];
      done = false;
    else
      done = true;
    end
    return
  end
%
%  Initialization for 0 < k.
%
  if ( done )

    a = 1 : k;
 
    done = false;
%
%  Next combination.
%
  else

    if ( a(k) < n )
      a(k) = a(k) + 1;
      return
    end

    for i = k : -1 : 2

      if ( a(i-1) < n-k+i-1 )

        a(i-1) = a(i-1) + 1;

        for j = i : k
          a(j) = a(i-1) + j - ( i-1 );
        end

        return

      end

    end

    done = true;

  end

  return
end
function A = rref_compute ( A )

%*****************************************************************************80
%
%% rref_compute() computes the reduced row echelon form of a matrix.
%
%  Discussion:
%
%    A rectangular matrix is in row reduced echelon form if:
%
%    * The leading nonzero entry in each row has the value 1.
%
%    * All entries are zero above and below the leading nonzero entry 
%      in each row.
%
%    * The leading nonzero in each row occurs in a column to
%      the right of the leading nonzero in the previous row.
%
%    * Rows which are entirely zero occur last.
%
%  Example:
%
%    M = 4, N = 7
%
%    Matrix A:
%
%     1    3    0    2    6    3    1
%    -2   -6    0   -2   -8    3    1
%     3    9    0    0    6    6    2
%    -1   -3    0    1    0    9    3
%
%    RREF(A):
%
%     1    3    0    0    2    0    0
%     0    0    0    1    2    0    0
%     0    0    0    0    0    1   1/3
%     0    0    0    0    0    0    0
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    27 December 2024
%
%  Author:
%
%    Original Python version by ChatGPT.
%    This version by John Burkardt.
%
%  Reference:
%
%    Charles Cullen,
%    An Introduction to Numerical Linear Algebra,
%    PWS Publishing Company, 1994,
%    ISBN: 978-0534936903,
%    LC: QA185.D37.C85.
%
%  Input:
%
%    real A(M,N), the matrix to be analyzed. 
%
%  Output:
%
%    real A(M,N), the reduced row echelon form of the matrix.
%
%    integer a_cols(*): the columns for which a pivot entry was found.
%

%
%  Initialize an empty list of independent columns.
%
  a_cols = [];
%
%  Get the dimensions of the matrix.
%
  [ rows, cols ] = size ( A );
%
%  Set a tolerance for the size of a pivot value.
%
  tol = sqrt ( eps );
%
%  Start from the first row.
%
  row = 1;
%
%  Seek a pivot for each column.
%
  for col = 1 : cols
%
%  Exit if we have run out of rows to examine.
%
    if ( rows < row )
      break;
    end
%
%  Find row index P of maximum element in subvector A(R:M,C).
%
    [ pivot_value, pivot_row ] = max ( abs ( A(row:rows,col) ) );
%
%  Adjust value of P as an entry in full vector A(1:M,C).
%
    pivot_row = pivot_row + row - 1;
%
%  If A(row:rows,col) was essentially zero, there is no pivot for this column.
%
    if ( abs ( pivot_value ) <= tol )
      A(pivot_row,col) = 0.0;
      continue;
    end
%
%  A pivot was found for column col.
%
    a_cols = [ a_cols, col ];
%
%  Swap current row and pivot row
%
    temp = A(row,col:cols);
    A(row,col:cols) = A(pivot_row,col:cols);
    A(pivot_row,col:cols) = temp;
%
%  Normalize row R so that A(row,col) = 1.
%
    A(row,col:cols) = A(row,col:cols) / A(row,col);
%
%  Eliminate nonzeros in column C using multiples of row R, except for A(R,C).
%
    for i = 1 : rows
      if ( i ~= row )
        A(i,col:cols) = A(i,col:cols) - A(i,col) * A(row,col:cols);
      end
    end
%
%  Move to the next row.
%
    row = row + 1;

  end

  return
end

