Perfect Squares, in depth: the minimum-squares DP, BFS over remainders and the Lagrange-Legendre O(sqrt n) shortcut

Why greedy fails

The tempting approach is greedy: subtract the largest square that fits, repeat. For 12 that takes 9, leaving 3, which needs 1 + 1 + 1, for a total of four squares. The optimum is three (4 + 4 + 4). Greedy fails because taking the biggest square first can leave a remainder that is expensive to finish, and nothing in the greedy rule looks ahead.

This is the same structure as making change with arbitrary coin denominations: the squares are the coins, every coin may be used any number of times, and we want the fewest coins. Greedy is correct only for special coin systems, and the squares are not one of them. When a local choice can be wrong, the standard remedy is to try every choice and remember the results of subproblems, which is dynamic programming.

The recurrence

Let best(i) be the minimum number of squares summing to i. Any optimal sum for i has some last square k*k with k*k at most i, and what remains, i - k*k, must itself be written optimally; if it were not, swapping in a better representation of the remainder would improve the whole. That optimal-substructure argument gives the recurrence:

best(0) = 0
best(i) = 1 + min over k >= 1 with k*k <= i of best(i - k*k)

The base case best(0) = 0 means the empty sum. Every i has at least one valid choice, k = 1, so best(i) is at most i (all ones), which makes i a safe initial value. The subproblems overlap heavily: best(3) is needed by 4, 7, 12 and many others. Computing each once in increasing order of i and storing it turns an exponential search tree into a linear table.

Bottom-up dynamic programming

Bottom-up evaluation fills an array from 0 to n. For each i, the inner loop walks the squares in increasing order and stops when the square exceeds i, so the loop for i runs about sqrt(i) times.

import math

def num_squares_dp(n: int) -> int:
    squares = [k * k for k in range(1, math.isqrt(n) + 1)]
    best = [0] + [n] * n            # best[i] <= i, so n is a safe "infinity"
    for i in range(1, n + 1):
        for s in squares:
            if s > i:
                break
            cand = best[i - s] + 1
            if cand < best[i]:
                best[i] = cand
    return best[n]

Two details matter. math.isqrt computes the integer square root exactly; int(math.sqrt(n)) works for small n but floating-point rounding can be off by one for very large values. And precomputing the list of squares keeps multiplication out of the inner loop. The same code in Java, where the work runs in tight primitive loops:

static int numSquares(int n) {
    int[] best = new int[n + 1];
    for (int i = 1; i <= n; i++) {
        int b = i;                       // all ones
        for (int k = 1; k * k <= i; k++) {
            b = Math.min(b, best[i - k * k] + 1);
        }
        best[i] = b;
    }
    return best[n];
}

Total work is the sum of sqrt(i) for i up to n, which is about (2/3) n^1.5 inner iterations, and memory is n + 1 integers. For n = 10,000 that is roughly 667,000 iterations; for n = 1,000,000 it is roughly 667 million, comfortable in Java or C++ and slow in pure Python.

Worked example: n = 12 by hand

Tracing the table by hand is the fastest way to trust the recurrence. Start with best[0] = 0. Then best[1] = best[0] + 1 = 1, best[2] = best[1] + 1 = 2, best[3] = 3, since only the square 1 fits below 4.

At i = 4 the square 4 fits: best[0] + 1 = 1 beats best[3] + 1 = 4, so best[4] = 1. Continuing: best[5] = min(best[4], best[1]) + 1 = 2, best[6] = min(best[5], best[2]) + 1 = 3, best[7] = min(best[6], best[3]) + 1 = 4. Seven is the first number that needs four squares. Then best[8] = min(best[7], best[4]) + 1 = 2 (4 + 4), best[9] = 1, best[10] = min(best[9], best[6], best[1]) + 1 = 2 (9 + 1), best[11] = min(best[10], best[7], best[2]) + 1 = 3 (9 + 1 + 1).

Finally best[12] = min(best[11], best[8], best[3]) + 1 = min(3, 2, 3) + 1 = 3. To recover the actual squares, store which square achieved the minimum at each cell and walk back: 12 chose 4 and went to 8, 8 chose 4 and went to 4, 4 chose 4 and went to 0. The decomposition is 4 + 4 + 4. Note that greedy's first move, 9, leads to cell 3, which costs 3 more: the table shows exactly why greedy loses.

