A friend posed the following puzzle (from a real life experience). He had 17 coins which he exchanged for exactly three quarters (laundry day).
What could those coins have been?
Another friend solved this using the following Mathematica™:
Reduce[{
25 q + 10 d + 5 n + p == 75,
q + d + n + p == 17,
q \[Element] Integers, d \[Element] Integers,
n \[Element] Integers, p \[Element] Integers,
q >= 0, d >= 0, n >= 0, p >= 0
}]
... and I'd like to be able to solve it (well, problems of this sort) using Python's SymPy if possible, or some other tool like pyDatalog or kanren or whatever.
However, I don't know quite how to connect the dots. I got (even before seeing the Mathematica solution) that it's a system of two equations:
#!python
import sympy
from sympy import Symbol, Eq
from sympy.solvers import solve
d = Symbol('d', integer=True) # Dimes
n = Symbol('n', integer=True) # Nickels
p = Symbol('p', integer=True) # Pennies
## Assuming he didn't trade quarters for more quarters
f1 = Eq(d + n + p, 17) # Number of coins he had
f2 = Eq(d * 10 + n * 5 + p, 75) # Total value of the coins
But I don't know how to go from there. When I call solve() on this I get:
#!python
solve((f1, f2), (d, n, p))
# ---> {n: -9*p/5 + 19, d: 4*p/5 - 2}
(For reference I know I should be getting a couple of different possible solutions: (d,n,p) => (2,10,5), (6,1,10) ... and there's another that's violating my constraint that he wouldn't exchange quarters for quarters).
I'm sure it's just that I don't understand the tool well enough.
So, how do I solve a system of equations in Sympy, with constraints that the results be non-negative integers? Alternatively how might I do this using a Python logical programming tool or framework?
The way so solve integer equations in sympy is with diophantine from sympy.solvers.diophantine.diophantine. It does not seem to support systems of equations but with a bit of substituting it can still be used:
Frist I reduce the problem to one equation:
term1 = d + n + p - 17
term2 = term1.subs(*solve(f2))
gives
n/2 + 9*p/10 - 19/2
But when I supply this to diophantine it complains it wants integer coefficients so I do
def make_coefs_integers(term):
factor = reduce(sympy.lcm, [x.q for x in term.as_coefficients_dict().values()])
return term*factor
term2 = make_coefs_integers(term2)
giving me
5*n + 9*p - 95
For this diophantine gives me the solution but not quite in the form I need. It returns a set with tuples where the elements of the tuples are solutions for the free sympols in the equation. Since we only have one element in the set I can convert it to a dictionary similar to what solve returns by zipping the solution tuple with the free symbols. I like that form a lot since it can be given to subs.
sol2 = dict(zip(sorted(term2.free_symbols, key=str), tuple(diophantine(make_coefs_integers(term2)))[0]))
which is
{n: 9*t_0 + 190, p: -5*t_0 - 95}
To also get d we put it in to term1 but use t instead of t0 to avoid confusion:
term3 = term1.subs(sol2).subs({t0:t})
again this is:
d + 4*t + 78
Solving once more:
sol3 = dict(zip(sorted(term3.free_symbols, key=str), tuple(diophantine(term3))[0]))
gives
{d: 4*t_0 - 78, t: -t_0}
Which we can combine to a complete solution:
sol4 = {k:v.subs({t0:-t0}) for k,v in sol2.items()} | {d:sol3[d]}
which is
{n: 190 - 9*t_0, p: 5*t_0 - 95, d: 4*t_0 - 78}
{k:v.subs({tuple(v.free_symbols)[0]:20}) for k,v in sol4.items()}
being
{n: 10, p: 5, d: 2}
and
{k:v.subs({tuple(v.free_symbols)[0]:21}) for k,v in sol4.items()}
being
{n: 1, p: 10, d: 6}