Warning: If you are looking for pizza delivery you have come to the wrong page. Bloody Google (other search engines are available).

Task
Dominoes
You are encouraged to solve this task according to the task description, using any language you may know.

Take a box of dominoes give them a good shuffle and then arrange them in diverse orientations such that they form a rectangle with 7 rows of 8 columns. Make a tableau of the face values e.g.

       05132231
       05505246
       43036620
       06235126
       11300245
       21433466
       64515414

Now torment your computer by making it identify where each domino is.

Do this for the above tableau and one of your own construction.

Extra credit:

How many ways are there to arrange dominoes in an 8x7 rectangle, first ignoring their values, then considering their values, and finally considering values but ignoring value symmetry, i.e. transposing 5 and 4.

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

const int32_t EMPTY = -1;

const std::vector<std::vector<int32_t>> tableau_one = { { 0, 5, 1, 3, 2, 2, 3, 1 },
												        { 0, 5, 5, 0, 5, 2, 4, 6 },
												        { 4, 3, 0, 3, 6, 6, 2, 0 },
	                                                    { 0, 6, 2, 3, 5, 1, 2, 6 },
	                                                    { 1, 1, 3, 0, 0, 2, 4, 5 },
	                                                    { 2, 1, 4, 3, 3, 4, 6, 6 },
	                                                    { 6, 4, 5, 1, 5, 4, 1, 4 } };

const std::vector<std::vector<int32_t>> tableau_two = { { 6, 4, 2, 2, 0, 6, 5, 0 },
	                                                    { 1, 6, 2, 3, 4, 1, 4, 3 },
														{ 2, 1, 0, 2, 3, 5, 5, 1 },
														{ 1, 3, 5, 0, 5, 6, 1, 0 },
														{ 4, 2, 6, 0, 4, 0, 1, 1 },
														{ 4, 4, 2, 0, 5, 3, 6, 3 },
														{ 6, 6, 5, 2, 5, 3, 3, 4 } };

class Domino {
public:
	Domino(const int32_t& aOne, const int32_t& aTwo) {
		one = std::min(aOne, aTwo);
		two = std::max(aOne, aTwo);
	}

	bool operator==(const Domino& other) const {
		return one == other.one && two == other.two;
	}

private:
	int32_t one, two;

};

class Point {
public:
	Point(const int32_t& aX, const int32_t& aY) : x(aX), y(aY) { }

	bool operator==(const Point& other) const {
		return x == other.x && y == other.y;
	}

	int32_t x, y;
};

struct Pattern {
	std::vector<std::vector<int32_t>> tableau;
	std::vector<Domino> dominoes;
	std::vector<Point> points;
};

int32_t first_empty_cell(const std::vector<std::vector<int32_t>>& tableau) {
	for ( uint64_t row = 0; row < tableau.size(); ++row ) {
		for ( uint64_t col = 0; col < tableau[0].size(); ++col ) {
			if ( tableau[row][col] == EMPTY ) {
				return row * tableau[0].size() + col;
			}
		}
	}
	return EMPTY;
}

void print_layout(const Pattern& pattern) {
	std::vector<std::string> output(
        2 * pattern.tableau.size(), std::string(2 * pattern.tableau[0].size() - 1, ' '));

	for ( uint64_t i = 0; i < pattern.points.size() - 1; i += 2 ) {
		const int x1 = pattern.points[i].x;
		const int y1 = pattern.points[i].y;
		const int x2 = pattern.points[i + 1].x;
		const int y2 = pattern.points[i + 1].y;
		const int n1 = pattern.tableau[x1][y1];
		const int n2 = pattern.tableau[x2][y2];
		output[2 * x1][2 * y1] = static_cast<char>('0' + n1);
		output[2 * x2][2 * y2] = static_cast<char>('0' + n2);
		if ( x1 == x2 ) {
			output[2 * x1][2 * y1 + 1] = '+';
		} else if ( y1 == y2 ) {
			output[2 * x1 + 1][2 * y1] = '+';
		}
	}

	for ( const std::string& line : output ) {
		std::cout << line << "\n";
	}
}

std::vector<Pattern> find_patterns(const std::vector<std::vector<int32_t>>& tableau) {
	const int32_t n_rows = tableau.size();
	const int32_t n_cols = tableau[0].size();
	const uint64_t domino_count = ( n_rows * n_cols ) / 2;
	std::vector<std::vector<int32_t>> empty_tableau =
        { static_cast<uint64_t>(n_rows), std::vector<int32_t>(n_cols, EMPTY) };
	std::vector<Pattern> patterns = { Pattern(empty_tableau, std::vector<Domino>(), std::vector<Point>()) };

	while ( true ) {
		std::vector<Pattern> next_patterns = { };
		for ( Pattern pattern : patterns ) {
			std::vector<std::vector<int32_t>> next_tableau = pattern.tableau;
			std::vector<Domino> dominoes = pattern.dominoes;
			std::vector<Point> points = pattern.points;

			int32_t index = first_empty_cell(next_tableau);
			if ( index == EMPTY ) {
				continue;
			}

			int32_t row = index / n_cols;
			int32_t col = index % n_cols;
			if ( row + 1 < n_rows && next_tableau[row + 1][col] == EMPTY ) {
				Domino domino(tableau[row][col], tableau[row + 1][col]);
				if ( std::find(dominoes.begin(), dominoes.end(), domino) == dominoes.end() ) {
					std::vector<std::vector<int32_t>> final_tableau = { };
					std::copy(next_tableau.begin(), next_tableau.end(), std::back_inserter(final_tableau));
					final_tableau[row][col] = tableau[row][col];
					final_tableau[row + 1][col] = tableau[row + 1][col];
					std::vector<Domino> next_dominoes(dominoes.begin(), dominoes.end());
					next_dominoes.emplace_back(domino);
					std::vector<Point> next_points(points.begin(), points.end());
					next_points.emplace_back(Point(row, col));
					next_points.emplace_back(Point(row + 1, col));
					next_patterns.emplace_back(Pattern(final_tableau, next_dominoes, next_points));
				}
			}

			if ( col + 1 < n_cols && next_tableau[row][col + 1] == EMPTY ) {
				Domino domino(tableau[row][col], tableau[row][col + 1]);
				if ( std::find(dominoes.begin(), dominoes.end(), domino) == dominoes.end() ) {
					next_tableau[row][col] = tableau[row][col];
					next_tableau[row][col + 1] = tableau[row][col + 1];
					dominoes.emplace_back(domino);
					points.emplace_back(Point(row, col));
					points.emplace_back(Point(row, col + 1));
					next_patterns.emplace_back(Pattern(next_tableau, dominoes, points));
				}
			}
		}

		if ( next_patterns.empty() ) {
			break;
		}
		patterns = next_patterns;
		if ( patterns[0].dominoes.size() == domino_count ) {
			break;
		}
	}

	return patterns;
}

int main() {
	for ( const std::vector<std::vector<int32_t>>& tableau : { tableau_one, tableau_two } ) {
		std::vector<Pattern> patterns = find_patterns(tableau);
		std::cout << "Layouts found: " << patterns.size() << "\n";
		print_layout(patterns[0]);
	}
}
Output:
Layouts found: 1
0+5 1+3 2 2+3 1
        +     +
0 5+5 0 5 2+4 6
+     +        
4 3 0 3 6+6 2 0
  + +       + +
0 6 2 3+5 1 2 6
+         +    
1 1 3 0+0 2 4 5
  + +       + +
2 1 4 3+3 4 6 6
+         +    
6 4+5 1+5 4 1+4
               
Layouts found: 2025
6 4 2 2 0 6+5 0
+ + + + +     +
1 6 2 3 4 1+4 3
               
2 1 0 2 3+5 5 1
+ + + +     + +
1 3 5 0 5 6 1 0
        + +    
4 2 6 0 4 0 1+1
+ + + +        
4 4 2 0 5 3 6 3
        + + + +
6+6 5+2 5 3 3 4        
// Dominoes: Nigel Galloway. November 17th., 2021.
let cP (n:seq<uint64 list * uint64>) g=seq{for y,n in n do for g in g do let l=n^^^g in if n+g=l then yield (g::y,l)}
let rec fG n g=match g with h::t->fG(cP n h)t |_->fst(Seq.head n)
let solve(N:int[])=let fG=let y=fG [([],0UL)]([for g in 0..47->((N.[g],N.[g+8]),(1UL<<<g)+(1UL<<<g+8))]@[for n in 0..6 do for g in n*8..n*8+6->((N.[g],N.[g+1]),(1UL<<<g)+(1UL<<<g+1))]
                                |>List.groupBy(fun((n,g),_)->(min n g,max n g))|>List.sort|>List.map(fun(_,n)->n|>List.map(fun(n,g)->g))) in (fun n g->if List.contains((1UL<<<n)+(1UL<<<g)) y then "+" else " ")
                   N|>Array.chunkBySize 8|>Array.iteri(fun n g->let n=n*8 in [0..6]|>List.iter(fun y->printf $"%d{g.[y]}%s{fG(n+y)(n+y+1)}"); printfn $"%d{g.[7]}"; [0..7]|>List.iter(fun g->printf $"%s{fG(n+g)(n+g+8)} "); printfn "")

solve [|0;5;1;3;2;2;3;1;
        0;5;5;0;5;2;4;6;
        4;3;0;3;6;6;2;0;
        0;6;2;3;5;1;2;6;
        1;1;3;0;0;2;4;5;
        2;1;4;3;3;4;6;6;
        6;4;5;1;5;4;1;4|]
Output:
0+5 1+3 2 2+3 1
        +     +
0 5+5 0 5 2+4 6
+     +
4 3 0 3 6+6 2 0
  + +       + +
0 6 2 3+5 1 2 6
+         +
1 1 3 0+0 2 4 5
  + +       + +
2 1 4 3+3 4 6 6
+         +
6 4+5 1+5 4 1+4
solve [|5;6;2;0;0;4;1;4;
        3;6;1;3;0;4;2;2;
        3;5;6;4;3;2;1;1;
        3;5;1;1;3;0;0;5;
        6;0;5;4;3;5;5;2;
        4;4;1;3;6;6;0;2;
        1;2;6;2;6;5;0;4|]
Output:
5+6 2+0 0 4 1+4
        + +
3+6 1 3 0 4 2+2
    + +
3 5 6 4 3+2 1 1
+ +         + +
3 5 1+1 3+0 0 5

6 0 5+4 3+5 5+2
+ +
4 4 1+3 6 6+0 2
        +     +
1+2 6+2 6 5+0 4
Translation of: Julia

Basic task

Const NROWS = 7
Const NCOLS = 8
Const MAXDOMINOS = 28
Const MAXSOLUTIONS = 10000

Type Domino
    a As Integer
    b As Integer
End Type

Type Move
    r1 As Integer
    c1 As Integer
    r2 As Integer
    c2 As Integer
End Type

Type Solution
    moves(0 To MAXDOMINOS-1) As Move
End Type

' Initialize tabla data
Dim Shared As Integer tableau1(NROWS-1, NCOLS-1) = { _
{0, 5, 1, 3, 2, 2, 3, 1}, _
{0, 5, 5, 0, 5, 2, 4, 6}, _
{4, 3, 0, 3, 6, 6, 2, 0}, _
{0, 6, 2, 3, 5, 1, 2, 6}, _
{1, 1, 3, 0, 0, 2, 4, 5}, _
{2, 1, 4, 3, 3, 4, 6, 6}, _
{6, 4, 5, 1, 5, 4, 1, 4} }

