Solving systems of polynomials in SymPy

[Maaz Muhammad]

Hi all! I'm a PhD student at the University of Toronto, working on polynomials (optimization, algorithms, and applications). 

As part of my work, I needed to solve systems of polynomials (inequalities and equalities), which SymPy currently lacks. I have therefore been working on such a solver for SymPy. The standard algorithm for this is cylindrical algebraic decomposition, which decomposes R^n into "cells" over which each polynomial is sign-invariant. You can then evaluate the system in each cell, which allows you to either get a satisfying point (or, a satisfying point from each true cell), or the algorithm returns that the system is inconsistent. Here is a demo of my code:

The implementation here returns a single satisfying point if it exists, and otherwise returns an empty list if no satisfying point exists. 

The algorithm relies heavily on the Hong projection operator. This was really the main thing I had to code. The Hong operator relies on taking resultants, reducta, and derivatives -- all of which can already be done in SymPy, so it was just a matter of putting it together (after understanding Hong's paper). 

I previously posted this on my Twitter and was asked by the SymPy Twitter account to post here. I would love to make this part of SymPy. Note that thus far, the only two implementations of CAD are in Mathematica or a software called QEPCAD. 

For now, the algorithm works with sample points from each cell (these sample points are generated during the running of the algorithm). I'm thinking that the poly_solve() function could have a parameter like 'return_type' which can take values either 'one' or 'all' -- the former would return a satisfying point from a single cell, the latter would return a satisfying point from all valid cells. 

Another consideration is that CAD can also be used to construct algebraic formulas that describe all satisfying points. This is called the solution formula. Mathematica and QEPCAD implement this. Unfortunately, this step is quite challenging, and I have not fully understood it yet. Pretty much all the information about this step is given in a PhD thesis from the 90s by Christopher Brown (who actually wrote QEPCAD). I'm studying it right now.

Overall, I have the following questions:

1. What would be the process to integrate this algorithm into SymPy?

2. Should I wait until I figure out the solution formula step to try to integrate into SymPy?

3. Are there any rules about me eventually writing a paper about my implementation in, for example, the INFORMS Journal of Computing (https://pubsonline.informs.org/page/ijoc/editorial-statement#Software%20Tools)? 

Thanks!