Rigorous packing of unit squares into a circle

This paper considers the task of finding the smallest circle into which one can pack a fixed number of non-overlapping unit squares that are free to rotate. Due to the rotation angles, the packing of unit squares into a container is considerably harder to solve than their circle packing counterparts. Therefore, optimal arrangements were so far proved to be optimal only for one or two unit squares. By a computer-assisted method based on interval arithmetic techniques, we solve the case of three squares and find rigorous enclosures for every optimal arrangement of this problem. We model the relation between the squares and the circle as a constraint satisfaction problem (CSP) and found every box that may contain a solution inside a given upper bound of the radius. Due to symmetries in the search domain, general purpose interval methods are far too slow to solve the CSP directly. To overcome this difficulty, we split the problem into a set of subproblems by systematically adding constraints to the center of each square. Our proof requires the solution of 6, 43 and 12 subproblems with 1, 2 and 3 unit squares respectively. In principle, the method proposed in this paper generalizes to any number of squares. Electronic supplementary material The online version of this article (10.1007/s10898-018-0711-5) contains supplementary material, which is available to authorized users.

). The proof relies on geometrical arguments and not on rigorous computations. Recent contributions in the packing of unit squares into a square include new bounds for the wasted area [6], the optimality proof for the cases n = 5, . . . , 10, 13 and 46 [1,10,23] and the optimality proof for n − 2 and n − 1 whenever n is a square [19]. Again, none of these contributions rely on computer-assisted proofs. For a dynamic survey on the packing of unit squares in a square, see [10].
The packing of unit squares into general containers received considerably less attention than the circle or the unit square packing into a square. For example, Friedman [9] maintains a list of proved and best-known values for the packing of unit squares into circles, triangles, L-shapes, and pentagons. In each case, only trivial arrangements are proved optimal. For the subject of interest in this paper, the packing of unit squares into a circle, the first open case is n = 3. For a list of figures of squares packed into a circle, see https://www2.stetson.edu/ efriedma/squincir/.

Contribution and outline
This paper introduces a computer-assisted method for finding rigorous enclosures for r in Problem (1) and the corresponding optimal arrangements. The method is of theoretical interest since it proves optimality instead of only presenting a feasible arrangement. Therefore, it is suitable for small values of n only. Our approach relies on the interval branch-and-bound framework. We implement the algorithm in C++ using the forward-backward constraint propagation [21] to reduce the search domain. Section 2 introduces the solver. The code is available at http://www.mat. univie.ac.at/~montanhe/publications/n3.zip. Section 3 formulates Problem (1) as a constraint satisfaction problem (CSP). This paper uses the concept of sentinels [4,18] to model non-overlapping conditions and the convexity of the circle to write containment constraints. Given an upper bound r n for r n , the CSP asks for every feasible arrangement satisfying r ≤ r n . Our software produces a list of small interval vectors with the property that every optimal arrangement of (1) belongs to at least one element in the list.
General purpose interval solvers are usually not capable of solving packing problems due to symmetries in the search domain. To overcome this difficulty, Sect. 4 shows how to split the original CSP into a set of subproblems by systematically adding constraints to the center of each square. We call them tiling constraints as the idea resembles the one proposed in [15,16,24]. The tiling divides the search domain into a set of isosceles triangles that must contain the center of at most one unit square. Then, one can replace the original CSP by a set of K n subproblems, where K is the number of triangles in the tiling. Our procedure iterates on the number of squares to avoid the exponential growth of subproblems. At the i-th iteration, we look at every possible combination of i triangles which can accommodate i unit squares into a circle with the radius at most r n . The rationale behind this strategy is twofold: (i) It allows us to discard a large number of hard subproblems by proving the infeasibility of more straightforward cases and (ii) It propagates the reduction on the search domain through the iterations. We also show that some combinations of triangles are symmetric by construction. Then one can discard them without any processing. This observation in addition to our iterative method reduces the number of hard cases considerably.
Section 5 illustrates the capabilities of our method. We find a mathematically rigorous enclosure for r 3 and the corresponding optimal arrangement. If one set r = 5 √ 17 16 as pointed by Friedman [9], the tiling produces 36 triangles. Our approach requires the solution of 6 subproblems with one square, 43 with two and only 12 subproblems with 3 squares to conclude the proof. It is < 1% of all possible 36 3 = 7140 combinations. The method could also be used to find optimal configurations for higher values of n (e.g., n = 4, 5, 6).

