Effective problem solving using SAT solvers
Abstract
In this article we demonstrate how to solve a variety of problems and puzzles using the builtin SAT solver of the computer algebra system Maple. Once the problems have been encoded into Boolean logic, solutions can be found (or shown to not exist) automatically, without the need to implement any search algorithm. In particular, we describe how to solve the queens problem, how to generate and solve Sudoku puzzles, how to solve logic puzzles like the Einstein riddle, how to solve the 15puzzle, how to solve the maximum clique problem, and finding GraecoLatin squares.
Keywords:
SAT solving Maple queens problem Sudoku Logic puzzles 15puzzle Maximum clique problem GraecoLatin squares1 Introduction
“…it is a constant source of annoyance when you come up with a clever special algorithm which then gets beaten by translation to SAT.”
—Chris Jefferson
The satisfiability (SAT) problem is to determine if a given Boolean expression can be satisfied—is there some way of assigning true and false to its variables that makes the whole formula true? Despite at first seeming disconnected from most of the kinds of problems that mathematicians care about we argue in this paper that it is in the interests of mathematicians to have a familiarity with SAT solving and encoding problems in SAT. An immense amount of effort over the past several decades has produced SAT solvers that are not only practical for many problems but are actually the fastest known way of solving an impressive variety of problems such as software and hardware verification problems [3]. They have also recently been used to resolve longstanding mathematical conjectures [12] and construct large combinatorial designs [7].
Since 2018, the computer algebra system Maple has included the awardwinning SAT solver MapleSAT [15] as its builtin SAT solver. This solver can be used through the Satisfy command of the Logic package. Satisfy returns a satisfying assignment of a given Boolean expression (if one exists) or NULL if no satisfying assignment exists. In this paper we demonstrate through a number of detailed examples how Satisfy can be an effective and efficient way of solving a variety of problems and puzzles.
Very little prerequisites are necessary to understand this paper; the main necessary background is a familiarity with Boolean logic which we outline in Section 2. We then present effective solutions to the queens problem (Section 3), logic puzzles like the Einstein riddle (Section 4), Sudoku puzzles (Section 5), Euler’s GraecoLatin square problem (Section 6), the maximum clique problem (Section 7), and the 15puzzle (Section 8). In each case we require no knowledge of any of the specialpurpose search algorithms that have been proposed to solve these problems; once the problems have been encoded into Boolean logic they are automatically solved using Maple’s Satisfy.
All of the examples discussed in this paper were implemented and run in Maple 2018 and Maple 2019. Due to space constraints we have not included our code in this paper, but Maple worksheets containing complete implementations have been made available online through the Maple Application Center [6].
2 Background
A basic understanding of Boolean logic is the only prerequisite necessary to understand the solutions described in this paper. One of the main advantages of Boolean logic (but also one of its main disadvantages) is its simplicity: each variable can assume only one of two values denoted by true and false. Boolean expressions consist of variables joined by Boolean operators. The most common Boolean operators (and the ones available in the Logic package of Maple) are summarized in Table 1.
Name  Symbol  Arity  Maple Syntax 