Dim Shared As Integer tableau2(NROWS-1, NCOLS-1) = { _
{6, 4, 2, 2, 0, 6, 5, 0}, _
{1, 6, 2, 3, 4, 1, 4, 3}, _
{2, 1, 0, 2, 3, 5, 5, 1}, _
{1, 3, 5, 0, 5, 6, 1, 0}, _
{4, 2, 6, 0, 4, 0, 1, 1}, _
{4, 4, 2, 0, 5, 3, 6, 3}, _
{6, 6, 5, 2, 5, 3, 3, 4} }

Function makeDomino(a As Integer, b As Integer) As Domino
    If a < b Then Return Type<Domino>(a, b) Else Return Type<Domino>(b, a)
End Function

Function dominoUsed(d As Domino, used() As Domino, count As Integer) As Boolean
    For i As Integer = 0 To count-1
        If used(i).a = d.a And used(i).b = d.b Then Return True
    Next
    Return False
End Function

Sub findFirstEmpty(grid() As Integer, Byref row As Integer, Byref col As Integer)
    For i As Integer = 0 To NROWS-1
        For j As Integer = 0 To NCOLS-1
            If grid(i*NCOLS+j) = -1 Then
                row = i : col = j : Exit Sub
            End If
        Next
    Next
    row = -1 : col = -1
End Sub

Sub printDominoLayout(moves() As Move, tableau() As Integer)
    Dim As Integer i, j, k
    Dim As String cell(NROWS-1, NCOLS-1), vplus(NROWS-2, NCOLS-1)
    
    For i = 0 To NROWS-1
        For j = 0 To NCOLS-1
            cell(i, j) = Str(tableau(i, j))
        Next
    Next
    For i = 0 To NROWS-2
        For j = 0 To NCOLS-1
            vplus(i, j) = " "
        Next
    Next
    For k = 0 To MAXDOMINOS-1
        Dim m As Move = moves(k)
        If m.r1 = -1 Then Exit For
        If m.r1 = m.r2 Then
            cell(m.r1, m.c1) = cell(m.r1, m.c1) + "+" + cell(m.r2, m.c2)
            cell(m.r2, m.c2) = ""
        Else
            vplus(m.r1, m.c1) = "+"
        End If
    Next
    
    ' Calculate the width of each cell (maximum 3: "a+b")
    Dim As Integer anchoCelda = 3
    
    For i = 0 To NROWS-1
        Dim As String linea = ""
        For j = 0 To NCOLS-1
            Dim As String s = Iif(cell(i, j) <> "", cell(i, j), "")
            ' Pad right for fixed width
            linea &= Left(s + Space(anchoCelda), anchoCelda)
            If j < NCOLS-1 Then linea &= " "
        Next
        linea = Rtrim(linea)
        Print linea
        
        If i < NROWS-1 Then
            Dim As String vline = ""
            For j = 0 To NCOLS-1
                vline &= vplus(i, j)
                ' Pad right for fixed width
                vline &= Space(anchoCelda - 1)
                If j < NCOLS-1 Then vline &= " "
            Next
            vline = Rtrim(vline)
            Print vline
        End If
    Next
    Print
End Sub

Sub solveRec(tableau() As Integer, grid() As Integer, used() As Domino, usedCount As Integer, moves() As Move, moveCount As Integer, solutions() As Solution, Byref num_solutions As Integer)
    Dim As Integer row, col
    findFirstEmpty grid(), row, col
    If row = -1 Then
        If num_solutions < MAXSOLUTIONS Then
            For i As Integer = 0 To MAXDOMINOS-1
                If i < moveCount Then
                    solutions(num_solutions).moves(i) = moves(i)
                Else
                    solutions(num_solutions).moves(i).r1 = -1
                End If
            Next
            num_solutions += 1
        End If
        Exit Sub
    End If
    
    ' Horizontal
    If col < NCOLS-1 Then
        If grid(row*NCOLS+col+1) = -1 Then
            Dim d As Domino = makeDomino(tableau(row, col), tableau(row, col+1))
            If dominoUsed(d, used(), usedCount) = 0 Then
                grid(row*NCOLS+col) = tableau(row, col)
                grid(row*NCOLS+col+1) = tableau(row, col+1)
                used(usedCount) = d
                moves(moveCount).r1 = row : moves(moveCount).c1 = col
                moves(moveCount).r2 = row : moves(moveCount).c2 = col+1
                solveRec tableau(), grid(), used(), usedCount+1, moves(), moveCount+1, solutions(), num_solutions
                grid(row*NCOLS+col) = -1
                grid(row*NCOLS+col+1) = -1
            End If
        End If
    End If
    
    ' Vertical
    If row < NROWS-1 Then
        If grid((row+1)*NCOLS+col) = -1 Then
            Dim d As Domino = makeDomino(tableau(row, col), tableau(row+1, col))
            If dominoUsed(d, used(), usedCount) = 0 Then
                grid(row*NCOLS+col) = tableau(row, col)
                grid((row+1)*NCOLS+col) = tableau(row+1, col)
                used(usedCount) = d
                moves(moveCount).r1 = row : moves(moveCount).c1 = col
                moves(moveCount).r2 = row+1 : moves(moveCount).c2 = col
                solveRec tableau(), grid(), used(), usedCount+1, moves(), moveCount+1, solutions(), num_solutions
                grid(row*NCOLS+col) = -1
                grid((row+1)*NCOLS+col) = -1
            End If
        End If
    End If
End Sub

Sub solve(tableau() As Integer, solutions() As Solution, Byref num_solutions As Integer)
    Dim As Integer grid(NROWS*NCOLS-1)
    For i As Integer = 0 To NROWS*NCOLS-1
        grid(i) = -1
    Next
    Dim As Domino used(MAXDOMINOS-1)
    Dim As Move moves(MAXDOMINOS-1)
    num_solutions = 0
    solveRec tableau(), grid(), used(), 0, moves(), 0, solutions(), num_solutions
End Sub

' Main progranm
Dim Shared As Solution solutions1(0 To MAXSOLUTIONS-1)
Dim As Integer num_solutions1
Dim Shared As Solution solutions2(0 To MAXSOLUTIONS-1)
Dim As Integer num_solutions2

Dim As Double t1 = Timer
solve(tableau1(), solutions1(), num_solutions1)
Dim As Double t2 = Timer
solve(tableau2(), solutions2(), num_solutions2)
Dim As Double t3 = Timer

If num_solutions1 > 0 Then printDominoLayout(solutions1(0).moves(), tableau1())
Print num_solutions1; " layouts found."; Iif(num_solutions1 > 1, " (first one shown)", "")
Print Using "Took #.##### seconds."; t2-t1
Print

If num_solutions2 > 0 Then printDominoLayout(solutions2(0).moves(), tableau2())
Print num_solutions2; " layouts found."; Iif(num_solutions2 > 1, " (first one shown)", "")
Print Using "Took #.##### seconds."; t3-t2
Print

Sleep
Output:
0+5     1+3     2   2+3     1
                +           +
0   5+5     0   5   2+4     6
+           +
4   3   0   3   6+6     2   0
    +   +               +   +
0   6   2   3+5     1   2   6
+                   +
1   1   3   0+0     2   4   5
    +   +               +   +
2   1   4   3+3     4   6   6
+                   +
6   4+5     1+5     4   1+4

 1 layouts found.
Took 0.00038 seconds.

6+4     2+2     0+6     5   0
                        +   +
1+6     2+3     4+1     4   3

2+1     0+2     3+5     5+1

1+3     5+0     5+6     1+0

4+2     6   0   4+0     1+1
        +   +
4+4     2   0   5   3+6     3
                +           +
6+6     5+2     5   3+3     4

 2025 layouts found. (first one shown)
Took 0.00581 seconds.

Extra credit task

Library: big_integer
#include once "big_int/big_integer.bi"

Function DominoTilingCount(m As Integer, n As Integer) As Bigint
    Dim As Double prod = 1.0
    Dim As Integer jmax = m \ 2, kmax = n \ 2
    
    For j As Integer = 1 To jmax
        For k As Integer = 1 To kmax
            Dim cj As Double = Cos(j * 4 * Atn(1) / (m + 1))
            Dim ck As Double = Cos(k * 4 * Atn(1) / (n + 1))
            Dim term As Double = 4 * cj * cj + 4 * ck * ck
            prod *= term
        Next
    Next
    
    ' Round only at the end and convert to Bigint
    Return Bigint(Clng(prod + 0.5))
End Function

Function FactBig(n As Integer) As Bigint
    Dim As Bigint res = 1
    For i As Integer = 2 To n
        res *= i
    Next
    Return res
End Function

' Main program
Dim As Double t4 = Timer

Dim As BigInt arrang = DominoTilingCount(7, 8)
Dim As BigInt perms = FactBig(28)
Dim As BigInt flips = 2^28

Print "Arrangements ignoring values: "; arrang
Print "Permutations of 28 dominos: "; perms
Print "Permuted arrangements ignoring flipping dominos: "; arrang * perms
Print "Possible flip configurations: "; flips
Print "Possible permuted arrangements with flips: "; flips * arrang * perms

Dim As Double t5 = Timer
Print Using !"\nTook #.##### seconds."; t5-t4

Sleep
Output:
Arrangements ignoring values: +1292698
Permutations of 28 dominos: +304888344611713860501504000000
Permuted arrangements ignoring flipping dominos: +394128553302873284042573217792000000
Possible flip configurations: +268435456
Possible permuted arrangements with flips: +105798077928477096112185665131382833152000000

Took 0.00247 seconds.
import java.awt.Point;
import java.util.ArrayList;
import java.util.Collections;
import java.util.List;
import java.util.stream.Collectors;
import java.util.stream.IntStream;

public final class Dominos {

	public static void main(String[] args) {
		List.of( tableauOne, tableauTwo ).forEach( tableau -> {
			List<Pattern> patterns = findPatterns(tableau);
			System.out.println("Layouts found: " + patterns.size() + System.lineSeparator());
			printLayout(patterns.getFirst());
		} );
	}
	