Interval notation
This paper is an application of the interval branch-and-bound framework [11,13]. We assume that the reader is familiar with concepts from interval analysis [20]. Let a, a ∈ R with a ≤ a. The set of nonempty compact real intervals is given by Let S ⊆ R be any set. Then the interval hull S of S is the smallest interval containing S. An interval vector (also called box) x := [x, x] is the Cartesian product of the closed real intervals We denote the set of all interval vectors of dimension n by IR n . We apply the width operator component wise on vectors. Therefore max(wid(x)) := max(wid(x 1 ), . . . , wid(x n )). Interval operations and functions are defined as in [13,20]. The absolute value of the interval a is given by Let a and b be two intervals. The maximum of a and b is defined by Let F : R n → R m be a function defined on x ∈ IR n and let f ∈ IR m . We denote the natural interval extension of the function F by F. A constraint satisfaction problem (CSP) is the task of finding every point satisfying We call x the search domain and the problem is said to be infeasible if there is no x ∈ x satisfying F(x) ∈ f. We also denote constraint satisfaction problems by the triplet (F, f, x).

The algorithm
This section describes the algorithm designed for solving the subproblems of form (1) using interval arithmetic [11,13,20] The solver consists of two components, the memory, and the reducer. The former manages the branch-and-bound tree while the latter is responsible for processing the current box. There is also a post-processing step called cluster builder to group close boxes in the solution list. The memory keeps the list of unprocessed boxes. It is also responsible for the box selector, and to split the boxes coming from the reducer that cannot be discarded or saved as a solution. In this paper, the selector is a depth-first search procedure while the splitter creates two boxes by dividing the input in the midpoint of the coordinate with maximum width.
The reducer contains a list of rigorous methods to reduce or discard boxes. This paper uses the forward-backward constraint propagation [21] and a feasibility verification method [7]. We consider a CSP of form (F, x, f) in the next paragraphs to overview each method.
The forward-backward constraint propagation decomposes F into a set of simple functions (like the exponential function or the sum of several elements) and displays the pieces in a graph. The forward step is a procedure to evaluate F(x) systematically. In this case, the data flow from the decision variable nodes of the graph to the constraint nodes F 1 , . . . , F m . At the end of this step, each constraint node contains an enclosure of F i (x) ∩ f i . The backward step acts reversely. It starts from the constraint nodes F(x) ∩ f and walks the graph applying inverse functions until reaching x 1 , . . . , x n . At the end of the backward step, we have a new box x ⊆ x with the reduced search domain. This paper employs the following feasibility verification method. Let x be a box and define the midpoint of x as x * . Then, we build a small box x * around the x * and check its feasibility. The box x * is a feasible if F(x * ) ⊆ f. We also save a box x as solution if it satisfies max(wid(x)) < x for a given x > 0.
The order into which we call the rigorous methods to process x may influence the efficiency of the branch-and-bound procedure. In this paper, the methods follow the finite state machine described in Table 1.
The parameter T > 0 is the threshold tolerance which controls the relative gain of the box x ⊆ x with the help of the following function It is clear that the input of G Rel at each iteration is the box x and the outcome of the rigorous method, x . After processing every box in the memory, we run a post-processing step to build clusters of solutions. This method supports the analysis of the solution list since it reduces the number of boxes on it. Given two intervals a and b, we define the gap between a and b by We save two boxes x, y ∈ IR n in the same group if where C is the cluster builder tolerance. After assigning a group to every box in the solution set, we return the interval hull of each group and conclude the procedure. Algorithm 1 summarizes the interval branch-and-bound method. We implement the algorithm in C++ using two interval arithmetic libraries, the Filib [14] and the Moore [17]. The user can choose any of these implementations in the verification of the proof. We report only results from the test with Filib in this paper. The supplementary material also reports the results utilizing Moore. They are consistent with each other.

Algorithm 1 Simplified solver
Input: The CSP (F, f, x). The acceptance, threshold and cluster builder tolerances x > 0, T > 0 and C > 0 respectively. The list M of rigorous methods to reduce the search space. Output: The list L of boxes satisfying the following property. If x ∈ x is a feasible point of (F, f, x) then there exists at least one box x ∈ L such that x ∈ x . 1: Run the forward-backward constraint propagation procedure on (F, f, x). Save the reduced domain in x ; 2: if x is proved to be infeasible then 3: Return ∅; 4: end if 5: Start the list of boxes to be processed with x ; 6: L ← ∅; 7: while has boxes to process do 8: Run the problem selector to obtain x; 9: i ← 1; 10: while i ≤ |M| do 11: Run the strategy M i on (F, f, x) to obtain the reduced box x ; 12: if G Rel (x, x ) > T then 13: i ← 1; 14: x ← x ; 15: continue; 16: end if 17: i ← i + 1; 18: x ← x ; 19: end while 20: if wid(x) ≤ x then 21: L ← x; 22: continue; 23: end if 24: Run the splitter on x and stack all subproblems in the memory; 25: end while 26: Run the cluster builder on L with the cluster tolerance C ; 27: Return L;