Negation  1  ¬  
Conjunction  ary  &and  
Disjunction  ary  &or  
Implication  2  &implies  
Biconditional  2  &iff  
Alternative denial  ary  &nand  
Joint denial  ary  &nor  
Exclusive disjunction  ary  &xor 
The (or), (and), and (not) operators have meanings based on their everyday English meanings: is true exactly when at least one is true, is true exactly when all are true, and is true exactly when is false. More generally, is true exactly when and have the same truth values, is false exactly when is true and is false, is true exactly when an odd number of are true, is true exactly when at least one is false, and is true exactly when all are false.
A literal is an expression of the form or where is a Boolean variable. A clause is an expression of the form where all are literals. A conjunctive normal form (CNF) expression is of the form where all are clauses. A standard theorem of Boolean logic is that any expression can be converted into an equivalent expression in conjunctive normal form where two expressions are said to be equivalent if they assume the same truth values under all variable assignments.
The current algorithms used in stateoftheart SAT solvers require that the input formula be given in conjunctive normal form. While this is convenient for the purposes of designing efficient solvers it is not convenient for the mathematician who wants to express their problem in Boolean logic—not all expressions are naturally expressed in conjunctive normal form. An advantage that Maple has over most current stateoftheart SAT solvers is that Maple does not require the input to be given in conjunctive normal form. Since the algorithms used by MapleSAT require CNF to work properly, Maple internally converts expressions into CNF automatically. This is done by using a number of equivalence transformations, e.g., the expression is rewritten as the clause .
Care has been taken to make the necessary conversion to CNF efficient. This is important because conversions that use the most straightforward equivalence rules generally require exponential time to complete. For example, the Maple command Normalize from the Logic package can be used to convert an expression into CNF. But watch out—many expressions explode in size following this conversion. For example, the expression when converted into CNF contains clauses. The main trick used to make the conversion into CNF efficient is the Tseitin transformation [23]. This transformation avoids the exponential blowup of the straightforward transformations by using additional variables to derive a new formula that is satisfiable if and only if the original formula is satisfiable. For example, the expression is rewritten as where is a new variable and is a CNF encoding of the formula , namely,
The transformation is then recursively applied to (the part of the formula not in CNF) until the entire formula is in CNF. The Maple command Tseitin of the Logic package can be applied to convert an arbitrary formula into CNF using this translation. Thus, Maple offers us the convenience of not requiring encodings to be in CNF while avoiding the inefficiencies associated with a totally unrestricted encoding.
3 The queens problem
The queens problem is to place chess queens on an chessboard such that no two queens are mutually attacking (i.e., in the same row, column, or diagonal). The problem was first proposed for by Bezzel in 1848 and the first solution for general was given by Pauls in 1874 [2]. The problem is solvable for all ; a solution for (found in 0.015 seconds using Satisfy) is shown in Figure 1.
The queens problem is a standard example of a constraint satisfaction problem [18]. The encoding that we use for this problem uses the Boolean variables with to denote if there is a queen on square . There are two kinds of constraints necessary for this problem: positive constraints that say that there are queens on the board and negative constraints that say that queens do not attack each other. A satisfying assignment of these constraints exists exactly when the queens problem is solvable.
Since there are rows and each row must contain a queen the positive constraints are of the form for . Similarly, each column must contain a queen; these constraints are of the form for . The negative constraints say that if contains a queen then all squares attacked by do not contain a queen. These constraints are represented in Boolean logic by where are the variables “attacked” by a queen on . In general this encoding uses constraints in order . Typically Satisfy is able to solve each order using slightly more time than the previous order and the last order it can solve in under a second is .
4 The Einstein riddle
The Einstein riddle is a logic puzzle apocryphally attributed to Albert Einstein and is often stated with the remark that it is only solvable by 2% of the world’s population. The true source of the puzzle is unknown, but a version of it appeared in the magazine Life International in 1962. In the puzzle there are five houses in a row with each house a different colour and each house owned by a man of a different nationality. Additionally, each of the owners have a different pet, prefer a different kind of drink, and smoke a different brand of cigarette. Furthermore, the following information is given:

The Brit lives in the red house.

The Swede keeps dogs as pets.

The Dane drinks tea.

The green house is next to the white house, on the left.

The owner of the green house drinks coffee.

The person who smokes Pall Mall rears birds.

The owner of the yellow house smokes Dunhill.

The man living in the centre house drinks milk.

The Norwegian lives in the first house.

The man who smokes Blends lives next to the one who keeps cats.

The man who keeps horses lives next to the man who smokes Dunhill.

The man who smokes Blue Master drinks beer.

The German smokes Prince.

The Norwegian lives next to the blue house.