	private static List<Pattern> findPatterns(int[][] tableau) {		
	    final int nRows = tableau.length;
	    final int nCols = tableau[0].length;
	    final int dominoCount = ( nRows * nCols ) / 2;
	    int[][] emptyTableau = IntStream.range(0, nRows)
	    	.mapToObj( i -> IntStream.range(0, nCols).map( j -> EMPTY ).toArray() ).toArray(int[][]::new);
	    List<Pattern> patterns = new ArrayList<Pattern>();
	    patterns.addLast( new Pattern(emptyTableau, new ArrayList<Domino>(), new ArrayList<Point>()) );
	    
	    while ( true ) {
	        List<Pattern> nextPatterns = new ArrayList<Pattern>();
	        for ( Pattern pattern : patterns ) {
	        	int[][] nextTableau = pattern.tableau;
	        	List<Domino> dominoes = pattern.dominoes;
	        	List<Point> points = pattern.points;	        	
	        	final int index = IntStream.range(0, nRows * nCols)
	        		.filter( i -> nextTableau[i / nCols][i % nCols] == EMPTY ).findFirst().orElse(EMPTY);	        
	        	if ( index == EMPTY ) {
	                continue;
	            }
	            
	            int row = index / nCols;
	            int col = index % nCols;            
	            if ( row + 1 < nRows && nextTableau[row + 1][col] == EMPTY ) {	               
        	        Domino domino = new Domino(tableau[row][col], tableau[row + 1][col]); 
        	        if ( ! dominoes.contains(domino) ) {        	        	
        	        	int[][] finalTableau = new int[nextTableau.length][];
        	        	for ( int i = 0; i < nextTableau.length; i++ ) {
        	        	      finalTableau[i] = new int[nextTableau[i].length];
        	        	      System.arraycopy(nextTableau[i], 0, finalTableau[i], 0, nextTableau[i].length);        	        		
        	        	}        	        	
        	        	finalTableau[row][col] = tableau[row][col];
        	        	finalTableau[row + 1][col] = tableau[row + 1][col];         	        	
        	        	List<Domino> nextDominoes = new ArrayList<Domino>(dominoes);
        	        	nextDominoes.addLast(domino);
        	        	List<Point> nextPoints = new ArrayList<Point>(points);
        	        	nextPoints.addAll(List.of( new Point(row, col), new Point(row + 1, col) ));
        	        	nextPatterns.addLast( new Pattern(finalTableau, nextDominoes, nextPoints) );                
        	        }
	            }	            
	            
	            if ( col + 1 < nCols && nextTableau[row][col + 1] == EMPTY ) {
	                Domino domino = new Domino(tableau[row][col], tableau[row][col + 1]);
	                if ( ! dominoes.contains(domino) ) {
	                	nextTableau[row][col] = tableau[row][col];
	                	nextTableau[row][col + 1] = tableau[row][col + 1];	                	                	
            	        dominoes.add(domino);
            	        points.addAll(List.of( new Point(row, col), new Point(row, col + 1) )); 
            	        nextPatterns.addLast( new Pattern(nextTableau, dominoes, points) );
	                }
	            }
	        }
	        
	        if ( nextPatterns.isEmpty() ) {  	
	        	break;
	        }
	        patterns = nextPatterns;
	        if ( patterns.getFirst().dominoes.size() == dominoCount ) {
	        	break;
	        }	                  
	    }
	    
	    return patterns;
	}	

	private static void printLayout(Pattern pattern) {	
		List<List<String>> output = IntStream.range(0, 2 * pattern.tableau.length)
		.mapToObj( i -> new ArrayList<String>(Collections.nCopies(2 * pattern.tableau[0].length - 1, " ")) )
		.collect(Collectors.toList());
		
	    for ( int i = 0; i < pattern.points.size() - 1; i += 2 ) {
	    	final int x1 = pattern.points.get(i).x;
	    	final int y1 = pattern.points.get(i).y;
	    	final int x2 = pattern.points.get(i + 1).x;
	    	final int y2 = pattern.points.get(i + 1).y;
	    	final int n1 = pattern.tableau[x1][y1];
	    	final int n2 = pattern.tableau[x2][y2];
	        output.get(2 * x1).set(2 * y1, String.valueOf(n1));
	        output.get(2 * x2).set(2 * y2, String.valueOf(n2));
	        if ( x1 == x2 ) {
	        	output.get(2 * x1).set(2 * y1 + 1, "+");
	        } else if ( y1 == y2 ) {
	        	output.get(2 * x1 + 1).set(2 * y1, "+");
	        }
	    }
	    
	    output.forEach( i -> System.out.println(i.stream().collect(Collectors.joining())) );
	} 
	
	private static class Domino {
	
		public Domino(int aOne, int aTwo) {
			one = Math.min(aOne, aTwo);
			two = Math.max(aOne, aTwo);
		}
		
		@Override
		public int hashCode() {
			return 2 * one + 3 * two;
		}
		
		@Override
		public boolean equals(Object other) {
			return switch ( other ) {
				case Domino domino -> one == domino.one && two == domino.two;
				case Object object -> false;
			};
		}
		
		private int one, two;		
	
	}
	
	private static final int EMPTY = -1;
	
	private static record Pattern(int[][] tableau, List<Domino> dominoes, List<Point> points) {}
	
	private static final int[][] tableauOne = { { 0, 5, 1, 3, 2, 2, 3, 1 },
												{ 0, 5, 5, 0, 5, 2, 4, 6 },
												{ 4, 3, 0, 3, 6, 6, 2, 0 },
	                                            { 0, 6, 2, 3, 5, 1, 2, 6 },
	                                            { 1, 1, 3, 0, 0, 2, 4, 5 },
	                                            { 2, 1, 4, 3, 3, 4, 6, 6 },
	                                            { 6, 4, 5, 1, 5, 4, 1, 4 } }; 
	
	private static final int[][] tableauTwo = { { 6, 4, 2, 2, 0, 6, 5, 0 },
	                                            { 1, 6, 2, 3, 4, 1, 4, 3 },
	                                            { 2, 1, 0, 2, 3, 5, 5, 1 },
	                                            { 1, 3, 5, 0, 5, 6, 1, 0 },
	                                            { 4, 2, 6, 0, 4, 0, 1, 1 },
	                                            { 4, 4, 2, 0, 5, 3, 6, 3 },
	                                            { 6, 6, 5, 2, 5, 3, 3, 4 } };	

}
Output:
Layouts found: 1

0+5 1+3 2 2+3 1
        +     +
0 5+5 0 5 2+4 6
+     +        
4 3 0 3 6+6 2 0
  + +       + +
0 6 2 3+5 1 2 6
+         +    
1 1 3 0+0 2 4 5
  + +       + +
2 1 4 3+3 4 6 6
+         +    
6 4+5 1+5 4 1+4
               
Layouts found: 2025

6 4 2 2 0 6+5 0
+ + + + +     +
1 6 2 3 4 1+4 3
               
2 1 0 2 3+5 5 1
+ + + +     + +
1 3 5 0 5 6 1 0
        + +    
4 2 6 0 4 0 1+1
+ + + +        
4 4 2 0 5 3 6 3
        + + + +
6+6 5+2 5 3 3 4

Extra Credit Task

import java.math.BigInteger;

public final class DominosExtraCredit {

	public static void main(String[] args) {
		BigInteger arrangements = BigInteger.valueOf(dominoTilingCount(7, 8));
		BigInteger permutations = factorial(BigInteger.valueOf(28));
		BigInteger flips = BigInteger.TWO.pow(28);
		
		System.out.println("Arrangements ignoring values: " + arrangements);
		System.out.println("Permutations of 28 dominos: " + permutations);
		System.out.println("Permuted arrangements ignoring flipping dominos: "
						   + permutations.multiply(arrangements));
		System.out.println("Possible flip configurations: " + flips);
		System.out.println("Possible permuted arrangements with flips: "
						   + permutations.multiply(flips).multiply(arrangements));
	}
	
	private static int dominoTilingCount(int rows, int cols) {
	    double product = 1.0;
	    for ( int i = 1; i <= ( rows + 1 ) / 2; i++ ) {
	    	for ( int j = 1; j <= ( cols + 1 ) / 2; j++ ) {
	    		final double cosRows = Math.cos(Math.PI * i / ( rows + 1 ));
	    		final double cosCols = Math.cos(Math.PI * j / ( cols + 1 ));
	    		product *= ( cosRows * cosRows + cosCols * cosCols ) * 4;
	    	}
	    }
	    return (int) product;
	}
	
	private static BigInteger factorial(BigInteger number) {
		if ( number.equals(BigInteger.TWO) ) {
			return number;
		}
		return number.multiply(factorial(number.subtract(BigInteger.ONE)));
	}

}
Output:
Arrangements ignoring values: 1292697
Permutations of 28 dominos: 304888344611713860501504000000
Permuted arrangements ignoring flipping dominos: 394128248414528672328712716288000000
Possible flip configurations: 268435456
Possible permuted arrangements with flips: 105797996085635281181632579889767907328000000


Works with: NodeJS 16.14.2
Translation of: Java
const EMPTY = -1;

const tableauOne = [
    [0, 5, 1, 3, 2, 2, 3, 1],
    [0, 5, 5, 0, 5, 2, 4, 6],
    [4, 3, 0, 3, 6, 6, 2, 0],
    [0, 6, 2, 3, 5, 1, 2, 6],
    [1, 1, 3, 0, 0, 2, 4, 5],
    [2, 1, 4, 3, 3, 4, 6, 6],
    [6, 4, 5, 1, 5, 4, 1, 4]
];

const tableauTwo = [
    [6, 4, 2, 2, 0, 6, 5, 0],
    [1, 6, 2, 3, 4, 1, 4, 3],
    [2, 1, 0, 2, 3, 5, 5, 1],
    [1, 3, 5, 0, 5, 6, 1, 0],
    [4, 2, 6, 0, 4, 0, 1, 1],
    [4, 4, 2, 0, 5, 3, 6, 3],
    [6, 6, 5, 2, 5, 3, 3, 4]
];

class Domino {
    constructor(aOne, aTwo) {
        this.one = Math.min(aOne, aTwo);
        this.two = Math.max(aOne, aTwo);
    }

    equals(other) {
        return this.one === other.one && this.two === other.two;
    }
}

class Point {
    constructor(aX, aY) {
        this.x = aX;
        this.y = aY;
    }

    equals(other) {
        return this.x === other.x && this.y === other.y;
    }
}

class Pattern {
    constructor(tableau, dominoes, points) {
        this.tableau = tableau;
        this.dominoes = dominoes;
        this.points = points;
    }
}

function firstEmptyCell(tableau) {
    for (let row = 0; row < tableau.length; ++row) {
        for (let col = 0; col < tableau[0].length; ++col) {
            if (tableau[row][col] === EMPTY) {
                return row * tableau[0].length + col;
            }
        }
    }
    return EMPTY;
}

function printLayout(pattern) {
    const output = Array(2 * pattern.tableau.length)
        .fill()
        .map(() => Array(2 * pattern.tableau[0].length - 1).fill(' '));

    for (let i = 0; i < pattern.points.length - 1; i += 2) {
        const x1 = pattern.points[i].x;
        const y1 = pattern.points[i].y;
        const x2 = pattern.points[i + 1].x;
        const y2 = pattern.points[i + 1].y;
        const n1 = pattern.tableau[x1][y1];
        const n2 = pattern.tableau[x2][y2];

        output[2 * x1][2 * y1] = String.fromCharCode('0'.charCodeAt(0) + n1);
        output[2 * x2][2 * y2] = String.fromCharCode('0'.charCodeAt(0) + n2);

        if (x1 === x2) {
            output[2 * x1][2 * y1 + 1] = '+';
        } else if (y1 === y2) {
            output[2 * x1 + 1][2 * y1] = '+';
        }
    }

    for (const line of output) {
        console.log(line.join(''));
    }
}