The standard model
This section introduces the mathematical model for the containment and the non-overlapping conditions of (1). We call the resulting model the standard constraint satisfaction problem since it is the same for every subproblem. We assume that the squares have side length s. The inequalities for the containment condition follow from the convexity of the circle. On the other hand, non-overlapping constraints rely on the concept of sentinels [4,18].

Containment
Let C r be the closed circle of radius r and centered at the origin. The convexity of the circle implies that c ∈ C r for any point c in the segment of line ab if a, b ∈ C r . Then, a given square belongs to C r if and only if its vertices belong to C r .
Let S 0,0 be the open square centered at the origin, with no rotation angle and side length s. Then We denote the closure of a set S by S. The set of vertices of S 0,0 is given by For any c ∈ R 2 and θ ∈ R, we define the displacement operator as where A θ is the rotation matrix The open square centered at c ∈ R 2 , with rotation angle θ ∈ [0, π 2 ) and side length s is the set given by The set of vertices of S c,θ , denoted by V c,θ , is the union of the following points Finally, we denote the circle of radius r and centered at the origin by C r . Then Proposition 1 Let g r (x) := x 2 1 + x 2 2 − r 2 and consider the following inequalities Then S c,θ ⊆ C r ⇔ (4) hold.  (4) hold. Conversely, since S c,θ is a bounded polytope, it is given by the convex hull of the elements of V c,θ . The result follows from the convexity of the circle.

Non-overlapping
This subsection shows that two squares S c 1 ,θ 1 and S c 2 ,θ 2 are non-overlapping if and only a set of nine points defined on S c 1 ,θ 1 do not belong to S c 2 ,θ 2 and vice-versa. We call such sets sentinels of a square. Figure 1 illustrates the need of the sentinels in the non-overlapping formulation.
The set of sentinels of S 0,0 is given by We denote the set of sentinels of S c,θ by T c,θ . This set is given by the union of the following points The next theorem states that the non-overlapping condition between two squares reduces to the containment verification of their sets of sentinels. It is a particular case of the sentinels theorem proved in [18].

Theorem 1 Let S c i ,θ i and S c j ,θ j be two squares defined by (3) and let T c i ,θ i and T c j ,θ j be their corresponding sets of sentinels. Then
In order to check conditions of form S c i ,θ i ∩ T c j ,θ j = ∅ numerically, we need the definition of the inverse of the displacement operator (2) Lemma 1 Let z ∈ R 2 and S c,θ be a square defined by (3). Then where h −1 1 and h −1 2 are the coordinates of the inverse operator. Proof If z ∈ S c,θ then there exists x ∈ S 0,0 such that x = h −1 (c, θ, z) and the implication follows immediately. Conversely, let x := h −1 (c, θ, z). The left hand side of the equivalence implies that x ∈ S 0,0 . If we let z := c + A θ x then z = c + A θ A T θ (z − c) = z. Therefore z ∈ S c,θ and the result follows.
Applying the inverse of the displacement operator of the square S c i ,θ i to the point T P Let c j,1 and c j,2 be the coordinates of the vector c j . Then the coordinates of (5) are given by The following proposition shows that the verification of S c i ,θ i ∩ T c j ,θ j = ∅ reduces to the evaluation of nine non-smooth functions.

Proposition 2 Let S c i ,θ i and T c j ,θ j be as in Theorem 1 and define the function
Proof Follows from the application of the Lemma 1 to the elements of T c j ,θ j .

The standard model
We conclude this section with the formal statement of the standard constraint satisfaction problem. Here and throughout we assume, without loss of generality, that the angle of the first square is always 0. This condition follows from the proper rotation of the remaining squares into the circle.

Definition 1 [SCSP]
Let r > 0 be an upper bound for the radius of the smallest circle into which one can pack n non-overlapping unit squares and s be a scaling factor. We denote the following problem by standard constraint satisfaction problem (SCSP) find (r , c 1 , θ 1 , . . . , c n , θ n ) Functions g r and u are given by Propositions 1 and 2 respectively.

