You can try the solvers I developed using the methods described below at pipip.
Recently, I was quite amused to see this GitHub project: sudoku-in-python-packaging. It's a Sudoku solver which uses Python packaging -- not Python, just the packaging dependency resolver -- to solve Sudoku puzzles. Basically you can formulate a Sudoku puzzle as a Python package dependency resolution problem and then call pip (or uv) to solve it.
I found this idea so amusing that I wanted to extend it. Package dependency resolution is NP-complete, meaning that any Boolean satisfiability (SAT) or (0/1) integer programming (IP) problem can be formulated as a package dependency resolution problem. This means that we can build a SAT and IP solver that purely relies on pip. No other libraries are needed except for the Python standard library. This of course means that you can solve a massive number of problems with pip, including classical problems like knapsack, subset sum, graph coloring, and so on. Below I will describe how this works.
In complexity theory, specifically when in the domain of NP-complete problems, all roads lead to SAT. So it was natural to start with SAT. To recap, a SAT problem in conjunctive normal form (CNF) is a conjunction (i.e., AND) of clauses, where each clause is a disjunction (i.e., OR) of literals. It is given as: \[ \bigwedge_{i=1}^m \bigvee_{j=1}^n \ell_{ij} \] where each \(\ell_{ij}\) is a literal, which is either a variable or its negation. An example of a SAT instance is: \[ (x_1 \vee \neg x_2 \vee x_4) \wedge (\neg x_1 \vee x_3) \wedge (x_2 \vee x_3 \vee \neg x_4) \wedge (\neg x_2 \vee \neg x_4) \] And the goal is to find an assignment of these Boolean variables that makes the expression true. If such an assignment exists, we say that the SAT instance is satisfiable (SAT), otherwise it is unsatisfiable (UNSAT).
The first thing to note is what sort of capabilities Python package dependency specification has. First for a given package, only a single version can be installed at a time. This allows us to encode Boolean variables as package versions, which ensures that each variable is either true or false. Second, it allows for the following comparison operators: \( \geq, >, =, <, \leq \). If multiple versions are given separated by a comma, the comma acts as a logical AND. Importantly, it does not allow for the logical OR operator, which means we'll have to work around that.
Given this, here is how we will encode a SAT instance as a package dependency resolution problem. We will make "dummy" packages for every variable and clause. For the packages that are variables, we will give them two possible versions: 1.0 for false and 2.0 for true. For the packages that are clauses, we will give them multiple versions: 1.0, 2.0, 3.0, and so on, until the number of literals in that clause. We specify the following dependencies for the clause packages: version 1.0 depends on the first literal being true, version 2.0 depends on the second literal being true, and so on. That is all!
Why does this work? A proof of the reduction is basically this: if the SAT instance is satisfiable, then assign the variable packages the appropriate versions. As it is satisfiable, at least one of the literals in each clause is true, meaning each clause package can be installed. Conversely, if the package installation is possible, then assign the variables to be true or false depending on the version installed. And, since each clause package is installed, then at least one of the literals in each clause must be true, so the SAT instance is satisfiable.
Here's an example. Consider the following SAT instance: \[ (x_1 \vee x_2 \vee \neg x_3) \wedge (x_1 \vee x_3 \vee x_4) \wedge (\neg x_2 \vee x_3 \vee x_4) \] From above procedure, we introduce the following packages:
x1, x2, x3, x4, each with versions 1.0 (false) and 2.0 (true)c1, c2, c3, each with versions 1.0, 2.0, and 3.0
c1==1.0 requires x1==2.0 # x_1 is true
c1==2.0 requires x2==2.0 # x_2 is true
c1==3.0 requires x3==1.0 # x_3 is false
c2==1.0 requires x1==2.0 # x_1 is true
c2==2.0 requires x3==2.0 # x_3 is true
c2==3.0 requires x4==2.0 # x_4 is true
c3==1.0 requires x2==1.0 # x_2 is false
c3==2.0 requires x3==2.0 # x_3 is true
c3==3.0 requires x4==2.0 # x_4 is true
x1==2.0x2==2.0x3==2.0x4==1.0c1==1.0c2==1.0c3==2.0c1==2.0, because we have \( x_3 \) being true. But all that matters is that the packages are installable.
I implemented this in pipip ("pip" + "ip", and I just thought it sounded funny). The file pipip/satsolver.py contains the SAT solver. It basically implements the above procedure. The extra stuff is in how it makes the dummy packages and dependencies. I definitely needed help from Claude to get this working. To make the packages, it makes wheel files, and then generates a requirements.in file to specify the dependencies. Then it calls pip-compile to compile the dependencies into a requirements.txt file. If this works, then it means the SAT instance is satisfiable. Otherwise, if it can't generate the file, then the SAT instance is unsatisfiable.
To use it you can either import it into your own project, or you can run it from the command line by passing a DIMACS file with a .cnf extension as follows: python satsolver.py sat01.cnf. It will output the solution if satisfiable, otherwise it will output unsatisfiable.
As I said before, all roads lead to SAT. And as 0/1 integer linear programming (IP) is also NP-complete, we will try to transform an IP into a SAT instance. Note that by IP here I mean the decision version, where the goal is to determine if there exists an integer solution that satisfies a given set of linear inequalities.
An IP instance is of the form: \[ Ax \leq b \] where \(A \in \mathbb{Z}^{m \times n}\) is a matrix of coefficients, \(x \in \{0,1\}^n\) is a column vector of variables, and \(b \in \mathbb{Z}^m\) is a column vector of constants. The variables are constrained to be 0 or 1 (we could relax this by introducing binary encodings, but let's keep it simple). We also assume that the entries of \(A\) and \(b\) are integers (this is not a restrictive assumption, because fixed-precision numbers are rational numbers, and we can convert them to integers by scaling). We seek to find a vector \(x\) that satisfies the system, i.e., to see if the system is feasible.
We will encode this into a SAT instance through a series of transformations and encodings. To build it up, let's first consider a simpler version, where we have a single linear inequality, with all positive coefficients (including the constant term).
In this special case, we have an inequality of the form \( \sum_{i=1}^n w_i x_i \leq b \), where all \(w_i\) and \(b\) are positive. Of course this is always feasible by setting all \(x_i\) to 0. The important part we will develop here is how to encode its feasibility as a SAT instance. With a slight abuse of notation, we will also use \(x_i\) to refer to the corresponding variable in the SAT instance.
The key idea to build our SAT instance is to introduce a series of new binary variables \( s_{i,j} \), where the variable \( s_{i,j} \) will be true if and only if some (possibly empty) subset of the first \(i\) variables sums exactly to \(j\). How many such variables are there? Well, we can take the sum of the weights \( w_i \) to get \( W = \sum_{i=1}^n w_i \). Then there are \( W+1 \) possible values for the sum, ranging from 0 to \( W \). So, we introduce \( (n+1) \times (W+1) \) variables, meaning the indices for \(s_{i,j}\) range from \(i = 0 ... n\) and \(j = 0 ... W\).
To encode the relationship between the \(s_{i,j}\) and the \(x_i\), we need to write a series of implications. It will be similar to a dynamic programming approach, as we have a recursive relationship between the \(s_{i,j}\) variables. The idea is, for each \(i, j\), to compare the weight of the \(i\)-th variable to the sum \(j\). We have two cases:
We now need to write the biimplications as CNF clauses. By applying a series of properties of logical operations (especially the fact that \(P \implies Q \iff \neg P \lor Q\)), we can convert the above implications into CNF clauses. I'll skip the details, but the reader can easily verify that the following clauses encode the above implications:
Lastly, we need to actually enforce the bound \( \leq b \). We do this by ensuring that \(s_{n,j}\) is false for all \(b < j \leq W\), i.e., the sum of the \(n\) variables should not exceed \(b\). This is done by conjuncting the following (single-literal) clauses \(\neg s_{n,j}\) for all \(b < j \leq W\).
Here's an example. Consider the following IP instance: \[ x_1 + 2x_2 \leq 2 \] Because the maximum sum \(W = 1 + 2 = 3\), we introduce \( s_{i,j} \) variables for \(i = 0, 1, 2\) and \(j = 0, 1, 2, 3\). These can be thought of as a dynamic programming table, with the biimplications we derived above. Below, \(i \) lies along the top and \(j\) lies along the left, and the entry is the equivalent for that \(s_{i,j}\). We start by initializing the base case:
| i=0 | i=1 | i=2 | |
|---|---|---|---|
| j=0 | 1 | ||
| j=1 | 0 | ||
| j=2 | 0 | ||
| j=3 | 0 |
Next, let's look at \(s_{1,0}\). Because \(w_1 = 1 > 0\), we have \(s_{1,0} \iff s_{0,0} \land \neg x_1\). The situation is different for \(s_{1,1}\), because \(w_1 = 1 \leq 1\), so we have \(s_{1,1} \iff (s_{0,1} \land \neg x_1) \lor (s_{0,0} \land x_1)\). We can keep going like this to fill up the table. The final table is as follows:
| i=0 | i=1 | i=2 | |
|---|---|---|---|
| j=0 | 1 | \(s_{0,0} \land \neg x_1\) | \(s_{1,0} \land \neg x_2\) |
| j=1 | 0 | \((s_{0,1} \land \neg x_1) \lor (s_{0,0} \land x_1)\) | \(s_{1,1} \land \neg x_2\) |
| j=2 | 0 | \((s_{0,2} \land \neg x_1) \lor (s_{0,1} \land x_1)\) | \((s_{1,2} \land \neg x_2) \lor (s_{1,0} \land x_2)\) |
| j=3 | 0 | \((s_{0,3} \land \neg x_1) \lor (s_{0,2} \land x_1)\) | \((s_{1,3} \land \neg x_2) \lor (s_{1,1} \land x_2)\) |
(I checked the table carefully, so I hope I didn't make any typos above.) Lastly, we encode the bounding constraint \( \leq 2 \) by adding the clause \(\neg s_{2,3}\), i.e., the sum of the variables should not exceed 2.
Note of course that we could propagate some of the known values and simplify the clauses. For example, because \(s_{0,0} = 1\), we can immediately infer that \(s_{1,0} = \neg x_1\), which makes sense because if \(x_1\) was true, then the sum would be 1, so the only way to make the sum 0 with the first \(1\) variables is if \(x_1\) is false. We won't do this here, and instead leave such simplification to the SAT solver (this is called "unit propagation" in SAT).
The biimplications above are converted into CNF clauses as above, and of course we also need to enforce the bound \( \leq 2 \). I won't write out all of the clauses here. Even for such a simple example, there are many variables and clauses. To be exact, we have \( n + (n+1) \times (W+1) \) variables, or in big-O notation, \(O(nW)\). To count the clauses, we have:
Now that we can encode a single inequality with all positive coefficients, we will progressively build up to the general case.
If some of the weights in \(\sum_{i=1}^n w_i x_i \leq b\) are negative, but \(b\) is still positive, we can't exactly do the same sum tracking trick, because including a variable in the sum can now detract from the sum. But, we can introduce for each variable \(x_i\) a new variable \(z_i = \neg x_i\), or, in arithmetic, \(z_i = 1 - x_i\). Then, for any \(i\) with negative weight \(w_i\), we can see that the term \(w_i x_i = w_i (1-z_i) = w_i - w_i z_i\). We subtract \(w_i\) from both sides, so we have \(-w_i z_i \leq b - w_i\). Because \(w_i\) is negative, \(b - w_i\) is positive, and \(-w_i z_i\) is a term with a positive weight, so we can apply the same sum tracking trick as before.
As an example, consider \(2x_1 - 3x_2 \leq 1\). We introduce \(z_1 = 1 - x_1\) and \(z_2 = 1 - x_2\). Then the inequality becomes \(2x_1 - 3(1-z_2) \leq 1\), which simplifies to \( 2x_1 - 3 + 3z_2 \leq 1 \), which is equivalent to \(2x_1 + 3z_2 \leq 4\). And now we can apply the same sum tracking trick as before.
If the bound \(b\) is negative, we can multiply the inequality by \(-1\) to get \(\sum_{i=1}^n -w_i x_i \geq -b\). The signs of the weights don't matter, as we have just addressed this above and can convert them into all positive weights. Now we can still do the same sum tracking trick, but the final bounding constraint is different. Instead of fixing \(s_{n,j}\) to false for all \(b < j \leq W\), we now require at least one of the \(s_{n,j}\) to be true for all \(b \leq j \leq W\). This can be done by adding the single clause \(s_{n,b} \lor s_{n,b+1} \lor \cdots \lor s_{n,W}\).
We now have all the tools to build a SAT instance for \(Ax \leq b\). First, we have the variables \(x_i\) and the new variables \(z_i\). We enforce that they are each others negation by adding the clauses \(x_i \lor z_i\) and \(\neg x_i \lor \neg z_i\). Then, for each inequality we do the sum tracking trick as above and introduce the bounding constraint. In total, the number of variables is \(O(m n \overline{W})\), where \(\overline{W}\) is the maximum sum of the absolute values of the weights of each row of \(A\) -- recall that \(m\) is the number of inequalities or the number of rows of \(A\). This is because we have \(2n\) variables (the \(x_i\) and \(z_i\)) plus \(O(nW)\) for each inequality, for a total of \(O(m n \overline{W})\). Similarly, the number of clauses for each inequality is \(O(nW)\), plus the \(2n\) clauses that define the relationship between \(x_i\) and \(z_i\), for a total of \(O(m n \overline{W})\).
I implemented this in pipip The file pipip/ipsolver.py contains the IP solver. It implements the above encoding, and then calls the SATsolver in pipip/satsolver.py to solve the SAT instance. It decodes the SAT solution to get the IP solution. You can either use it by importing into your own project, and pass it \(A\) as a list of lists, and \(b\) as a list, or you can run it from the command line by passing a plain-text file that contains \(A\) and \(b\) in the following format:
a_11 a_12 ... a_1n
a_21 a_22 ... a_2n
...
a_m1 a_m2 ... a_mn
b_1 b_2 ... b_m
1 2 3
2 3 4
ipsolver.py also contains functions for converting it into a CNF file, so you can also inspect the encodings themselves if you are interested.
There are several examples contained in the examples/ folder, which contains two folders: sat and ip.
The file examples/sat/001.cnf is an easy instance:
\[
(x_1 \lor x_2 \lor \neg x_3) \land (x_1 \lor x_3 \lor x_4) \land (\neg x_2 \lor x_3 \lor x_4)
\]
Running python satsolver.py examples/sat/001.cnf will output that it is satisfiable, and a solution is \(x_1 = 1, x_2 = 1, x_3 = 0, x_4 = 1\), which is easily verified to be correct.
The file examples/sat/002.cnf is an easy unsatisfiable instance:
\[
x_1 \land \neg x_1
\]
Running python satsolver.py examples/sat/002.cnf will output that it is unsatisfiable.
I took two benchmarks from the SATLIB benchmarks. I took a satisfiable instance and an unsatisfiable instance, in examples/sat/004.cnf and examples/sat/005.cnf, respectively, both with 50 variables and 218 clauses. They are 3-SAT instances, meaning each clause has exactly 3 literals.
Here I started running into some issues. The satisfiable instance was returning unsatisfiable, after around 5 minutes, and I couldn't figure out why. After doing some digging around, I realized that I was running into the time limits of pip-compile. There are some ways to influence this, e.g., by changing the max-rounds option, but I didn't have luck with this. I then thought about using uv as the package dependency manager instead, which is known to be much faster. I believe that uv actually uses a SAT solver in the background, so you might argue this whole exercise has been rather circular. Anyway, I added the argument --uv to satsolver.py, so you can call it like python satsolver.py examples/sat/004.cnf --uv. Now it correctly returns satisfiable, and it took 0.572 seconds, and the unsatisfiable instance took 1.514 seconds, on my Macbook Air M3.
The file examples/ip/001.txt is an easy IP instance containing \(2x_1 + 3x_2 \leq 1\). This is easily feasible by setting both to \(0\) and indeed this is the solution found by the solver. It is important to note here that with pip-compile, the solver returned infeasible, but again this is likely due to the computational limits of pip-compile. With the --uv flag, it correctly returns feasible. The file examples/ip/002.txt is an infeasible IP instance with \(2x_1 + 3x_2 \leq -1\), and the solver correctly returns infeasible. A slightly more interesting example is the file examples/ip/003.txt, which contains the following IP instance:
\[
x_1 + x_2 - x_3 \leq 2
\]
\[
-x_1 + x_2 + x_3 \leq 1
\]
The solver correctly returns feasible, and a solution is \(x_1 = 1, x_2 = 1, x_3 = 1\).
In examples/ip/004.txt, there is a subset sum problem with the numbers \( \{3, 9, 20, 25, 30, 45\} \) and the target sum is 76. By inspection this is clearly infeasible. To be explicit, the IP instance is:
\[
3x_1 + 9x_2 + 20x_3 + 25x_4 + 30x_5 + 45x_6 \leq 76
\]
\[
-3x_1 - 9x_2 - 20x_3 - 25x_4 - 30x_5 - 45x_6 \leq -76
\]
The solver correctly returns infeasible in 2 min 2.12 sec.
The file examples/ip/005.txt contains a set-covering problem with the following sets:
\[
S_1 = \{1, 2, 3\},
S_2 = \{3, 4\},
S_3 = \{2, 5\},
S_4 = \{1, 4, 5\}
\]
And we ask for a set cover of size at most 3. The solver correctly returns feasible in 0.609 seconds, with the solution \(S_1, S_2, S_3\).
I hope this gives a sense of the power of NP-completeness -- any problem that is NP-complete can be transformed into any other NP-complete problem. If you're stranded on an island with nothing but a standard Python install, you can solve a wide variety of problems, no Gurobi required! A fun exercise (depending on who you ask?) might be to try solving some Leetcode problems using this method as many of them are basically combinatorial problems that can be encoded as SAT or IP.