Sudoku Puzzle with boxes containing square numbers

This is where I would use a SMT solver. They are a lot more powerful than people give credit for. If the best algorithm you can think of is essentially bruteforce, try a solver instead. Simply listing your constraints and running it gives your unique answer in a couple seconds:

278195436
695743128
134628975
549812763
386457291
721369854
913286547
862574319
457931682

The code used (and reference image for coordinates):

import z3

letters = "ABCDEFGHI"
numbers = "123456789"
boxes = """
A1 A2 A3
B1 B2 C2 C3
C1 D1 D2
E1 E2 F2
F1 G1
H1 I1
G2 H2 G3 H3 H4
I2 I3 I4
B3 B4 C4
D3 E3 F3
A4 A5 B5
C5 B6 C6
G5 H5 I5 I6
A6 A7
B7 C7
D7 D8 D9
E7 E8 F7 F8
G7 H7
I7 I8
A8 B8 C8
G8 H8
A9 B9 C9
E9 F9
G9 H9 I9
"""
positions = [letter + number
             for letter in letters
             for number in numbers]
S = {pos: z3.Int(pos) for pos in positions}

solver = z3.Solver()

# Every symbol must be a number from 1-9.
for symbol in S.values():
    solver.add(z3.Or([symbol == i for i in range(1, 10)]))

# Every row value must be unique.
for row in numbers:
    solver.add(z3.Distinct([S[col + row] for col in letters]))

# Every column value must be unique.
for col in letters:
    solver.add(z3.Distinct([S[col + row] for row in numbers]))

# Every block must contain every value.
for i in range(3):
    for j in range(3):
        solver.add(z3.Distinct([S[letters[m + i * 3] + numbers[n + j * 3]]
                                for m in range(3)
                                for n in range(3)]))

# Colored boxes.
for box in boxes.split("\n"):
    box = box.strip()
    if not box: continue
    boxsum = z3.Sum([S[pos] for pos in box.split()])
    solver.add(z3.Or([boxsum == 1, boxsum == 4, boxsum == 9,
                      boxsum == 16, boxsum == 25, boxsum == 36]))

# Print solutions.
while solver.check() == z3.sat:
    model = solver.model()
    for row in numbers:
        print("".join(model.evaluate(S[col+row]).as_string()
                    for col in letters))
    print()

    # Prevent next solution from being equivalent.
    solver.add(z3.Or([S[col+row] != model.evaluate(S[col+row])
                      for col in letters
                      for row in numbers]))