Tiling
General purpose interval branch-and-bound procedures cannot solve the SCSP in a reasonable amount of time even for small values of n due to symmetries in the search space. This section introduces a tiling method to split (6) into a set of subproblems suitable for the Algorithm 1.
We employ the Matlab-like notation g := a : s : b to denote the array with k := b−a s +1 elements where g i := a + is for i = 0, . . . , k − 1. In addition, we denote the array with the midpoints of g by g c . Then, Let r > 0 be an upper bound for the SCSP. Then, the step length In the same way, we write the elements of C Let ABC be the triangle with vertices A, B, C ∈ R 2 . Then, we define the following triangles for 0 Here, T , L, D and R stand for top, left, down and right respectively. Figure 2 shows that the definition aims to split the square with vertices v i, j , v i+1, j , v i+1, j+1 and v i, j+1 into four triangles. One can easily verify that the triangles can be written as

For
The proof is similar for o ∈ {L, D, R}. 6. Let S i, j be the closed square with vertices The result follows by noting that We also assign a label to each triangle in the tiling. It helps us to easily identify a specific triangle during the proof of the case n = 3 in Sect. 5. Triangles of form T i, j receive an index that is divisible by 4. In the same way, we assign labels to the left, down and right triangles with the congruence classes 1, 2 and 3 modulo 4, respectively. We denote the triangle with label i by T i . Figure 3-Left shows the tiling for the best known upper bound of r 3 .
We show now that each triangle of form o i, j contains the center of at most one unit square.

Lemma 3 The minimal distance between the centers of two non-overlapping unit squares is 1.
Proof Assume the contrary, let pq be a line segment of the centers with lower than 1. Let C p and C q the circles of radius 1 2 drawn into the squares. Then C p and C q intersect. But then since the squares are supersets of C p and C q , respectively, they also intersect. A contradiction.
Proof Lemma 2-1 shows that l < 1 and Lemma 2-5 gives that the base length of o i, j is l while its legs have length l √ 2 2 . The result follows from Lemma 3.
Let K := 4 p 2 be the number of triangles in the tiling. Proposition 3 states that we can split the SCSP into a set of K n subproblems. In each subproblem, we enforce that the center of each square belongs to a given triangle. For example, one can define the subproblem T 0 T 19 T 33 in the same tiling displayed in Fig. 3-Left. In this case, we set the standard constraint satisfaction problem defined in (1) and add to the model the linear inequalities given by Eqs. (8)-(11) for T 0,0 , R 1,1 and L 2,2 respectively. We conclude this subsection by proving that several subproblems can be discarded without any processing due to symmetries in V and C. Proposition 4 allows us to discard subproblems that are symmetric by rotations or reflections. For example, let r 3 = 5 √ 17 16 and r 3 , S c 1 ,θ 1 , S c 2 ,θ 2 , S c 3 ,θ 3 be a feasible arrangement for (6) with c 1 ∈ T 7 , c 2 ∈ T 12 and c 3 ∈ T 22 . Then, Proposition 4 ensures that there exists a feasible arrangement r 3 , S c 1 ,θ 1 , S c 2 ,θ 2 , S c 3 ,θ 3 satisfying c 1 ∈ T 19 , c 2 ∈ T 12 and c 3 ∈ T 22 . Moreover, since T 19 T 12 T 22 is obtained by a reflection around the y axis of T 7 T 12 T 22 , we know that c i = f y (c i ) for i = 1, 2, 3.
The tiling produced by Algorithm 2 suffices if one wants to use a complete global optimization approach for the packing problem. On the other hand, it is not suitable for a rigorous approach since the elements in V and C are floating point vectors subject to rounding errors. To overcome this problem, we introduce a scaled tiling. In this case, we ensure that the points at V and C are integer vectors to the cost of working with squares that are not unit but have the side length contained in a small interval s. Algorithm 3 produces the scaled vertices for the tiling as well as the interval s.
The elements in V and C are integer vectors by construction. Then, the Eqs. Markót and Csendes [24] propose tiling constraints for the circle packing problem based on rectangles. The same idea could be used for the packing of squares into a circle. On the other hand for the case n = 3, one would need to split the search domain in 144 squares instead of 36 as proposed in this paper.