BFS over remainders

A second view: treat each remainder as a node and draw an edge from r to r - k*k for every square that fits. The answer is the length of the shortest path from n to 0. All edges have weight one, so breadth-first search finds it, and BFS can stop at the first level where it reaches zero instead of filling the whole table.

from collections import deque
import math

def num_squares_bfs(n: int) -> int:
    squares = [k * k for k in range(1, math.isqrt(n) + 1)]
    seen = {n}
    frontier = deque([n])
    depth = 0
    while frontier:
        depth += 1
        for _ in range(len(frontier)):
            r = frontier.popleft()
            for s in squares:
                if s > r:
                    break
                if s == r:
                    return depth          # r itself is a square: reached 0
                nxt = r - s
                if nxt not in seen:
                    seen.add(nxt)
                    frontier.append(nxt)
    raise AssertionError("unreachable: 1 is always a square")

Because the answer is never more than four (next section), BFS explores at most four levels. In practice it touches far fewer than n nodes for most inputs, but the worst case still visits a large share of the remainders, and the visited set costs memory per node. The level loop that processes exactly one depth at a time is the important pattern; forgetting it and incrementing depth per node is the classic bug.

The number theory shortcut

Number theory removes the search altogether. Lagrange's four-square theorem (1770) says every natural number is a sum of four integer squares, allowing zeros, so the answer is always between 1 and 4. Legendre's three-square theorem says exactly which numbers need all four: n is not a sum of three squares if and only if n = 4^a (8b + 7) for non-negative integers a and b. That gives a decision procedure:

import math

def is_square(x: int) -> bool:
    r = math.isqrt(x)
    return r * r == x

def num_squares_math(n: int) -> int:
    if is_square(n):
        return 1
    m = n
    while m % 4 == 0:          # strip factors of 4
        m //= 4
    if m % 8 == 7:
        return 4               # Legendre: 4^a(8b+7)
    for a in range(1, math.isqrt(n) + 1):
        if is_square(n - a * a):
            return 2
    return 3

The order of tests matters. Check for one square first, then the Legendre form, then try every first square for a two-square sum; anything left over must be three. Checking 7: no square, 7 mod 8 = 7, so four (4 + 1 + 1 + 1). Checking 28: strip one factor of 4 to get 7, so four as well (25 + 1 + 1 + 1). Checking 43: not a square, 43 mod 8 = 3, and no a makes 43 - a*a square, so three (25 + 9 + 9). The cost is O(sqrt n) time and O(1) memory, which handles n in the trillions where the table approach cannot even allocate its array.

Failure modes

The bugs that show up in reviews and interview rounds are predictable:

The last row bites the Java version once n approaches the int limit: k * k can wrap negative and still satisfy k * k <= i. The array would not fit in memory at that size anyway, which is itself a signal to use the theorem.

Trade-offs

The DP earns its keep when you answer many queries: build the table once up to the maximum n and every later lookup is O(1). It also generalises to variants the theorems do not cover, such as allowing only some squares, using cubes, or counting the number of representations; those are all the same unbounded-knapsack loop with a different item list or a sum instead of a min. The theorem-based method is the fastest by orders of magnitude but specific to squares, and it relies on results you must state correctly; a reviewer who does not know Legendre's theorem will want a comment citing it.

In an interview, present the DP first, since it shows the reasoning, mention BFS as the shortest-path view, then offer the theorem as the optimisation. In production, use the theorem and keep the DP as a test oracle, exactly as was done for this article.

What to do next

1. Write the DP from memory, then trace best[0..13] by hand and compare with the table above.
2. Add a choice array and return the actual squares, not just the count.
3. Implement the BFS version and confirm it agrees with the DP for every n up to 10,000.
4. Implement the Legendre-based solution, then use the DP as an oracle in a property test.
5. Change the item list to cubes and see which approaches still apply (only the DP and BFS do).
6. Revisit coin change in its two variants and the knapsack family to see the same loop with other items.
7. Read BFS on implicit graphs and dynamic programming fundamentals to generalise both views.