Optimizing an algebraic assignment puzzle

I am going to do a math competition which allows writing programs to solve problems.

The code attempts to solve problem 25 from the Purple Comet Spring 2004 Mathematics Meet:

   W H I T E
+ W A T E R
-----------
P I C N I C


each distinct letter represents a different digit. Find the number represented by the answer PICNIC.

No imports are allowed (tqdm was just a progress bar). What I tried is below. At the rate at which my computer goes, it is not fast enough for the timed competition because it needs to range over 10 to the power of digits there is.

I'm not good at coding, so is there a clear solution to this problem that a beginner could understand? Sorry for the loop tower.

from tqdm import tqdm

def find_solution():
for W in tqdm(range(10)):
for H in tqdm(range(10), desc='H'):
for I in tqdm(range(10), desc='I'):
for T in tqdm(range(10)):
for E in tqdm(range(10)):
for A in (range(10)):
for R in (range(10)):
for P in (range(1, 10)): # P cannot be 0
for C in (range(10)):
for N in (range(10)):
white = W * 10000 + H * 1000 + I * 100 + T * 10 + E
water = W * 10000 + A * 1000 + T * 100 + E * 10 + R
picnic = P * 100000 + I * 10000 + C * 1000 + N * 100 + I * 10 + C

if white + water == picnic:
return {'W': W, 'H': H, 'I': I, 'T': T, 'E': E, 'A': A, 'R': R, 'P': P, 'C': C, 'N': N}

return None

solution = find_solution()

if solution:
print("Solution found:")
print(solution)
else:
print("No solution found.")


I was expecting it to be faster. I originally posted this question on Stack Overflow

• Is this a duplicate of math.stackexchange.com/q/4865353/54983 ? Commented Mar 3 at 3:10
• We can help you write the code better, but for purposes of the competition, you might be better off asking this on Puzzling. While this is solvable via script, there are more efficient ways. Also, you might find that Puzzling would better help you lay out the parameters. For example, P must be 1 and W must be at least 5. It's also worth noting that these puzzles usually don't allow duplicates. So if P is 1, then the other nine can't be. Commented Mar 3 at 15:48

First some context. I understand that you want to practice "competition programming" with no libraries, but:

• you've already violated this yourself;
• that approach deprives you of learning how this problem may be approached in real life with no artificial restrictions;
• by its nature, code submitted to competitions tends to be of lower quality than code written with sufficient care and access to libraries; and
• Code Review considers any insightful observation to be on-topic.

With this in mind, let's discuss the problem and a family of solutions that are efficient and require any of several popular optimisation libraries.

This is an assignment problem. To linearise it, we can use

• $$\0 \le i < m\$$ as the unique-alphabet index,
• $$\0 \le j < n\$$ as the power index (0 being least-significant, meaning that we read right-to-left),
• $$\0 \le k < 10\$$ as the digit index,
• $$\0 \le s_{i,k} \le 1 \quad \forall i,k \$$ as a binary selector matrix for every letter and digit,
• $$\a_j\$$, $$\b_j\$$ and $$\t_j\$$ map vectors of power position to alphabet index.

The addend side (WHITE, WATER) is calculated as $$\sum_j 10^j \sum_k k (s_{{a_j}k} + s_{{b_j}k})$$

and any valid solution must have this equal to the sum side (PICNIC), calculated as $$\sum_j 10^j \sum_k k s_{{t_j}k}$$

You must also constrain the selectors to be unique across both the alphabet and digit dimensions: $$1 = \sum_i s_{i,k} \quad\forall k$$ $$1 = \sum_k s_{i,k} \quad\forall i$$

When implemented in PuLP, this is efficient, beating the accepted answer on your other question in speed by a factor of about 10 (though not as efficient as the runtime claims for the other manual-loop answer).

You can gain a further (minor) improvement in performance by replacing the grand-total form with a per-power form and carry variables $$\c_j\$$: $$10 c_{j+1} + \sum_k k s_{{t_j}k} = c_j + \sum_k k (s_{{a_j}k} + s_{{b_j}k}) \quad \forall j$$

Implicit (grand-total) form

import pulp

a_str = 'WHITE'
b_str = 'WATER'
total_str = 'PICNIC'
letters = sorted({*a_str, *b_str, *total_str})
front_letters = {
seq[0] for seq in (a_str, b_str, total_str)
}

assigns = {
letter: [
# The front digit in each sequence cannot be zero
0 if digit == 0 and letter in front_letters
else (digit == 1) if letter == total_str[0]
else pulp.LpVariable(
name=f'{letter}_{digit}', cat=pulp.LpBinary,
)
for digit in range(10)
]
for letter in letters
}

a, b, total = (
pulp.lpSum(
10**i * pulp.lpDot(
assigns[letter], range(10),
)
for i, letter in enumerate(seq[::-1])
)
for seq in (a_str, b_str, total_str)
)