function findPatterns(tableau) {
    const nRows = tableau.length;
    const nCols = tableau[0].length;
    const dominoCount = (nRows * nCols) / 2;

    const emptyTableau = Array(nRows)
        .fill()
        .map(() => Array(nCols).fill(EMPTY));

    let patterns = [new Pattern(emptyTableau, [], [])];

    while (true) {
        const nextPatterns = [];

        for (const pattern of patterns) {
            const nextTableau = pattern.tableau.map(row => [...row]);
            const dominoes = [...pattern.dominoes];
            const points = [...pattern.points];

            const index = firstEmptyCell(nextTableau);
            if (index === EMPTY) continue;

            const row = Math.floor(index / nCols);
            const col = index % nCols;

            // Check down
            if (row + 1 < nRows && nextTableau[row + 1][col] === EMPTY) {
                const domino = new Domino(tableau[row][col], tableau[row + 1][col]);
                if (!dominoes.some(d => d.equals(domino))) {
                    const finalTableau = nextTableau.map(row => [...row]);
                    finalTableau[row][col] = tableau[row][col];
                    finalTableau[row + 1][col] = tableau[row + 1][col];

                    const nextDominoes = [...dominoes, domino];
                    const nextPoints = [...points, new Point(row, col), new Point(row + 1, col)];

                    nextPatterns.push(new Pattern(finalTableau, nextDominoes, nextPoints));
                }
            }

            // Check right
            if (col + 1 < nCols && nextTableau[row][col + 1] === EMPTY) {
                const domino = new Domino(tableau[row][col], tableau[row][col + 1]);
                if (!dominoes.some(d => d.equals(domino))) {
                    nextTableau[row][col] = tableau[row][col];
                    nextTableau[row][col + 1] = tableau[row][col + 1];

                    const nextDominoes = [...dominoes, domino];
                    const nextPoints = [...points, new Point(row, col), new Point(row, col + 1)];

                    nextPatterns.push(new Pattern(nextTableau, nextDominoes, nextPoints));
                }
            }
        }

        if (nextPatterns.length === 0) break;
        patterns = nextPatterns;
        if (patterns[0].dominoes.length === dominoCount) break;
    }

    return patterns;
}

// Main execution
for (const tableau of [tableauOne, tableauTwo]) {
    const patterns = findPatterns(tableau);
    console.log(`Layouts found: ${patterns.length}`);
    if (patterns.length > 0) {
        printLayout(patterns[0]);
    }
}
Output:
Layouts found: 1
0+5 1+3 2 2+3 1
        +     +
0 5+5 0 5 2+4 6
+     +        
4 3 0 3 6+6 2 0
  + +       + +
0 6 2 3+5 1 2 6
+         +    
1 1 3 0+0 2 4 5
  + +       + +
2 1 4 3+3 4 6 6
+         +    
6 4+5 1+5 4 1+4
               
Layouts found: 2025
6 4 2 2 0 6+5 0
+ + + + +     +
1 6 2 3 4 1+4 3
               
2 1 0 2 3+5 5 1
+ + + +     + +
1 3 5 0 5 6 1 0
        + +    
4 2 6 0 4 0 1+1
+ + + +        
4 4 2 0 5 3 6 3
        + + + +
6+6 5+2 5 3 3 4
               



const tableau = [
    0 5 1 3 2 2 3 1;
    0 5 5 0 5 2 4 6;
    4 3 0 3 6 6 2 0;
    0 6 2 3 5 1 2 6;
    1 1 3 0 0 2 4 5;
    2 1 4 3 3 4 6 6;
    6 4 5 1 5 4 1 4
]

const dominoes = [(i, j) for i in 0:size(tableau)[1]-1, j in 0:size(tableau)[2]-1 if i <= j]
sorted(dom) = first(dom) > last(dom) ? reverse(dom) : dom

"""
`patterns` contains solution(s), each containing a partially completed grid, the
dominos used, and steps taken to get to that point in the grid. Proceed via iterating
through possible tile placements from upper left to lower right, adding horizontal and
vertical tile placements, dropping those that require more than one of the same domino.
Consolidate in `patterns`` the newly lengthened layouts each step as moves are added.
"""
function findlayouts(tab = tableau, doms = dominoes)
    nrows, ncols = size(tab)
    patterns = [(zero(tab) .- 1, Tuple{Int, Int}[], Int[])]
    while true
        newpat = empty(patterns)
        for (ut, ud, up) in patterns
            pos = findfirst(x -> x == -1, ut)
            pos == nothing && continue
            row, col = Tuple(pos)
            if row < nrows && ut[row + 1, col] == -1 &&
               !(sorted((tab[row, col], tab[row + 1, col])) in ud)
                newut = copy(ut)
                newut[row:row+1, col] .= tab[row:row+1, col]
                push!(newpat, (newut, [ud; sorted((tab[row, col], tab[row + 1, col]))],
                   [up; [row, col, row + 1, col]]))
            end
            if col < ncols && ut[row, col + 1] == -1 &&
               !(sorted((tab[row, col], tab[row, col + 1])) in ud)
                newut = copy(ut)
                newut[row, col:col+1] .= tab[row, col:col+1]
                push!(newpat, (newut, [ud; sorted((tab[row, col], tab[row, col + 1]))],
                   [up; [row, col, row, col + 1]]))
            end
        end
        isempty(newpat) && break
        patterns = newpat
        length(last(first(patterns))) == length(doms) && break
    end
    return patterns
end

function printlayout(pattern)
    tab, dom, pos = pattern
    bytes = [[UInt8(' ') for _ in 1:(size(tab)[2] * 2 - 1)] for _ in 1:size(tab)[1]*2]
    for idx in 1:4:length(pos)-1
        x1, y1, x2, y2 = pos[idx:idx+3]
        n1, n2 = tab[x1, y1], tab[x2, y2]
        bytes[x1 * 2 - 1][y1 * 2 - 1] = Char(n1 + '0')
        bytes[x2 * 2 - 1][y2 * 2 - 1] = Char(n2 + '0')
        if x1 == x2 # horizontal
            bytes[x1 * 2 - 1][y1 * 2] = Char('+')
        elseif y1 == y2 # vertical
            bytes[x1 * 2][y1 * 2 - 1] = Char('+')
        end
    end
    println(join(String.(bytes), "\n"))
end


for pat in findlayouts()
    printlayout(pat)
end
@time findlayouts()

const t2 = [
    6 4 2 2 0 6 5 0;
    1 6 2 3 4 1 4 3;
    2 1 0 2 3 5 5 1;
    1 3 5 0 5 6 1 0;
    4 2 6 0 4 0 1 1;
    4 4 2 0 5 3 6 3;
    6 6 5 2 5 3 3 4
]
@time lays = findlayouts(t2, dominoes)
printlayout(first(lays))
println(length(lays), " layouts found.")
Output:
0+5 1+3 2 2+3 1
        +     +
0 5+5 0 5 2+4 6
+     +
4 3 0 3 6+6 2 0
  + +       + +
0 6 2 3+5 1 2 6
+         +
1 1 3 0+0 2 4 5
  + +       + +
2 1 4 3+3 4 6 6
+         +
6 4+5 1+5 4 1+4

  0.000507 seconds (6.06 k allocations: 1.715 MiB)

  0.023503 seconds (92.66 k allocations: 35.817 MiB)
6 4 2 2 0 6+5 0
+ + + + +     +
1 6 2 3 4 1+4 3

2 1 0 2 3+5 5 1
+ + + +     + +
1 3 5 0 5 6 1 0
        + +
4 2 6 0 4 0 1+1
+ + + +
4 4 2 0 5 3 6 3
        + + + +
6+6 5+2 5 3 3 4

2025 layouts found.

Extra credit task

""" From https://en.wikipedia.org/wiki/Domino_tiling#Counting_tilings_of_regions
The number of ways to cover an m X n rectangle with m * n / 2 dominoes, calculated
independently by Temperley & Fisher (1961) and Kasteleyn (1961), is given by
"""
function dominotilingcount(m, n)
    return BigInt(
        floor(
            prod([
                prod([
                    big"4.0" * (cospi(j / (m + 1)))^2 + 4 * (cospi(k / (n + 1)))^2 for
                    k in 1:(n+1)÷2
                ]) for j in 1:(m+1)÷2
            ]),
        ),
    )
end

arrang = dominotilingcount(7, 8)
perms = factorial(big"28")
flips = 2^28

