Side length of unit cube, given one needs to use N unit cubes to tile/fill a 3D space of given dimension

Question

One has a rectangular-prism 3 dimensional space, e.g. {xmin, xmax}, {ymin, ymax}, {zmin, zmax}. The goal is to use a pre-defined "grid" of at most N cubes to maximally tile the 3 dimensional space. It will use no more than N cubes, if N cannot evenly be tiled into the space. Note cubes are located at the cube centerpoints (i.e. it is a grid), not the cube edges.

For example, in the simple case of N=30, and {0,3}, {0,3}, {0,3}. The correct answer is to use cubes of side length 1, in which case 27 cubes will be used in 3 rows of 3 columns of 3 stacks.

In the more difficult case of N=1000, and {0, 100}, {0, 0}, {0, 50}, obviously the y-dimension uses only 1 row of cubes (it can not use any more), so the z and x dimensions must use the 1000 cubes arranged at a ratio of 2:1. 100/L * 50/L = 1000. So, side length is probably around 2.2361, because 22 x 1 x 45 = 990 is the closest we're going to get to 1000 (alternatively, we've minimized the expanse that we could not cover with our cube centers). Of course the sides must be an integer number of the cubes.

---

I'm trying to work out a system of equations to solve this, but I have a weird feeling that this is a satisfaction (minimization) problem and thus there needs to be some iteration... Here is an attempt at the system of equations anyway:

L, xe, ye, ze ∊ Positive Reals
N, x, y, z ∊ Natural Numbers
x*y*z ≤ N
(x-1)*L ≤ xe
(y-1)*L ≤ ye
(z-1)*L ≤ ze
x,y,z ≥ 1
L = min(i) { (xe*ye*ze) - ((x-1)*i*(y-1)*i*(z-1)*i) }

I made it so that it is x-1, y-1 etc., to mean that they can go outside the width by exactly and at most 1 cube size. I can't decide whether it is correct to do this in all dimensions, or specify that it must go outside in at most one dimension (the smallest?).

Answer 1

General idea

The intuition here is that the amount of cubes in each direction is going to be proportional to the size of the space in that direction, because cubes are regular. Another hint to this is that the solution to the problem "how many cubes in each direction" does only depend on the ratios between sizes, not on their actual values: if you multiply each dimension by the same dilating factor, your solution will still be the same.

Simplifying the problem

Suppose your space is of size dx, dy, dz, all non-zero, and you want to use up to N cubes.

Let us express our constraints as:

find Nx, Ny, Nz, L such that
- Nx * Ny * Nz ≤ N
- Nx * L ≤ dx
- Ny * L ≤ dy
- Nz * L ≤ dz
Minimizing the "lost volume": (dx * dy * dz - L^3 * Nx * Ny * Nz)

You may want to maximize Nx * Ny * Nz instead. Note also that minimizing the volume that isn't tiled is equivalent to maximizing L^3 * Nx * Ny * Nz.

A final observation is that if your solution is optimal in any way (cube number or tiled volume), at least one of the dimensions will have the equality N@ * L == d@, otherwise we may increase L.

Let us suppose without loss of generality that dx is this dimension and let us divide all lengths by dx, thus with a@ = d@ / dx for @ being y and z (and ax = 1 if needed). Then we may rewrite the problem:

find Nx, Ny, Nz, L such that
- Nx * Ny * Nz ≤ N
- Nx = dx / L
- Ny / Nx ≤ ay
- Nz / Nx ≤ az
Maximizing the cubes' volume: dx^3 * (Ny / Nx) * (Nz / Nx)

Hence your problem is how to get ratios of numbers of cubes on each side as close as possible to the ratios of the sizes, while maintaining Nx * Ny * Nz ≤ N -- the same way you instinctively approached your N = 1000 case.

Algorithm suggestion: You could, for example, iteratively increase Nx and set N@ = floor(a@ * Nx) for @ in {y,z} while Nx * Ny * Nz ≤ N, keeping the best solution.

Optimal solutions

Note that in the case that your ratios ay and az are rational, as in your examples, an exact solution to the minimization might be found rather easily.

Let us write a@ = n@ / m@ for @ in {y,z}, then we can find an exact solution with by setting
Nx = lcm(my, mz) - the least common multiple of the denominators. Note that in the case where dx, dy, dz are all integer, this means Nx = dx / gcd(dx, dy, dz).

Thus if Nx^3 * ax * ay = lcm(my, mz) ^3 * ny * nz / (my * mz) ≤ N, this solution also fits the constraint on N. To write it out in full:

Nx = lcm(my, mz)
L  = dx / Nx
Ny = Nx * ay
Nz = Nx * az

In your examples, for the 3D case, you have ay = az = 1 and a solution with Nx = Ny = Nz = 1, and for the 2D case we have only dimensions x, z with az = 1/2, this you have a solution with Nx = 2 and Nz = 1

Now you can add more cubes because you have a lot of slack in the Nx * Ny * Nz ≤ N constraint. In particular, you can multiply each coordinate by floor(cbrt(N / (Nx * Ny * Nz))) which is floor(cbrt(30 / 1)) = 3 in your first case, and floor(sqrt(1000 / 2)) = 22 in your 2D case.

Summing up results to your examples:

space={3,3,3} N=30 : Nx=3, Ny=3, Nz=3, L=1
space={100, 50} N=1000 : Nx=44, Nz=22, L=2.72727272...

Heuristic(s)

If we can't find an exact answer, we can make an educated guess. In order to not iterate over Nx values, or just to have a sensible initial guess.

Note that the bigger the denominator, the smaller the intervals between two fractions with this same denominator. This doesn't mean that given a real number your want to approximate, the precision of your approximation is proportional to the size of the denominator, but that the average precision when approximating a lot of real numbers would be better (hence, that you've got a better probability in the single case).

We are trying to maximize (Ny / Nx) * (Nz / Nx), thus from our constraints on ay and az we get (N@ + 1) / Nx ≥ a@, hence N@ ≥ a@ * Nx - 1. Now, Nx * Ny * Nz ≤ N thus a bound on Nx is:
Nx * (Nx * ay - 1) * (Nx * az - 1) ≤ N

But I'm too lazy to solve that so let's even simplify some more: N@ ≥ a@ * Nx - 1 ≥ a@ * Nx
Thus:Nx * (Nx * ay) * (Nx * az) ≤ N

So you may try as a first guess:

Nx = floor(crbt(N / (ay * az)))
L  = dx / Nx
Ny = floor(Nx * ay)
Nz = floor(Nx * az)

or (exactly the same thing, without the a@):

Nx = floor(crbt(N * dx * dx / (dy * dz)))
L  = dx / Nx
Ny = floor(dy / L)
Nz = floor(dz / L)

If you want to make sure you have the optimal value, you might want to try to see if the simplifications restrained Nx too much, so increment it until Nx * Ny * Nz ≥ N and keep the best configuration.

Similarly, different values of Nx might allow to get Nx * ay and Nx * az closer to their integer part, thus wasting less space. You could decrement Nx until it is small. Note that you can skip any value that divides a value of Nx that you already tried.

Choice of "filled dimension" x

Remains the choice of which direction will be x -- instinctively I would say take the dimension in which your space is the biggest but this is again only a heuristic, and if you want the optimal answer you will have to try all directions.

On exact solutions though, where side size ratios are rational, picking any side as x will do since all dimensions will be filled.