The man who smokes Blends has a neighbour who drinks water.
The puzzle is: Who owns the fish?
To solve this riddle using Maple, we label the houses 1 to 5 and use the variables where and a is an attribute (one of the colours, nationalities, pets, drinks, or cigarette brands). For example, if is a colour then is in the set and similarly for the other attribute types; there are five distinct possible attributes for each type of attribute. In total there are possible values for and variables .
We know that each attribute is not shared among the five houses or their owners. Since there are exactly five houses, each attribute must appear exactly once among the five houses. The knowledge that each attribute appears at least once can be encoded as the clauses for each attribute and the knowledge that each attribute is not shared can be encoded as where is a house index not equal to and is an attribute. Additionally, the fact that each house has some colour is encoded as for each house index and the knowledge that each house cannot have two colours can be encoded as where is a house index and and are two distinct colours (and similarly for the other kinds of attributes).
The known facts can be encoded into logic fairly straightforwardly; for example, the first fact can be encoded into logic as for house indices and the last fact can be encoded as for and . Using Satisfy on these constraints produces the unique satisfying solution (that includes the equations and , thereby solving the puzzle) in under 0.01 seconds.
5 Sudoku puzzles
Sudoku is a popular puzzle that appears in many puzzle books and newspapers. Given a 9 by 9 grid whose squares are either blank or contain a number between 1 and 9, the objective is to fill in the blank squares in such a way that each row and column contains exactly one digit between 1 and 9. Additionally, each of the nine 3 by 3 subgrids which compose the grid (called blocks) must also contain exactly one digit between 1 and 9. Figure 2 contains a Sudoku puzzle designed by mathematician Arto Inkala and claimed to be the world’s hardest Sudoku [9].
It is known that Sudoku can be modelled as a SAT problem [16] or a constraint satisfaction problem [20]. A straightforward encoding uses variables with where is true exactly when the square contains the digit . The rules of Sudoku state that each square must be filled with a digit between 1 and 9 and that the same digit cannot appear twice in the same row, column, or block. The first constraint has the form for all and the second constraint has the form for all where does not equal but is in the same row, column, or block as .
One can also include the constraints for all with that say that each square can contain at most one digit. However, these constraints are unnecessary since they are logically implied by the first two constraints. In our tests, including these additional constraints slightly decreased the performance of Satisfy. Without the additional constraints the “world’s hardest Sudoku” was solved in 0.25 seconds and with them it was solved in 0.33 seconds.
Additionally, we developed a method of generating Sudoku puzzles with a unique solution using Satisfy. This allowed us to write an interactive Sudoku game where random puzzles can automatically be generated on command. To begin, Satisfy is used to find a solution to the above Sudoku constraints with an empty grid (no starting clues) and the produced satisfying solution generates a completed Sudoku grid . A random seed is passed to the SAT solver so that a different solution is generated each time; this is done by passing the following solveroptions in the call to Satisfy:
Let denote the the entry in the solution where . The are randomly ordered and the first 50 entries are selected as the potential starting configuration of a Sudoku puzzle. This puzzle has the solution by construction, though other solutions may also exist.
To verify that the generated solution is unique, we rerun Satisfy with the additional 50 unit clauses corresponding to the starting configuration along with the constraint which blocks the solution . If Satisfy returns another solution then we start over and find a new to try. Otherwise the starting configuration forms a legal Sudoku puzzle.
Additionally, it may be the case that we can use fewer than 50 entries and still obtain a Sudoku puzzle with a unique solution. To estimate how many entries need to be assigned using only a few extra calls to the SAT solver we use a variant of binary search, letting and be lower and upper bounds on how many entries we will define in the puzzle. Next, we let and repeat the first step except using only the first entries . If the resulting SAT instance is satisfiable then we need to use strictly more than entries to ensure that a unique solution exists and if the resulting SAT instance is unsatisfiable then we can perhaps use strictly fewer than entries. Either way, we improve the bounds on how many entries to assign (in the former case we can update to and in the latter case we can update to ) and this step can be repeated a few times to find more precise bounds on how many entries need to be assigned to ensure a unique solution exists.
6 Euler’s GraecoLatin square problem
A Latin square of order is an matrix containing integer entries between and such that every row and every column contains each entry exactly once. Two Latin squares are orthogonal if the superposition of one over the other produces all distinct pairs of integers between and . A pair of orthogonal Latin squares was called a GraecoLatin square by the mathematician Leonhard Euler who in 1782 used Latin characters to represent the entries of the first square and Greek characters to represent the entries of the second square [11]. Figure 3 contains a visual representation of a GraecoLatin square.
Euler studied the orders for which GraecoLatin squares exist and found methods for constructing them when was odd or a multiple of 4. Since such squares do not exist for and he was unable to find a solution for he conjectured that GraecoLatin squares do not exist when . Euler’s conjecture became famous as he was not able to resolve it in his lifetime.
The first progress on the conjecture did not come until over a hundred years later when in 1900 Tarry showed that GraecoLatin squares of order do not exist [22]. This gave credence to Euler’s conjecture and many mathematicians thought the conjecture was true—in fact, three independent proofs of the conjecture were published in the early 20th century [17, 19, 24]. In 1959–1960 Bose, Shrikhande, and Parker [4, 5] made explosive news (even appearing on the front page of the New York Times) by showing that these proofs were invalid by giving explicit constructions for GraecoLatin squares in all orders except two and six. As it turns out, a lot of time could have been saved if Euler had a copy of Maple—we now show that Euler’s conjecture can be automatically disproven in Maple. With Satisfy we are able to construct small GraecoLatin squares without any knowledge of search algorithms or construction methods.
Our encoding for the GraecoLatin square problem of order uses the variables and with . The variables will be true exactly when the th entry of the Latin square is and will be true exactly when the th entry of the Graeco square is .
There are three kinds of constraints that specify that is a GraecoLatin square: Those that specify that every entry of and is an integer between and , those that specify that the rows and columns of and contain no duplicate entries, and those that specify that and are orthogonal. Additionally, there are constraints that are not logically necessary but help cut down the search space. Some work has previously been done using SAT solvers to search for special kinds of GraecoLatin squares [25]. The encoding we use is similar but takes advantage of the fact that Maple does not require constraints to be specified in conjunctive normal form.
First, we specify that the entries of are welldefined, i.e., consist of a single integer between and . The constraints that say that each entry of contains at least one integer are of the form for each index pair and the constraints that say that each entry of contains at most one integer are of the form for each index pair and integer . Similar constraints are also used to specify that the entries of are welldefined.
Second, we specify that is a Latin square, i.e., all columns and rows contain distinct entries. These have the form where and but is in the same column or row as . Similarly, we also specify that is a Latin square.
Third, we specify that and are orthogonal, i.e., for every pair there exists some pair such that holds. These constraints are of the form for each pair .
Lastly, we include some “symmetry breaking” constraints. These constraints are not strictly necessary but they shrink the search space and thereby make the search more efficient. In general, when a search space splits into symmetric subspaces it is beneficial to add constraints that remove or “break” the symmetry. GraecoLatin squares have a number of symmetries, in particular, a row or column permutation simultaneously applied to and produces another GraecoLatin square. Also, any permutation of may be applied to the entries of either or .
The result of these symmetries is that any GraecoLatin square can be transformed into one where the first row and column of has entries in ascending order (by permuting rows/columns) and the first row of has entries in ascending order (by renaming the entries of ). Thus, we can assume the constraint . Altogether this encoding uses constraints.
Using this encoding the orders up to eight can be solved in 25 total seconds (including 14 seconds to show that no GraecoLatin squares exist in order six), a GraecoLatin square of order nine can be found in about 45 minutes, and a GraecoLatin square of order ten can be found in about 23 hours, thereby disproving Euler’s GraecoLatin square conjecture.
7 The maximum clique problem
The maximum clique problem is to find a clique of maximum size in a given graph. A clique of a graph is a subset of its vertices that are all mutually connected (see Figure 4). The decision version of this problem (does a graph contain a clique of size ?) is in NP, meaning that it is easy to verify the correctness of a solution if one can be found. By the Cook–Levin theorem [10] the problem can be encoded into a SAT instance in polynomial time. However, the reduction involves simulating the computation of a machine that solves the maximum clique problem and is therefore not very convenient to use in practice. Thus, we provide a simpler encoding into Boolean logic.
Suppose the given graph has vertices labelled , , and we want to find a clique of size in . Our encoding uses the variables , , where represents that the vertex appears in the clique we are attempting to find. We need to enforce a constraint that says that if and are true (for any distinct vertices ) then the edge exists in the graph . Equivalently, if the edge does not exist in the graph then the variables and cannot both be true (for any vertices ). In other words, for every edge in the complement of we use the clause .
Additionally, we need a way to enforce that the found clique is of size . The most naive way to encode this is as a disjunction over all conjunctions of length on the variables , , . However, this encoding is very inefficient in practice. A cleverer encoding uses Boolean counter variables (where and ) that represent that at least of the variables , , are assigned to true. We know that will be false for and that will be true for . Additionally, we know that is true exactly when is true or is true and is true. This is represented by the formulas
or in conjunctive normal form by the clauses , , , and . To enforce that the found clique contains at least vertices we also assign to true.
To solve the maximum clique problem for a given graph we initialize to (assuming the graph has at least one edge, otherwise the problem is trivial) and use the above encoding to search for a clique of size . If such a clique exists we increase by and repeat the search in this manner until a clique of size does not exist. The last satisfying assignment found then provides a maximum clique of using an encoding with variables and clauses.
Benchmark  SAT time (sec)  Maple 2018 (sec)  Vertices  Edges  Clique size 

