Project Euler Solution 85: Counting Rectangles

In Problem 85: Counting Rectangles we are asked to compute the number of rectangles on a grid.

As you can see in the image in the problem statement, there are many different ways to put a rectangle onto a finite grid of points. For a grid of size 3×2 there are 18 possible rectangles. We are asked to find the grid which contains as close as two million rectangles.

The key here is to count the number of rectangles in a grid efficiently. It seems pretty obvious: Iterate through all the sizes that can fit on the grid and then multiply these by the number of shifts that can be done in the down and in the right direction. This gives a straightforward function to count the rectangles:

def rectangles_in_area(height: int, width: int) -> int:
    result = 0
    for h in range(1, height + 1):
        for w in range(1, width + 1):
            result += (height - h + 1) * (width - w + 1)
    return result

We can make sure that it works using the example as a test case:

def test_rectangles_in_area() -> int:
    assert rectangles_in_area(3, 2) == 18

Since we want to stay close to two million rectangles on the grid, we could try to manipulate height and width of the grid such that we stay close to the goal. Then we could make it an almost one-dimensional optimization problem by following this line on the height-width-plane.

The numbers are so low that one can just try everything on the grid and just take the one with the minimum distance. This is our solution:

def solution() -> int:
    goal = 2_000_000
    results = {}
    for width in range(1, 100):
        for height in range(1, 100):
            num = rectangles_in_area(height, width)
            if num > 2_100_000:
                break
            results[abs(goal - num)] = (height, width, height * width)
    return min(results.items())[1][2]

It takes 509 ms to compute, so that's fast enough.

Factorizing the count of rectangles

There are a couple of trivial things to improve, though. First of all, let us take a look at the number of rectangles on a grid of size $(H, W)$: $$ N(H, W) = \sum_{h=1}^H \sum_{w=1}^W (H - h + 1)(W - w + 1) \,. $$

This sum can be factored into two separate sums: $$ N(H, W) = \sum_{h=1}^H (H - h + 1) \cdot \sum_{w=1}^W(W - w + 1) \,. $$

Both parts are the same, so we can just define $$ N(L) = \sum_{l=1}^L (L - l + 1) $$ and then rewrite it as $$ N(H, W) = N(H) \cdot N(W) \,. $$

Representing this in code we can split this into two functions as well:

def rectangles_along_axis(length: int) -> int:
    return sum((length - l + 1) for l in range(1, length + 1))


def rectangles_in_area(height: int, width: int) -> int:
    return rectangles_along_axis(height) * rectangles_along_axis(width)

This change along brings down the time to 29 ms.

And there is a lot more which can be done. We can cache the results from rectangles_along_axis to make it faster. We then only need 2.2 ms.

@functools.cache
def rectangles_along_axis(length: int) -> int:
    return sum((length - l + 1) for l in range(1, length + 1))

But we don't need to compute both all the time, we can reuse one part of it. Writing the solution like this makes it even faster:

def solution() -> int:
    goal = 2_000_000
    results = {}
    for width in range(1, 100):
        num_along_width = rectangles_along_axis(width)
        for height in range(1, 100):
            num_along_height = rectangles_along_axis(height)
            num = num_along_width * num_along_height
            if num > 2_100_000:
                break
            results[abs(goal - num)] = height * width
    return min(results.items())[1]

And that brings us down to 1.3 ms. That is a ton faster than the 500 ms that we had with the initial solution.

Bisection

We can do much better, still. We don't need to search the whole grid. For each width we can find the nearest two candidates in height via bisection. We precompute the numbers along one axis, and then we can bisect in them.

This is the bisection, one has to be really careful with the indices in this one.

def solution() -> int:
    goal = 2_000_000
    results = {}
    numbers = [rectangles_along_axis(length) for length in range(1, 100)]
    for width, num_width in enumerate(numbers, 1):
        j = bisect.bisect_left(numbers, goal // num_width)
        for height in [j + 1, j + 2]:
            if 0 <= height - 1 < len(numbers):
                num_height = numbers[height - 1]
                num = num_height * num_width
                results[abs(goal - num)] = height * width
    return min(results.items())[1]

That brings down the computation time to 0.064 ms, such that we are already four orders of magnitude faster than before.

This solution is the best that I can come up with, so that has to be sufficient for this problem.