println("Arrangements ignoring values: $arrang")
println("Permutations of 28 dominos: ", perms)
println("Permuted arrangements ignoring flipping dominos: ", arrang * perms)
println("Possible flip configurations: $flips")
println("Possible permuted arrangements with flips: ", flips * arrang * perms)
Output:
Arrangements ignoring values: 1292697
Permutations of 28 dominos: 304888344611713860501504000000
Permuted arrangements ignoring flipping dominos: 394128248414528672328712716288000000   
Possible flip configurations: 268435456
Possible permuted arrangements with flips: 105797996085635281181632579889767907328000000
ClearAll[VisualizeState]
VisualizeState[sol_List, tab_List] := Module[{rects},
  rects = Apply[Rectangle[#1 - {1, 1} 0.5, #2 + {1, 1} 0.5] &, sol[[All, 2]], {1}];
  Graphics[{FaceForm[], EdgeForm[Black], rects, MapIndexed[Text[Style[#1, 14, Black], #2] &, tab, {2}]}]
]
ClearAll[FindSolutions]
FindSolutions[tab_] := Module[{poss, possshort, possshorti, posssets, posss, sols},
  poss = Catenate[MapIndexed[#2 &, tab, {2}]];
  possshort = Thread[poss -> Range[Length[poss]]];
  possshorti = Thread[Range[Length[poss]] -> poss];
  posssets = Select[Subsets[poss, {2}], Apply[ManhattanDistance]/*EqualTo[1]];
  posssets = {# /. possshort, Sort[Extract[tab, #]]} & /@ posssets;
  posss = GatherBy[posssets, Last];
  posss = #[[1, 2]] -> #[[All, 1]] & /@ posss;
  posss //= SortBy[Last/*Length];
  
  sols = {};
  ClearAll[RecursePlaceDomino];
  RecursePlaceDomino[placed_List, left_List] := Module[{newplaced, sortedleft, newleft, next},
    If[Length[left] == 0,
     AppendTo[sols, placed];
     ,
     sortedleft = SortBy[left, Last/*Length];
     next = sortedleft[[1]];
     Do[
      newplaced = Append[placed, next[[1]] -> n];
      newleft = Drop[sortedleft, 1];
      newleft[[All, 2]] = newleft[[All, 2]] /. {___, Alternatives @@ n, ___} :> Sequence[];
      If[AnyTrue[newleft[[All, 2]], Length/*EqualTo[0]], Continue[]];
      RecursePlaceDomino[newplaced, newleft]
      ,
      {n, next[[2]]}
      ];
     ]
    ];
  RecursePlaceDomino[{}, posss];
  sols[[All, All, 2]] = sols[[All, All, 2]] /. possshorti;
  sols
]

tab = {{6, 2, 1, 0, 4, 0, 0}, {4, 1, 1, 6, 3, 5, 5}, {5, 4, 3, 2, 0, 
    5, 1}, {1, 3, 0, 3, 3, 0, 3}, {5, 3, 0, 5, 6, 5, 2}, {4, 4, 2, 1, 
    6, 2, 2}, {1, 6, 4, 2, 2, 4, 3}, {4, 6, 5, 6, 0, 6, 1}};
sols = FindSolutions[tab];
Length[sols]
VisualizeState[sols[[1]], tab]
tab = {{6, 4, 4, 1, 2, 1, 6}, {6, 4, 2, 3, 1, 6, 4}, {5, 2, 6, 5, 0, 
    2, 2}, {2, 0, 0, 0, 2, 3, 2}, {5, 5, 4, 5, 3, 4, 0}, {3, 3, 0, 6, 
    5, 1, 6}, {3, 6, 1, 1, 5, 4, 5}, {4, 3, 1, 0, 1, 3, 0}};
sols = FindSolutions[tab];
Length[sols]
VisualizeState[sols[[1]], tab]
Output:
1
[Graphical output of the solution]
2025
[Graphical output of the first solution]

Extra credit task

ClearAll[DominoTilingCount]
DominoTilingCount[m_, n_] := Module[{},
  Round[Product[
    4 Cos[(Pi j)/(m + 1)]^2 + 4 Cos[(Pi k)/(n + 1)]^2,
     {j, Ceiling[m/2]},
     {k, Ceiling[n/2]}
   ]]
]
arrangements = DominoTilingCount[7, 8] // Round;
permutations = 28!;
flips = 2^28;
Print["Arrangements ignoring values: ", arrangements]
Print["Permutations of 28 dominos: ", permutations]
Print["Permuted arrangements ignoring flipping dominos: ", arrangements*permutations]
Print["Possible flip configurations: ", flips]
Print["Possible permuted arrangements with flips: ", flips*arrangements*permutations]
Output:
Arrangements ignoring values: 1292697
Permutations of 28 dominos: 304888344611713860501504000000
Permuted arrangements ignoring flipping dominos: 394128248414528672328712716288000000
Possible flip configurations: 268435456
Possible permuted arrangements with flips: 105797996085635281181632579889767907328000000
Translation of: Wren
import std/[monotimes, sequtils, strformat, strutils, times]

const
  None = -1
  NoPosition = (None, None)

type
  Value = None..6
  Domino = (Value, Value)
  Position = (int, int)
  Tableau = array[7, array[8, Value]]
  Pattern = (seq[seq[Value]], seq[Domino], seq[Position])

const

  Tableau1 = [[Value 0, 5, 1, 3, 2, 2, 3, 1],
              [Value 0, 5, 5, 0, 5, 2, 4, 6],
              [Value 4, 3, 0, 3, 6, 6, 2, 0],
              [Value 0, 6, 2, 3, 5, 1, 2, 6],
              [Value 1, 1, 3, 0, 0, 2, 4, 5],
              [Value 2, 1, 4, 3, 3, 4, 6, 6],
              [Value 6, 4, 5, 1, 5, 4, 1, 4]]

  Tableau2 = [[Value 6, 4, 2, 2, 0, 6, 5, 0],
              [Value 1, 6, 2, 3, 4, 1, 4, 3],
              [Value 2, 1, 0, 2, 3, 5, 5, 1],
              [Value 1, 3, 5, 0, 5, 6, 1, 0],
              [Value 4, 2, 6, 0, 4, 0, 1, 1],
              [Value 4, 4, 2, 0, 5, 3, 6, 3],
              [Value 6, 6, 5, 2, 5, 3, 3, 4]]


func domino(v1, v2: Value): Domino =
  if v1 > v2: (v2, v1) else: (v1, v2)

func findLayouts(tab: Tableau; domCount: Positive): seq[Pattern] =
  let nrows = tab.len
  let ncols = tab[0].len
  var m = newSeqWith(nrows, repeat(Value None, ncols))
  result = @[(m, newSeq[Domino](), newSeq[Position]())]
  while true:
    var newpats: seq[Pattern]
    for pat in result:
      let (ut, ud, up) = pat
      var pos = NoPosition
      block Outer:
        for j in 0..<ncols:
          for i in 0..<nrows:
            if ut[i][j] == None:
              pos = (i, j)
              break Outer
      if pos == NoPosition:
        continue
      let (row , col) = pos
      if row < nrows - 1 and ut[row+1][col] == None:
        let dom = domino(tab[row][col], tab[row+1][col])
        if dom notin ud:
          var newUt = ut
          newut[row][col] = tab[row][col]
          newut[row+1][col] = tab[row+1][col]
          newpats.add (newut,
                       ud & domino(tab[row][col], tab[row+1][col]),
                       up & @[(row, col), (row+1, col)])
      if col < ncols - 1 and ut[row][col+1] == None:
        let dom = domino(tab[row][col], tab[row][col+1])
        if dom notin ud:
          var newUt = ut
          newut[row][col] = tab[row][col]
          newut[row][col+1] = tab[row][col+1]
          newpats.add (newut,
                       ud & domino(tab[row][col], tab[row][col+1]),
                       up & @[(row, col), (row, col+1)])
    if newPats.len == 0: break
    result = move(newPats)
    if result[0][2].len == domCount:
      break

proc printLayout(pattern: Pattern) =
  let (tab, _, pos) = pattern
  var output = newSeqWith(2 * tab.len, repeat(' ', tab[0].len * 2 - 1))
  var idx = 0
  while idx < pos.len - 1:
    let
      (x1, y1) = pos[idx]
      (x2, y2) = pos[idx+1]
    let
      n1 = tab[x1][y1]
      n2 = tab[x2][y2]
    output[x1*2][y1*2] = chr(48 + n1)
    output[x2*2][y2*2] = chr(48 + n2)
    if x1 == x2:
      output[x1*2][y1*2+1] = '+'
    elif y1 == y2:
      output[x1*2+1][y1*2] = '+'
    inc idx, 2
  for i in 0..output.high:
    echo output[i]

var domCount: Positive
for j in 0..Tableau1[0].high:
  for i in 0..Tableau1.high:
    if i <= j:
      inc domCount

for t in [Tableau1, Tableau2]:
  let start = getMonoTime()
  let lays = t.findLayouts(domCount)
  lays[0].printLayout()
  let lc = lays.len
  let pl = if lc > 1: "s" else: ""
  let fo = if lc > 1: " (first one shown)" else: ""
  echo &"{lc} layout{pl} found{fo}."
  echo &"Took {(getMonoTime() - start).inMilliseconds} ms.\n"
Output:
0+5 1+3 2 2+3 1
        +     +
0 5+5 0 5 2+4 6
+     +        
4 3 0 3 6+6 2 0
  + +       + +
0 6 2 3+5 1 2 6
+         +    
1 1 3 0+0 2 4 5
  + +       + +
2 1 4 3+3 4 6 6
+         +    
6 4+5 1+5 4 1+4
               
1 layout found.
Took 1 ms.

6 4 2 2 0 6+5 0
+ + + + +     +
1 6 2 3 4 1+4 3
               
2 1 0 2 3+5 5 1
+ + + +     + +
1 3 5 0 5 6 1 0
        + +    
4 2 6 0 4 0 1+1
+ + + +        
4 4 2 0 5 3 6 3
        + + + +
6+6 5+2 5 3 3 4
               
2025 layouts found (first one shown).
Took 42 ms.

Extra credit task

Translation of: Wren
Library: Nim-Integers
import std/[math, monotimes, strutils, times]
import integers

func dominoTilingCount(m, n: Positive): int =
  var prod = 1.0
  for j in 1..((m + 1) div 2):
    for k in 1..((n + 1) div 2):
      let cm = cos(PI * (j / (m + 1)))
      let cn = cos(PI * (k / (n + 1)))
      prod *= (cm * cm + cn * cn) * 4
  result = int(prod)

let
  start = getMonoTime()
  arrang = dominoTilingCount(7, 8)
  perms = factorial(28)
  flips = 1 shl 28

echo "Arrangements ignoring values: ", insertSep($arrang)
echo "Permutations of 28 dominos: ", insertSep($perms)
echo "Permuted arrangements ignoring flipping dominos: ", insertSep($(perms * arrang))
echo "Possible flip configurations: ", insertSep($flips)
echo "Possible permuted arrangements with flips: ", insertSep($(perms * flips * arrang))
echo "\nTook $# µs.".format((getMonoTime() - start).inMicroseconds)
Output:
Arrangements ignoring values: 1_292_697
Permutations of 28 dominos: 304_888_344_611_713_860_501_504_000_000
Permuted arrangements ignoring flipping dominos: 394_128_248_414_528_672_328_712_716_288_000_000
Possible flip configurations: 268_435_456
Possible permuted arrangements with flips: 105_797_996_085_635_281_181_632_579_889_767_907_328_000_000

Took 107 µs.
Library: ntheory
use v5.36;

# NB: not 'blank' lines, need to be full-width white-space
my $grid1 = <<END;
0 5 1 3 2 2 3 1
               
0 5 5 0 5 2 4 6
               
4 3 0 3 6 6 2 0
               
0 6 2 3 5 1 2 6
               
1 1 3 0 0 2 4 5
               
2 1 4 3 3 4 6 6
               
6 4 5 1 5 4 1 4
END

$grid2 = <<END;
0 0 0 1 1 1 0 2
               
1 2 2 2 0 3 1 3
               
2 3 3 3 0 4 1 4
               
2 4 3 4 4 4 0 5
               
1 5 2 5 3 5 4 5
               
5 5 0 6 1 6 2 6
               
3 6 4 6 5 6 6 6
END

my $grid3 = <<END;
0 0 1 1
       
0 2 2 2
       
1 2 0 1
END

sub find ($rows, $cols, $x, $y, $try, $orig) {
    state $solution;
    my $sum = $rows + $cols;
    my $gap = qr/(.{$sum}) (.{$sum})/s;

    if( $x > $y ) {
        $x = 0;
        $solution = $orig |. $try and return if ++$y == $rows;  # solved!
    }

    while ( $try =~ /(?=(?|$x $gap $y|$y $gap $x))/gx ) {       # vertical
        my $new = $try;
        substr $new, $-[0], 2*($rows+$cols)+3, " $1+$2 ";
        find($rows, $cols, $x + 1, $y, $new, $orig );
    }

    while ( $try =~ /(?=$x $y|$y $x)/g ) {                      # horizontal
        my $new = $try;
        substr $new, $-[0], 3, ' + ';
        find($rows, $cols, $x + 1, $y, $new, $orig );
    }

    $solution
}

say find(7, 8, 0, 0, $grid1, $grid1 ) . "\n=======\n\n";
say find(7, 8, 0, 0, $grid2, $grid2 ) . "\n=======\n\n";
say find(3, 4, 0, 0, $grid3, $grid3 ) . "\n=======\n\n";

use constant PI => 2*atan2(1,0);
use ntheory 'factorial';

sub comma { reverse ((reverse shift) =~ s/(.{3})/$1,/gr) =~ s/^,//r }

my($m,$n, $arrangements)  = (7,8, 1);
for my $j (1 .. $m/2) {
  for my $k (1 .. $n/2) {
    $arrangements *= 4*cos(PI*$j/($m+1))**2 + 4*cos(PI*$k/($n+1))**2
  }
}

printf "%32s:%60s\n", 'Arrangements ignoring values',     comma    $arrangements;
printf "%32s:%60s\n", 'Permutations of 28 dominos',       comma my $permutations = factorial 28;
printf "%32s:%60s\n", 'Flip configurations',              comma my $flips = 2**28;
printf "%32s:%60s\n", 'Permuted arrangements with flips', comma    $flips * $permutations * $arrangements;
Output:
0+5 1+3 2 2+3 1
        +     +
0 5+5 0 5 2+4 6
+     +        
4 3 0 3 6+6 2 0
  + +       + +
0 6 2 3+5 1 2 6
+         +    
1 1 3 0+0 2 4 5
  + +       + +
2 1 4 3+3 4 6 6
+         +    
6 4+5 1+5 4 1+4

=======

0+0 0+1 1+1 0+2
               
1 2+2 2 0+3 1+3
+     +        
2 3+3 3 0 4 1+4
        + +    
2+4 3+4 4 4 0+5
               
1 5 2+5 3+5 4+5
+ +            
5 5 0+6 1+6 2 6
            + +
3+6 4+6 5+6 6 6

=======

0+0 1+1

0+2 2+2

1+2 0+1

=======

    Arrangements ignoring values:                                                   1,292,697
      Permutations of 28 dominos:                     304,888,344,611,713,860,501,504,000,000
             Flip configurations:                                                 268,435,456
Permuted arrangements with flips: 105,797,996,085,635,281,181,632,579,889,767,907,328,000,000
with javascript_semantics

function domino_set()
    sequence set = {}
    for i=0 to 6 do
        for j=i to 6 do
            set = append(set,{i,j})
        end for
    end for
    return set
end function

sequence set = domino_set(),
        used = repeat(repeat(false,7),7),
        tags = shuffle(tagset(length(set))),
        grid

function unpack(sequence s)
    s = split(s,' ')
    s = apply(true,join,{s,'?'})
    s = join(s,"\n? ? ? ? ? ? ? ?\n")
    s = split(s,'\n')
    return s
end function

function clear(integer r, c, sequence s)
    if grid[r][c]!='?' then ?9/0 end if
    grid[r][c]='+'
    sequence res = {{r,c}}
    for i=1 to length(s) do
        {r,c} = s[i]
        if r>=1 and r<=13
        and c>=1 and c<=15 then
            integer prev = grid[r][c]
            if prev='?' then
                grid[r][c] = ' '
                res = append(res,{r,c})
            elsif prev!=' ' then
                ?9/0
            end if
        end if
    end for
    return res
end function

procedure restore(sequence s)
    for i=1 to length(s) do
        integer {r,c} = s[i]
        grid[r][c] = '?'
    end for
end procedure

function rand_grid(integer rem=28)
    if rem=0 then
        for r=1 to 13 by 2 do
            for c=2 to 14 by 2 do
                grid[r][c] = '?'
            end for
        end for
        for r=2 to 12 by 2 do
            for c=1 to 15 by 2 do
                grid[r][c] = '?'
            end for
        end for
        return grid
    end if
    for r=1 to 13 by 2 do
        for c=1 to 15 by 2 do
            bool flat = (c<15 and grid[r,c+1]='?'),
                 vert = (r<13 and grid[r+1,c]='?')
            sequence res = {},
                     opt = {}
            if flat then
                opt = append(opt,{{r,c+2,r,c+1},{{r,c-1},{r,c+3},{r-1,c},{r-1,c+2},{r+1,c},{r+1,c+2}}})
            end if
            if vert then
                opt = append(opt,{{r+2,c,r+1,c},{{r-1,c},{r,c-1},{r+2,c-1},{r,c+1},{r+2,c+1},{r+3,c}}})
            end if
            opt = shuffle(opt)
            for i=1 to length(opt) do
                integer {r2,c2,r3,c3} = opt[i][1]
                sequence tile = shuffle(set[tags[rem]])
                grid[r][c] = tile[1]+'0'
                grid[r2][c2] = tile[2]+'0'
                sequence reset = clear(r3,c3,opt[i][2])
                res = rand_grid(rem-1)
                if length(res) then return res end if
                restore(reset)
            end for
            if flat or vert then
                return {}
            end if
        end for
    end for
    return {}
end function

string soln1 = ""
string solnn = ""

function solve(integer rem=28)
    if rem=0 then
        solnn = join(grid,'\n')&"\n\n\n"
        if soln1 = "" then soln1 = solnn end if
        return 1
    end if
    for r=1 to 13 by 2 do
        for c=1 to 15 by 2 do
            bool flat = (c<15 and grid[r,c+1]='?'),
                 vert = (r<13 and grid[r+1,c]='?')
            integer count = 0
            if flat then
                integer {R,C} = sort({grid[r][c]-'0'+1,grid[r][c+2]-'0'+1})
                if not used[R][C] then
                    used[R][C] = true
                    sequence reset = clear(r,c+1,{{r,c-1},{r,c+3},{r-1,c},{r-1,c+2},{r+1,c},{r+1,c+2}})
                    count += solve(rem-1)
                    restore(reset)
                    used[R][C] = false
                end if
            end if
            if vert then
                integer {R,C} = sort({grid[r][c]-'0'+1,grid[r+2][c]-'0'+1})
                if not used[R][C] then
                    used[R][C] = true
                    sequence reset = clear(r+1,c,{{r-1,c},{r,c-1},{r+2,c-1},{r,c+1},{r+2,c+1},{r+3,c}})
                    count += solve(rem-1)
                    restore(reset)
                    used[R][C] = false
                end if
            end if
            if flat or vert then
                return count -- (may still be 0)
            end if
        end for
    end for
    return 0
end function

procedure test(sequence g)
    grid = g
    g = {}
    soln1 = ""
    solnn = ""
    atom t0 = time()
    integer n = solve()
    puts(1,soln1)
    if n>1 then
        if n>2 then puts(1,"...\n\n\n") end if
        puts(1,solnn)
    end if
    printf(1,"%d solution%s found (%s)\n\n\n",{n,iff(n=1?"":"s"),elapsed(time()-t0)})
end procedure
test(unpack("05132231 05505246 43036620 06235126 11300245 21433466 64515414"))
test(rand_grid())
Output:
0+5 1+3 2 2+3 1
        +     +
0 5+5 0 5 2+4 6
+     +
4 3 0 3 6+6 2 0
  + +       + +
0 6 2 3+5 1 2 6
+         +
1 1 3 0+0 2 4 5
  + +       + +
2 1 4 3+3 4 6 6
+         +
6 4+5 1+5 4 1+4


1 solution found (0s)


6+4 2+2 0+6 5 0
            + +
1+6 2+3 4+1 4 3

2+1 0+2 3+5 5+1

1+3 5+0 5+6 1+0

4+2 6 0 4+0 1+1
    + +
4+4 2 0 5 3+6 3
        +     +
6+6 5+2 5 3+3 4


...


6 4 2 2 0 6+5 0
+ + + + +     +
1 6 2 3 4 1+4 3

2 1 0 2 3+5 5 1
+ + + +     + +
1 3 5 0 5 6 1 0
        + +
4 2 6 0 4 0 1+1
+ + + +
4 4 2 0 5 3 6 3
        + + + +
6+6 5+2 5 3 3 4


2025 solutions found (0.0s)

Note that 2025 is not the maximum number of solutions or anything like that, just a higher than average result.

Extra credit

Translation of: Wren
with javascript_semantics
include mpfr.e

function domino_tiling_count(integer m=7, n=8)
    atom prod = 1
    for j=1 to ceil(m/2) do
        for k=1 to ceil(n/2) do
            atom cm = cos(PI * (j / (m + 1))),
                 cn = cos(PI * (k / (n + 1)))
            prod *= (cm*cm + cn*cn) * 4
        end for
    end for
    return floor(prod)
end function
 
atom start = time()
integer arrang = domino_tiling_count(),
        flips  = power(2,28)
mpz perms = mpz_init()
mpz_fac_ui(perms,28)
 
printf(1,"Arrangements ignoring values: %,d\n", arrang)
printf(1,"Permutations of 28 dominos: %s\n", mpz_get_str(perms,10,true))
mpz_mul_si(perms,perms,arrang)
printf(1,"Permuted arrangements ignoring flipping dominos: %s\n", mpz_get_str(perms,10,true))
printf(1,"Possible flip configurations: %,d\n", flips)
mpz_mul_si(perms,perms,flips)
printf(1,"Possible permuted arrangements with flips: %s\n", mpz_get_str(perms,10,true))
printf(1,"Took %s\n",elapsed(time()-start))
Output:
Arrangements ignoring values: 1,292,697
Permutations of 28 dominos: 304,888,344,611,713,860,501,504,000,000
Permuted arrangements ignoring flipping dominos: 394,128,248,414,528,672,328,712,716,288,000,000
Possible flip configurations: 268,435,456
Possible permuted arrangements with flips: 105,797,996,085,635,281,181,632,579,889,767,907,328,000,000
Took 0.1s
from typing import Final, Generator, Iterator, Optional
import random
from dataclasses import dataclass

# Use a cryptographically secure random number generator
rnd = random.SystemRandom()

# Sentinel value to represent an empty/unfilled cell in the grid during solving
EMPTY: Final = -1

# Predefined 8x8 grid of integers (0–6), representing a domino puzzle board
tableau = [
    [0, 5, 1, 3, 2, 2, 3, 1],
    [0, 5, 5, 0, 5, 2, 4, 6],
    [4, 3, 0, 3, 6, 6, 2, 0],
    [0, 6, 2, 3, 5, 1, 2, 6],
    [1, 1, 3, 0, 0, 2, 4, 5],
    [2, 1, 4, 3, 3, 4, 6, 6],
    [6, 4, 5, 1, 5, 4, 1, 4]
]
# Randomly generated 8x8 grid using values 0–6 for testing alternative puzzles
customTableau = [[rnd.randint(0, 6) for _ in range(8)] for _ in range(8)]

# # Represents a domino tile with two ends (a and b). 
# It is treated as unordered (i.e., (1,2) == (2,1)) via custom hash and equality.
@dataclass(frozen=True)
class Domino:
    a: int
    b: int
    
    # Hash based on sorted tuple so (a,b) and (b,a) are considered the same domino
    def __hash__(self) -> int:
        return hash(tuple(sorted((self.a, self.b))))
    
    # String representation for debugging/printing
    def __repr__(self) -> str:
        return f"({self.a},{self.b})"

# Represents a coordinate (x, y) on the grid
@dataclass
class Point:
    x: int
    y: int
    
    def __repr__(self) -> str:
        return f"({self.x},{self.y})"

# Bundles a solved layout: the filled grid, list of placed dominoes, and their positions
@dataclass
class Pattern:
    tableau: list[list[int]]  # Final filled grid (copied from solver state)
    dominoes: list[Domino]    # List of dominoes used (in placement order)
    points: list[Point]       # Paired list of points indicating domino placements

def findPatterns(
    source: list[list[int]],
    *,
    maxSolutions: Optional[int] = None,
    asGenerator: bool = False
) -> Iterator[Pattern] | list[Pattern]:
    """
    Solves a domino tiling puzzle where each domino must be unique (unordered pair),
    and every cell must be covered exactly once.

    Args:
        source: 2D grid of integers (0-6) representing the puzzle.
        maxSolutions: Optional limit on number of solutions to find.
        asGenerator: If True, returns a generator yielding solutions one by one;
                     otherwise returns a list of all found solutions (up to maxSolutions).

    Returns:
        Either a list of Pattern objects or a generator yielding them.
    """
    nRows = len(source)
    assert nRows > 0
    nCols = len(source[0])
    # Ensure all rows have the same length
    for row in source:
        if len(row) != nCols:
            raise ValueError("All rows must have same length")
    
    # Total number of dominoes needed to cover the board
    dominoGoal = (nRows * nCols) // 2
    # If total cells is odd, tiling is impossible
    if (nRows * nCols) % 2 != 0:
        return []
    
    # Working grid: EMPTY (-1) means not yet covered
    grid = [[EMPTY for _ in range(nCols)] for _ in range(nRows)]
    usedDominoes: set[Domino] = set()    # For O(1) uniqueness checks
    usedDominoesList: list[Domino] = []  # To preserve order of placement
    points: list[Point] = []             # Stores paired coordinates of domino placements
    solutions: list[Pattern] = []        # Accumulates found solutions
    
    def findFirstEmpty(startIndex: int = 0) -> Optional[int]:
        """Find the first empty cell in row-major order starting from `startIndex`."""
        total = nRows * nCols
        for i in range(startIndex, total):
            r = i // nCols
            c = i % nCols
            if grid[r][c] == EMPTY:
                return i
        return  # No empty cell found
    
    def collectSolution() -> Pattern:
        """Create a deep copy of the current solution state."""
        tableauCopy = [row[:] for row in grid]
        return Pattern(tableauCopy, usedDominoesList.copy(), points.copy())
    
    def dfs(startIndex: int = 0):
        """
        Depth-first search to place dominoes recursively.
        Stops early if `maxSolutions` is reached.
        """
        if maxSolutions is not None and len(solutions) >= maxSolutions:
            return
        
        index = findFirstEmpty(startIndex)
        if index is None:
            # All cells filled — check if we used exactly the right number of dominoes
            if len(usedDominoesList) == dominoGoal:
                solutions.append(collectSolution())
            return
        
        r = index // nCols
        c = index % nCols
        
        # Try placing a vertical domino (downwards)
        if r + 1 < nRows and grid[r+1][c] == EMPTY:
            d = Domino(source[r][c], source[r+1][c])
            if d not in usedDominoes:
                # Place domino
                grid[r][c] = source[r][c]
                grid[r+1][c] = source[r+1][c]
                usedDominoes.add(d)
                usedDominoesList.append(d)
                points.extend([Point(r, c), Point(r+1, c)])
                
                dfs(index+1)
                
                # Backtrack
                points.pop()
                points.pop()
                usedDominoesList.pop()
                usedDominoes.remove(d)
                grid[r][c] = EMPTY
                grid[r+1][c] = EMPTY
                
                # Early exit if enough solutions found
                if maxSolutions is not None and len(solutions) >= maxSolutions:
                    return
        
        # Try placing a horizontal domino (rightwards)
        if c + 1 < nCols and grid[r][c+1] == EMPTY:
            d = Domino(source[r][c], source[r][c+1])
            if d not in usedDominoes:
                # Place domino
                grid[r][c] = source[r][c]
                grid[r][c+1] = source[r][c+1]
                usedDominoes.add(d)
                usedDominoesList.append(d)
                points.extend([Point(r, c), Point(r, c+1)])
                
                dfs(index+1)
                
                # Backtrack
                points.pop()
                points.pop()
                usedDominoesList.pop()
                usedDominoes.remove(d)
                grid[r][c] = EMPTY
                grid[r][c+1] = EMPTY
                
                # Early exit if enough solutions found
                if maxSolutions is not None and len(solutions) >= maxSolutions:
                    return
    
    # If requested, return a generator instead of collecting all solutions upfront
    if asGenerator:
        def gen() -> Generator:
            def dfsGen(startIndex: int = 0) -> Generator:
                """Generator version of DFS that yields solutions as they are found."""
                index = findFirstEmpty(startIndex)
                if index is None:
                    if len(usedDominoesList) == dominoGoal:
                        yield collectSolution()
                    return
                r = index // nCols
                c = index % nCols
                
                # Vertical placement
                if r + 1 < nRows and grid[r+1][c] == EMPTY:
                    d = Domino(source[r][c], source[r+1][c])
                    if d not in usedDominoes:
                        grid[r][c] = source[r][c]
                        grid[r+1][c] = source[r+1][c]
                        usedDominoes.add(d)
                        usedDominoesList.append(d)
                        points.extend([Point(r, c), Point(r+1, c)])
                        
                        yield from dfsGen(index+1)
                        
                        # Backtrack
                        points.pop()
                        points.pop()
                        usedDominoesList.pop()
                        usedDominoes.remove(d)
                        grid[r][c] = EMPTY
                        grid[r+1][c] = EMPTY
                
                # Horizontal placement
                if c + 1 < nCols and grid[r][c+1] == EMPTY:
                    d = Domino(source[r][c], source[r][c+1])
                    if d not in usedDominoes:
                        grid[r][c] = source[r][c]
                        grid[r][c+1] = source[r][c+1]
                        usedDominoes.add(d)
                        usedDominoesList.append(d)
                        points.extend([Point(r, c), Point(r, c+1)])
                        
                        yield from dfsGen(index+1)
                        
                        # Backtrack
                        points.pop()
                        points.pop()
                        usedDominoesList.pop()
                        usedDominoes.remove(d)
                        grid[r][c] = EMPTY
                        grid[r][c+1] = EMPTY
                        
            yield from dfsGen(0)
        
        return gen()
    
    # Run standard DFS and return collected solutions
    dfs(0)
    return solutions

def printLayout(pattern: Pattern):
    """
    Pretty-prints a solved domino layout with ASCII art:
    - Numbers represent tile values.
    - '-' connects horizontally adjacent domino halves.
    - '|' connects vertically adjacent domino halves.
    """
    nRows = len(pattern.tableau)
    nCols = len(pattern.tableau[0]) if nRows > 0 else 0
    
    # Create a character grid large enough to show connections
    output = [[" " for _ in range(nCols*3-1)] for _ in range(nRows*2-1)]
    
    # Place numbers in the output grid
    for i in range(nRows):
        for j in range(nCols):
            val = pattern.tableau[i][j]
            ch = "?" if val == EMPTY else str(val)
            output[i*2][j*3] = ch
    
    # Draw connections between paired points (domino halves)
    for k in range(0, len(pattern.points), 2):
        if k + 1 >= len(pattern.points):
            break
        
        p0 = pattern.points[k]
        p1 = pattern.points[k+1]
        
        # Horizontal domino: same row, adjacent columns
        if p0.x == p1.x and p0.y + 1 == p1.y:
            output[p0.x*2][p0.y*3+1] = "-"
            output[p0.x*2][p0.y*3+2] = "-"
        # Vertical domino: same column, adjacent rows
        elif p0.y == p1.y and p0.x+1 == p1.x:
            output[p0.x*2+1][p0.y*3] = "|"
    
    # Print the final layout line by line
    for line in output:
        print("".join(line))

# Entry point: solve both the predefined and random tableaus, print first solution if found
if __name__ == "__main__":
    for t in [tableau, customTableau]:
        sols = findPatterns(t, maxSolutions=1)
        print(f"Layouts found: {len(sols)}")
        if len(sols) > 0:
            printLayout(sols[0])
Output:
Layouts found: 1
0  5  1  3  2  2  3--1 
|  |  |  |  |  |       
0  5  5  0  5  2  4  6 
                  |  | 
4  3  0  3  6--6  2  0 
|  |  |  |             
0  6  2  3  5  1  2  6 
            |  |  |  | 
1--1  3  0  0  2  4  5 
      |  |             
2  1  4  3  3  4  6  6 
|  |        |  |  |  | 
6  4  5--1  5  4  1  4 
Layouts found: 0

Since the custom tableau is randomly-generated, the program will sometimes show the solution for it, or will sometimes not.

Layouts found: 1
6--0  2  0  5--6  2  6 
      |  |        |  | 
1--5  4  6  0--2  6  6 
                       
6  2  5  2--1  2--3  4 
|  |  |              | 
1  2  5  0  3  6--2  0 
         |  |          
1  3  1  3  0  6  4--2 
|  |  |        |       
2  1  1  2--0  5  3--2 
                       
6  6  3  1  2  0  0  1 
|  |  |  |  |  |  |  | 
4  3  5  3  5  0  1  6 
sub domino (Str $s) { $s.comb.sort.join }

sub solve ( UInt $rows, UInt $cols, @tab ) {
    die unless @tab.elems == $rows*$cols;
    die unless @tab.elems %% 2;

    my @orientations = 'A' xx @tab.elems; # [A]vailable,[R]ight,[D]own,[L]eft,[U]p

    my SetHash $hand .= new: map &domino, [X~] @tab.unique xx 2;

    return gather {
        my sub place_domino_at_first_available_after ( UInt $pos_last = 0 ) {
            my $p1 = $pos_last + @orientations.skip($pos_last).first(:k, 'A');

            for (False, $p1+1    , <L R>),
                (True , $p1+$cols, <D U>) -> ($is_down, $p2, @two_letters) {

                next if $p2 > @tab.end or (!$is_down && $p2 %% $cols); # Board boundaries.

                my @two_cells := @orientations[$p1, $p2]; # Bind to candidate location

                next if @two_cells !eqv <A A>; # Both cells must be available for placement.

                my $dom = domino( @tab[$p1, $p2].join );

                $hand{$dom}-- or next;          # Take domino from hand, iff present
                @two_cells = @two_letters;      # Assign placement

                # Either the solution is complete, or we must recurse.
                if !$hand { take @orientations.clone }
                else      { place_domino_at_first_available_after($p1) }

                @two_cells = <A A>;             # Unassign placement
                $hand{$dom}++;                  # Restore domino into hand
            }
        }
        place_domino_at_first_available_after();
    }
}
sub bonus ( UInt $rows, UInt $cols ) {
    die "Both are odd, so no tiling possible" if $rows !%% 2 and $cols !%% 2;

    my $half = $rows * $cols div 2;

    my sub T ($N) { map { ( pi * $_/($N+1) ).cos² * 4 }, 1..($N/2).ceiling }

    my $arrangements = [*] T($rows) X+ T($cols);
    my $flips        =     2 ** $half;
    my $permutations = [*] 1 .. $half;
    my $product      = [*] ($arrangements, $flips, $permutations);
    return :$arrangements, :$flips, :$permutations, :$product;
}
my @tests =
    (7, 8, '05132231055052464303662006235126113002452143346664515414'),
    (7, 8, '64220650162341432102355113505610426040114420536366525334'),
    (3, 4, '001102221201'),
;
sub solution_works ( UInt $rows, UInt $cols, @tab, @sol --> Bool ) {
    die unless @sol.all eq <L R U D>.any;
    my BagHash $b .= new;
    for @sol.keys -> $p1 {
        given @sol[$p1] {
            when 'L' { my $p2 = $p1+1;     die unless @sol[$p2] eq 'R'; $b{ domino( @tab[$p1,$p2].join ) }++; }
            when 'D' { my $p2 = $p1+$cols; die unless @sol[$p2] eq 'U'; $b{ domino( @tab[$p1,$p2].join ) }++; }
        }
    }
    die unless $b.values.max == 1 and $b.elems == $rows * $cols / 2;
    return True;
}
for @tests -> ( $rows, $cols, $flat_table ) {
    my @solutions = solve($rows, $cols, $flat_table.comb);
    say 'Solutions found: ', +@solutions;

    # Confirm that every solution is valid.
    solution_works($rows, $cols, $flat_table.comb, $_) or die for @solutions;

    # Print first solution
    say .join.trans( [<L R D U>] => [<╼ ╾ ╽ ╿>]) for @solutions[0].batch($cols);

    if ++$ == 1 {
        say 'Bonus:';
        say .value.fmt("%45d "), .key.tclc for bonus($rows, $cols);
        say '';
    }
}
Output:
Solutions found: 1
╼╾╼╾╽╼╾╽
╽╼╾╽╿╼╾╿
╿╽╽╿╼╾╽╽
╽╿╿╼╾╽╿╿
╿╽╽╼╾╿╽╽
╽╿╿╼╾╽╿╿
╿╼╾╼╾╿╼╾
Bonus:
                                      1292697 Arrangements
                                    268435456 Flips
               304888344611713860501504000000 Permutations
105797996085635619122449324709940305294524416 Product

Solutions found: 2025
╼╾╼╾╼╾╽╽
╼╾╼╾╼╾╿╿
╼╾╼╾╼╾╼╾
╼╾╼╾╼╾╼╾
╼╾╽╽╼╾╼╾
╼╾╿╿╽╼╾╽
╼╾╼╾╿╼╾╿
Solutions found: 4
╼╾╼╾
╼╾╼╾
╼╾╼╾

When you have built-in path everything looks like a path-solving problem...

Try it in Uiua Pad!

GridOne ← ↯7_∞"05132231055052464303662006235126113002452143346664515414"
GridTwo ← ↯7_∞"64220650162341432102355113505610426040114420536366525334" # 2025 solutions

Offsets       ← [1_0 0_1]
FirstFreeCell ← ⊢⬚1▽⤚(¬°⊚⨂)⊙◌↯∞_2°⊡⊃(⋅∘|∘)
# Get pair of cells for down- and right- facing dominos.
CandDoms ← ≡⌞⊂⟜(≡˜⊟⊸+Offsets˙⊟ FirstFreeCell)
# Is the tail-cell of each domino in bounds and not visited?
ValidDoms ← ↧⊃(¬˜∊°⊂⇌|≠@_⬚@_⊡⊣)
# Has this domino already been 'played'?
NotDupeDom ← ⨬(¬˜∊°⊂≡⍆↯∞_2⇌⊡|⋅1) ⊸(=2⧻)
# Find two possible next doms, check them, and step.
Nexts ← ▽⤚≡NotDupeDom▽◡≡ValidDoms⊙¤ ⊸CandDoms
Done  ← =∩⧻⊙♭

# Add spaces to a solution grid, for MarkDominos to fill.
SpacedGrid  ← ↘1↯[∞ ∘]÷2⧻⊸⊢≡⊂˜↯@\s⊸△≡/(⊂$"_ ")
MarkDominos ← ∧(⍜⊡⋅@+)≡(/+↯∞_2) ↯∞_4

Solve ← (
  ⟜SpacedGrid
  °□⊣⊢path(≡□Nexts°□|Done°□)□[] # Show first solution.
  MarkDominos
)
&p"Only solution for first grid:"
&s Solve GridOne
Count ← ⧻path(≡□Nexts°□|Done°□)□[] # Count all solutions.
&p"Count of solutions for second grid:"
Count GridTwo
Output:
Only solution for first grid:
╭─                   
╷ "0+5 1+3 2 2+3 1"  
  "        +     +"  
  "0 5+5 0 5 2+4 6"  
  "+     +        "  
  "4 3 0 3 6+6 2 0"  
  "  + +       + +"  
  "0 6 2 3+5 1 2 6"  
  "+         +    "  
  "1 1 3 0+0 2 4 5"  
  "  + +       + +"  
  "2 1 4 3+3 4 6 6"  
  "+         +    "  
  "6 4+5 1+5 4 1+4"  
                    ╯
Count of solutions for second grid:
2025
Translation of: Julia

Basic task

var tableau = [
    [0, 5, 1, 3, 2, 2, 3, 1],
    [0, 5, 5, 0, 5, 2, 4, 6],
    [4, 3, 0, 3, 6, 6, 2, 0],
    [0, 6, 2, 3, 5, 1, 2, 6],
    [1, 1, 3, 0, 0, 2, 4, 5],
    [2, 1, 4, 3, 3, 4, 6, 6],
    [6, 4, 5, 1, 5, 4, 1, 4]
]

var tableau2 = [
    [6, 4, 2, 2, 0, 6, 5, 0],
    [1, 6, 2, 3, 4, 1, 4, 3],
    [2, 1, 0, 2, 3, 5, 5, 1],
    [1, 3, 5, 0, 5, 6, 1, 0],
    [4, 2, 6, 0, 4, 0, 1, 1],
    [4, 4, 2, 0, 5, 3, 6, 3],
    [6, 6, 5, 2, 5, 3, 3, 4]
]

var dominoes = []
for (j in 0...tableau[0].count) {
    for (i in 0...tableau.count) if (i <= j) dominoes.add([i, j])
}

var containsDom = Fn.new { |l, m, n|  // assumes m <= n
    for (i in 0...l.count) {
        var d = l[i]
        if (d[0] == m && d[1] == n) return true
    }
    return false
}

var copyTab = Fn.new { |t|
    var c = List.filled(t.count, null)
    for (r in 0...t.count) c[r] = t[r].toList
    return c
}

var sorted = Fn.new { |dom| (dom[0] > dom[1]) ? [dom[1], dom[0]] : dom }

var findLayouts = Fn.new { |tab, doms|
    var nrows = tab.count
    var ncols = tab[0].count
    var m = List.filled(nrows, null)
    for (i in 0...nrows) m[i] = List.filled(ncols, -1)
    var patterns = [ [m, [], []] ]
    var count = 0
    while (true) {
        var newpat = []
        for (pat in patterns) {
            var ut = pat[0]
            var ud = pat[1]
            var up = pat[2]
            var pos = null
            for (j in 0...ncols) {
                var breakOuter = false
                for (i in 0...nrows) {
                   if (ut[i][j] == -1) {
                       pos = [i, j]
                       breakOuter = true
                       break
                   }
                }
                if (breakOuter) break
            }
            if (!pos) continue
            var row = pos[0]
            var col = pos[1]
            if (row < nrows - 1 && ut[row+1][col] == -1) {
                var dom = sorted.call([tab[row][col], tab[row+1][col]])
                if (!containsDom.call(ud, dom[0], dom[1])) {
                    var newut = copyTab.call(ut)
                    newut[row][col] = tab[row][col]
                    newut[row+1][col] = tab[row+1][col]
                    newpat.add([newut, ud + [sorted.call( [tab[row][col], tab[row+1][col]])],
                        up + [row, col, row+1, col]])
                }
            }
            if (col < ncols - 1  && ut[row][col+1] == -1) {
                var dom = sorted.call([tab[row][col], tab[row][col+1]])
                if (!containsDom.call(ud, dom[0], dom[1])) {
                    var newut = copyTab.call(ut)
                    newut[row][col] = tab[row][col]
                    newut[row][col+1] = tab[row][col+1]
                    newpat.add([newut, ud + [sorted.call([tab[row][col], tab[row][col+1]])],
                        up + [row, col, row, col+1]])
                }
            }
        }
        if (newpat.count == 0) break
        patterns = newpat
        if (patterns[0][-1].count == doms.count) break
    }
    return patterns
}

var printLayout = Fn.new { |pattern|
    var tab = pattern[0]
    var dom = pattern[1]
    var pos = pattern[2]
    var bytes = List.filled(tab.count*2, null)
    for (i in 0...bytes.count) bytes[i] = List.filled(tab[0].count*2 - 1, " ")
    var idx = 0
    while (idx < pos.count-1) {
        var p = pos[idx..idx+3]
        var x1 = p[0]
        var y1 = p[1]
        var x2 = p[2]
        var y2 = p[3]
        var n1 = tab[x1][y1]
        var n2 = tab[x2][y2]
        bytes[x1*2][y1*2] = String.fromByte(48+n1)
        bytes[x2*2][y2*2] = String.fromByte(48+n2)
        if (x1 == x2) { // horizontal
            bytes[x1*2][y1*2 + 1] = "+"
        } else if (y1 == y2) { // vertical
            bytes[x1*2 + 1][y1*2] = "+"
        }
        idx = idx + 4
    }

    for (i in 0...bytes.count) {
        System.print(bytes[i].join())
    }
}

for (t in [tableau, tableau2]) {
    var start = System.clock
    var lays = findLayouts.call(t, dominoes)
    printLayout.call(lays[0])
    var lc = lays.count
    var pl = (lc > 1) ? "s" : ""
    var fo = (lc > 1) ? " (first one shown)" : ""
    System.print("%(lays.count) layout%(pl) found%(fo).")
    System.print("Took %(System.clock - start) seconds.\n")
}
Output:
0+5 1+3 2 2+3 1
        +     +
0 5+5 0 5 2+4 6
+     +        
4 3 0 3 6+6 2 0
  + +       + +
0 6 2 3+5 1 2 6
+         +    
1 1 3 0+0 2 4 5
  + +       + +
2 1 4 3+3 4 6 6
+         +    
6 4+5 1+5 4 1+4
               
1 layout found.
Took 0.014597 seconds.

6 4 2 2 0 6+5 0
+ + + + +     +
1 6 2 3 4 1+4 3
               
2 1 0 2 3+5 5 1
+ + + +     + +
1 3 5 0 5 6 1 0
        + +    
4 2 6 0 4 0 1+1
+ + + +        
4 4 2 0 5 3 6 3
        + + + +
6+6 5+2 5 3 3 4
               
2025 layouts found (first one shown).
Took 0.217176 seconds.

Extra credit (Cli)

Library: Wren-big
Library: Wren-fmt
import "./big" for BigInt
import "./fmt" for Fmt

var dominoTilingCount = Fn.new { |m, n|
    var prod = 1
    for (j in 1..(m/2).ceil) {
        for (k in 1..(n/2).ceil) {
            var cm = (Num.pi * (j / (m + 1))).cos
            var cn = (Num.pi * (k / (n + 1))).cos
            prod = prod * ((cm*cm + cn*cn) * 4)
        }
    }
    return prod.floor
}

var start  = System.clock
var arrang = dominoTilingCount.call(7, 8)
var perms  = BigInt.factorial(28)
var flips  = 2.pow(28)

Fmt.print("Arrangements ignoring values: $,i", arrang)
Fmt.print("Permutations of 28 dominos: $,i", perms)
Fmt.print("Permuted arrangements ignoring flipping dominos: $,i", perms * arrang)
Fmt.print("Possible flip configurations: $,i", flips)
Fmt.print("Possible permuted arrangements with flips: $,i", perms * flips * arrang)
System.print("\nTook %(System.clock - start) seconds.")
Output:
Arrangements ignoring values: 1,292,697
Permutations of 28 dominos: 304,888,344,611,713,860,501,504,000,000
Permuted arrangements ignoring flipping dominos: 394,128,248,414,528,672,328,712,716,288,000,000
Possible flip configurations: 268,435,456
Possible permuted arrangements with flips: 105,797,996,085,635,281,181,632,579,889,767,907,328,000,000

Took 0.00046 seconds.

Extra credit (Embedded)

Library: Wren-gmp

This is just to give what will probably be a rare outing to the Mpf class though (despite their usage in the Julia example) we don't need 'big floats' here, just 'big ints'. Slightly slower than the Wren-cli example as a result.

import "./gmp" for Mpz, Mpf
import "./fmt" for Fmt

var dominoTilingCount = Fn.new { |m, n|
    var prec = 128
    var prod = Mpf.from(1, prec)
    for (j in 1..(m/2).ceil) {
        for (k in 1..(n/2).ceil) {
            var cm = Mpf.pi(prec).mul(Mpf.from(j / (m + 1))).cos.square
            var cn = Mpf.pi(prec).mul(Mpf.from(k / (n + 1))).cos.square
            prod.mul(cm.add(cn).mul(4))
        }
    }
    return Mpz.from(prod.floor)
}

var start  = System.clock 
var arrang = dominoTilingCount.call(7, 8)
var perms  = Mpz.new().factorial(28)
var flips  = 2.pow(28)
 
Fmt.print("Arrangements ignoring values: $,i", arrang)
Fmt.print("Permutations of 28 dominos: $,i", perms)
Fmt.print("Permuted arrangements ignoring flipping dominos: $,i", perms * arrang)
Fmt.print("Possible flip configurations: $,i", flips)
Fmt.print("Possible permuted arrangements with flips: $,i", perms * flips * arrang)
System.print("\nTook %(System.clock - start) seconds.")
Output:
Arrangements ignoring values: 1,292,697
Permutations of 28 dominos: 304,888,344,611,713,860,501,504,000,000
Permuted arrangements ignoring flipping dominos: 394,128,248,414,528,672,328,712,716,288,000,000
Possible flip configurations: 268,435,456
Possible permuted arrangements with flips: 105,797,996,085,635,281,181,632,579,889,767,907,328,000,000

Took 0.00058 seconds.