prob = pulp.LpProblem(name='p25')

# Each letter can only get one assignment
for letter, digits in assigns.items():
digit_total = pulp.lpSum(digits)
if len(digit_total) > 0:
name=f'excl_{letter}', constraint=1 == digit_total,
)
# Each digit can only get one assignment
for digit in range(10):
name=f'excl_{digit}', constraint=1 == pulp.lpSum(
assigns[letter][digit]
for letter in letters
)
)

prob.addConstraint(name='total', constraint=a + b == total)

print(prob)
prob.solve()
assert prob.status == pulp.LpStatusOptimal
a_desc, b_desc, total_desc = (
' '.join(
next(
str(digit)
for digit, asn in enumerate(assigns[letter])
if (value := pulp.value(asn)) > 0.5
)
for letter in seq
)
for seq in (a_str, b_str, total_str)
)

print(f'{a_desc:>12}')
print(f'+ {b_desc:>10}')
print('------------')
print(f'{total_desc:>12}')


outputs

Problem MODEL has 24 rows, 90 columns and 403 elements
Enumerated nodes:               6554
Total iterations:               49049
Time (CPU seconds):             0.39
Time (Wallclock seconds):       0.41

8 3 6 4 2
+  8 5 4 2 7
------------
1 6 9 0 6 9


Explicit (per-power) form

import pulp

a_str = 'WHITE'
b_str = 'WATER'
total_str = 'PICNIC'
letters = sorted({*a_str, *b_str, *total_str})
front_letters = {
seq[0] for seq in (a_str, b_str, total_str)
}

assigns = {
letter: [
# The front digit in each sequence cannot be zero
0 if digit == 0 and letter in front_letters
else int(digit == 1) if letter == total_str[0]
else pulp.LpVariable(
name=f'{letter}_{digit}', cat=pulp.LpBinary,
)
for digit in range(10)
]
for letter in letters
}

prob = pulp.LpProblem(name='p25')

# Each letter can only get one assignment
for letter, digits in assigns.items():
name=f'excl_{letter}', constraint=1 == pulp.lpSum(digits),
)
# Each digit can only get one assignment
for digit in range(10):
name=f'excl_{digit}', constraint=1 == pulp.lpSum(
assigns[letter][digit]
for letter in letters
)
)

# This is the part that varies from the implicit method. We go digit-by-digit.
carries = [
0,
*pulp.LpVariable.matrix(
name='c', cat=pulp.LpContinuous, indices=range(1, len(total_str)),
),
]
a_digits = (*(assigns[l] for l in a_str[::-1]), (0,)*10)
b_digits = (*(assigns[l] for l in b_str[::-1]), (0,)*10)
t_digits = [assigns[l] for l in total_str[::-1]]
for i, (carry, a, b, tot) in enumerate(zip(carries, a_digits, b_digits, t_digits)):
decs = range(10)
total_dec_in = (
carry
+ pulp.lpDot(decs, a)
+ pulp.lpDot(decs, b)
)

if i < len(t_digits) - 1:
next_carry = carries[i + 1]
else:
next_carry = 0

total_dec_out = 10*next_carry + pulp.lpDot(decs, tot)

print(prob)
prob.solve()
assert prob.status == pulp.LpStatusOptimal
a_desc, b_desc, total_desc = (
' '.join(
next(
str(digit)
for digit, asn in enumerate(assigns[letter])
if pulp.value(asn) > 0.5
)
for letter in seq
)
for seq in (a_str, b_str, total_str)
)

print('', ' '.join(f'{pulp.value(c):.0f}' for c in carries[::-1]))
print(f'{a_desc:>12}')
print(f'+ {b_desc:>10}')
print('------------')
print(f'{total_desc:>12}')


outputs

Enumerated nodes:               187
Total iterations:               4206
Time (CPU seconds):             0.36
Time (Wallclock seconds):       0.37

1 0 1 0 0 0
8 3 6 4 2
+  8 5 4 2 7
------------
1 6 9 0 6 9

• Why "claims"? Proof is there, click on the link, then on execute, see for yourself. Commented Mar 4 at 2:59
• Among other reasons, I haven't done an apples-to-apples benchmark on my own hardware Commented Mar 4 at 3:00

Since each letter corresponds to a distinct digit, there are only 10! ≈ 3.6M letter mappings instead of 10^10. These can be generated with itertools.permutations or done manually.

The advantage of doing it manually with loops is that you can quickly exclude large parts of the search space by doing simple checks along the way. For example, once you assign E and R, there can only be one value for C. Then once you assign T, you only have one value for I. Eliminating steps from right-to-left will cut down the search space dramatically and you should be able to brute-force from there. With more simple constraints like P=1, you can probably solve it completely by hand.