Packing 3 unit squares
Friedman [9] gives an upper bound for the case n = 3, r 3 = 5 Moreover, the parameters of S c 1 ,θ 1 , S c 2 ,θ 2 and S c 3 ,θ 3 belong to the boxes in Table 6.
Proof We perform the computational part of the proof in a core i7 processor with a frequency of 2.6 GHz, 6 GB of RAM and Windows 10. We compiled the code using the g++ 7.3 compiler with the option −O3. A supplementary material for the proof, containing the statistics and the log files for each subproblem is available in http://www.mat.univie.ac.at/~montanhe/ publications/n3.zip We prove the theorem in three phases. At the i-th iteration, we consider instances of form (6) and define subproblems by adding tiling constraints of form (8)-(11) accordingly.
The proof considers the scaled version of the problem to ensure the mathematical certainty of our statements. Therefore, the CSPs in this section are of form (6) with the constant s replaced by the interval s in Eq. (12). We obtain the unscaled interval for r 3 and Table 6 by dividing every box found in the last iteration by s.
We also assume the labeling scheme for the triangles introduced in Sect. 4 and displayed on Fig. 3-Right. Therefore, the subproblem T 7 T 12 T 22 refers to the SCSP with the interval scaling factor s and such that c 1 ∈ T 7 := R 0,1 , c 2 ∈ T 12 := T 1,0 and c 3 ∈ T 22 := D 1,2 . Phase 1 In this iteration, we are interested in reducing the search domain of each subproblem and finding triangles which can contain the center of squares with no rotation. The   Since we are assuming that θ 1 = 0 in the optimal configuration for n = 3, we only have to check the combinations containing at least one of these triangles. Phase 2 This phase aims to discard as many instances as possible to reduce the number of hard subproblems in the last iteration. There are 630 possible combinations of 36 triangles taken 2 by 2. After eliminating symmetric and previously discarded cases, we obtain 43 instances. We also propagate any reduction in the search domain in the first phase to the subproblems in the second phase. Again, we remove the condition θ 1 = 0 from Problem (6).
We run the Algorithm 1 with T = 10 −1 , C = 10 −11 , x = 10 −13 and time limit of 3600 s. We stop the algorithm as soon as the feasibility verification method described in Sect. 2 succeeds in finding a feasible point. The supplementary material contains the list of all instances discarded without processing. Table 4 gives the statistics for the 43 processed instances.
We conclude the second phase with 22 infeasible subproblems. Again, any case in the next phase containing a combination found infeasible in this step can be discarded without any processing.
Phase 3 In this phase we set the full model in Problem (6), including the constraint θ 1 = 0. Table 3 shows that c 1 ∈ T 7 or c 1 ∈ T 16 . Therefore, after removing the cases where one of these conditions do not hold and eliminating symmetric and already proved infeasible subproblems, we obtain 12 instances of the 7140 possible ones.
If an instance contains both triangles T 7 and T 16 , we denote by T 7 * T 16 T x the case where we enforce the angle of the square centered in T 7 to be zero. In the same way, we write T 7 T 16 * T x for the instances where the square centered in T 16 has no rotation angle.
For the last phase, we run Algorithm 1 with T = 10 −1 , C = 10 −11 , x = 10 −13 and no time limit. Table 5 provides the statistics of the processed instances. Moreover, Table 5 shows that it is the only instance containing the optimal configurations for n = 3. Figure 4-Right shows an approximation of the center of each square in the optimal case. Algorithm 1 produces 4 clusters for the instance T 7 T 12 T 22 . The maximum width of a cluster is 6.23 * 10 −13 . The precision is smaller than x due to the cluster builder procedure described in Sect. 2. Table 6 gives the unscaled clusters.

Conclusion
This paper presents a framework for the rigorous optimization of the packing of unit squares into a circle. We express the question as the standard constraint satisfaction problem stated Table 6 Enclosures of the optimal arrangement for n by Definition 1. The model considers the concept of sentinels to formulate non-overlapping constraints and the convexity of the squares and the circle to describe containment conditions. General purpose rigorous optimization solvers cannot achieve the solution of the standard constraint satisfaction problem due to symmetries in the search domain. To overcome this difficulty, we propose a tiling method that splits the search space related to the center of each unit square into isosceles triangles. Our tiling divides the original problem into a set of subproblems that are suitable for the interval branch-and-bound approach. We also ensure that the parameters in each subproblem are free of rounding errors by introducing a proper scaling of the search domain.
To show the capabilities of our approach, we solve the first open case reported in the literature, n = 3. We implement the interval branch-and-bound in the C++ and the code is publicly available. We perform the proof on an ordinary laptop with 6 GB of RAM and a core i7 processor.
The proof of the case n = 3 requires the solution of 6 subproblems with one square, 43 with two and only 12 with three squares. We discard most subproblems without processing due to symmetries in the tiling. Among the 61 subproblems, just 6 require more than 100 s to conclude the search. At the end of the process, we obtained 4 boxes with the following properties 1. The maximum width of any coordinate of the resulting boxes is 6.23 * 10 −13 . 2. If one disregard symmetries, every solution of (1) is contained in at least one of the 4 boxes The method proposed in this paper could, in principle, be used to find the optimal arrangement for higher values of n (e.g., n = 4, 5, 6.).