brock200_2 
23.52  22.19  200  9876  12 
cfat2001 
0.71  0.03  200  1534  12 
cfat2002 
2.32  0.13  200  3235  24 
cfat2005 
9.07  20.63  200  8473  58 
cfat5001 
4.78  0.13  500  4459  14 
cfat5002 
11.12  0.82  500  9139  26 
cfat5005 
42.46  132.84  500  23191  64 
cfat50010 
134.42  Timeout  500  46627  126 
hamming62 
0.64  51.10  64  1824  32 
hamming64 
0.04  0.02  64  704  4 
hamming82 
58.51  Timeout  256  31616  128 
hamming84 
7.96  3393.26  256  20864  16 
johnson824 
0.01  0.01  28  210  4 
johnson844 
0.23  7.80  70  1855  14 
johnson1624 
5.62  642.24  120  5460  8 
keller4 
7.76  414.80  171  9435  11 
MANN_a9 
0.13  226.18  45  918  16 
p_hat3001 
11.81  3.30  300  10933  8 
p_hat5001 
308.87  34.60  500  31569  9 
p_hat7001 
1281.62  169.68  700  60999  11 
This implementation was tested on the maximum clique problems from the second DIMACS implementation challenge [13]. Additionally, it was compared with Maple’s MaximumClique function from the GraphTheory package that uses a branchandbound backtracking algorithm [14]. Of the 80 benchmarks, the SAT method solved 18 in under 3 minutes and the branchandbound method solved 15. The SAT approach was faster in over half of the solved benchmarks and in one case solved a benchmark in 8 seconds that MaximumClique required 57 minutes to solve (see Table 2).
The SAT method has been made available in Maple 2019 by using the method=sat option of MaximumClique. By default the MaximumClique function in Maple 2019 will run the previous method used by Maple and the SAT method in parallel and return the answer of whichever method finishes first. This hybrid approach is the method of choice especially when more than a single core is available.
8 The 15puzzle
The 15puzzle is a classic “sliding tile” puzzle that was first designed in 1880 and became very popular in the 1880s [21]. It consists of a grid containing tiles numbered through along with one missing tile (see Figure 5). The objective of the puzzle is to arrange the tiles so that they are in ascending order when read from left to right and top to bottom and to end with the blank tile in the lower right. The only moves allowed are those that slide a tile adjacent to the blank space into the blank space. Half of the possible starting positions are solvable [1] and the hardest legal starting positions require eighty moves to complete [8].
Our encoding of the 15puzzle is more complicated because unlike the other problems we’ve considered, a solution to the 15puzzle is not static. In other words, our encoding must be able to deal with the state of the puzzle changing over time. To do this we use the variables to denote that the entry at contains tile at timestep . Here , (we let denote the blank tile), and since each instance requires at most moves to complete.
The board does not permit two tiles to occupy the same location at the same time; these constraints are of the form for each tile numbers , valid indices and , and valid timesteps . Next we need to generate constraints that tell the SAT solver how the state of the board can change from time to time . There are two cases to consider, depending on if a square (or an adjacent square) contains the blank tile.
The easier case is when a square and none of the squares adjacent to that square are blank. In that case, the rules of the puzzle imply that the tile in square does not change. We define the function to be the constraint that says that the tile in square does not change at time ; these constraints are of the form . We also use the function to denote the squares adjacent to and define to be . The static transition constraints are of the form
for all valid squares and timesteps .
The harder transition case is when a square contains the blank tile. In this case we need to encode the fact that the blank tile will switch positions with exactly one of the squares adjacent to square . If the tile on square switches positions with the square at time this can be encoded as the constraint . We also need to enforce that all squares adjacent to other than do not change; these constraints are of the form where is adjacent to but not equal to . Let denote the conjunction of the above two constraints. Then the slide transition constraints are of the form
for all valid squares and timesteps .
The constraint that says the board is solved at timestep can be encoded as . For efficiency reasons we only start looking for solutions with at most 5 moves; if no solution is found then we look for solutions using at most 10 moves and continue in this manner until a solution is found. In other words, we call Satisfy with the given starting constraints, the constraints of the puzzle as described above, and the constraint where is initialized to and then increased by every time no solution is found.
This method was applied to the starting configuration from Figure 5. It found that no solutions with at most 5 moves exist in 1.5 seconds, no solutions with at most 10 moves exist in 2.8 seconds, and found a solution with 15 moves in 6.3 seconds. It was also able to solve puzzles requiring up to 40 moves in 20 minutes. While this is not competitive with dedicated solvers for the 15puzzle, it requires no knowledge beyond the rules of the game and makes an interesting example of how to push SAT solvers to their limits.
9 Conclusion
In this paper we’ve demonstrated how to solve a variety of problems and puzzles using the computer algebra system Maple and its SAT solver MapleSAT [15]. We discussed a number of encodings and ways for improving those encodings, e.g., by using symmetry breaking (as in Section 6) or by using auxiliary variables (as in Section 7). We also took advantage of Maple’s ability to solve SAT problems not encoded in conjunctive normal form in Sections 3, 6, and 8. Maple code for all the examples covered in this paper (including code to read the output of Satisfy and generate the figures included in this paper) are available for download from the Maple Application Center [6].
The implementations presented in this paper can be considered examples of declarative programming where the programmer focuses on describing the problem but not the solution—the computer automatically decides the best way to solve the problem. This is in contrast to imperative programming where a programmer needs to describe precisely how the computation is to take place. An advantage of declarative programming is that the programmer does not need to worry about specifying a potentially complicated search algorithm. However, a disadvantage of declarative programming is that it lacks the kind of detailed control over the method of solution that can be required for optimally efficient solutions. Furthermore, not all problems are naturally expressed in a declarative way.
As we saw in the maximum clique problem, sometimes declarative solutions can outperform imperative solutions. This also occurs in the graph colouring (or chromatic number) problem of colouring the vertices of a graph using the fewest number of colours subject to the constraint that adjacent vertices are coloured differently. For example, prior to Maple 2018 the ChromaticNumber function required several hours to find a minimal colouring of the queens graph but a SAT encoding can solve this problem in under 10 seconds [6]. The SAT approach is available in Maple 2019 using the method=sat option of ChromaticNumber.
This somewhat unconventional manner of using Maple is not applicable to all problems but we hope the examples in this paper have convinced the reader that SAT solvers are more useful and powerful than they might at first appear. With its extensive logic functionality and convenient method of expressing logical constraints, Maple is an ideal tool for experimenting with SAT solvers and logical programming.
References
 [1] Archer, A.F.: A modern treatment of the 15 puzzle. The American Mathematical Monthly 106(9), 793–799 (1999)
 [2] Bell, J., Stevens, B.: A survey of known results and research areas for queens. Discrete Mathematics 309(1), 1–31 (2009)
 [3] Biere, A., Heule, M.J.H., van Maaren, H., Walsh, T.: Handbook of satisfiability, Frontiers in Artificial Intelligence and Applications, vol. 185. IOS press (2009)
 [4] Bose, R.C., Shrikhande, S.S., Parker, E.T.: Further results on the construction of mutually orthogonal Latin squares and the falsity of Euler’s conjecture. Canadian Journal of Mathematics 12, 189–203 (1960)
 [5] Bose, R.C., Shrikhande, S.S.: On the falsity of Euler’s conjecture about the nonexistence of two orthogonal Latin squares of order . Proceedings of the National Academy of Sciences of the United States of America 45(5), 734–737 (1959)
 [6] Bright, C.: Maple applications by Curtis Bright. https://www.maplesoft.com/applications/Author.aspx?mid=345070
 [7] Bright, C., Kotsireas, I., Ganesh, V.: A SAT+CAS method for enumerating Williamson matrices of even order. In: McIlraith, S., Weinberger, K. (eds.) ThirtySecond AAAI Conference on Artificial Intelligence. pp. 6573–6580. AAAI Press (2018)
 [8] Brüngger, A., Marzetta, A., Fukuda, K., Nievergelt, J.: The parallel search bench ZRAM and its applications. Annals of Operations Research 90, 45–63 (1999)
 [9] Collins, N.: World’s hardest sudoku: can you crack it? https://www.telegraph.co.uk/news/science/sciencenews/9359579/Worldshardestsudokucanyoucrackit.html (2012)
 [10] Cook, S.A.: The complexity of theoremproving procedures. In: Proceedings of the third annual ACM symposium on Theory of computing. pp. 151–158. ACM (1971)
 [11] Euler, L.: Recherches sur un nouvelle espéce de quarrés magiques. Verhandelingen uitgegeven door het zeeuwsch Genootschap der Wetenschappen te Vlissingen pp. 85–239 (1782)
 [12] Heule, M.J.H., Kullmann, O., Marek, V.W.: Solving very hard problems: Cubeandconquer, a hybrid SAT solving method. In: Sierra, C. (ed.) Proceedings of the TwentySixth International Joint Conference on Artificial Intelligence, IJCAI17, pp. 4864–4868. IJCAI (2017)
 [13] Johnson, D.S., Trick, M.A.: Cliques, coloring, and satisfiability: second DIMACS implementation challenge, October 11–13, 1993, vol. 26. American Mathematical Soc. (1996)
 [14] Kreher, D., Stinson, D.: Combinatorial Algorithms: Generation, Enumeration, and Search. Discrete Mathematics and Its Applications, Taylor & Francis (1998)
 [15] Liang, J.H., Govind V. K., H., Poupart, P., Czarnecki, K., Ganesh, V.: An empirical study of branching heuristics through the lens of global learning rate. In: International Conference on Theory and Applications of Satisfiability Testing. pp. 119–135. Springer (2017), https://ece.uwaterloo.ca/maplesat/
 [16] Lynce, I., Ouaknine, J.: Sudoku as a SAT problem. In: 9th International Symposium on Artificial Intelligence and Mathematics (2006)
 [17] MacNeish, H.F.: Euler squares. Annals of Mathematics 23(2), 221–227 (1922)
 [18] Nadel, B.A.: Representation selection for constraint satisfaction: A case study using queens. IEEE Intelligent Systems 5(3), 16–23 (1990)
 [19] Peterson, J.: Les 36 officieurs. Annuaire des Mathématiciens pp. 413–427 (1902)
 [20] Russell, S.J., Norvig, P.: Artificial Intelligence: A Modern Approach. Pearson Education (2010)
 [21] Slocum, J., Sonneveld, D.: The 15 Puzzle: How It Drove the World Crazy. Slocum Puzzle Foundation (2006)
 [22] Tarry, G.: Le problème des 36 officiers. Association Française pour l’Avancement des Sciences: Compte Rendu de la 29 session en Paris 1900 pp. 170–203 (1901)
 [23] Tseitin, G.S.: On the complexity of derivation in propositional calculus. In: Slisenko, A.O. (ed.) Studies in Constructive Mathematics and Mathematical Logic. pp. 115–125 (1970)
 [24] Wernicke, P.: Das problem der 36 offiziere. Jahresbericht der Deutschen MathematikerVereinigung 19, 264–267 (1910)
 [25] Zaikin, O., Kochemazov, S.: The search for systems of diagonal Latin squares using the SAT@home project. International Journal of Open Information Technologies 3(11), 4–9 (2015)