Generic Cylindrical Algebraic Decomposition

Given a list of polynomials in $x_1,…,x_n$, GCAD separates $ℝ^n$ into cells of constant sign of all polynomials, each given as

$$\bigwedge_{1≤i≤n} l_i(x_1,…,x_{i-1}) < x_i < h_i(x_1,…,x_{i-1}),$$

where the boundaries $l_i$ and $h_i$ are either ±∞, or real roots of polynomials.

Conventionally, the input polynomials are given as inequalities, and only the cells that satisfy them are returned.

How to use

To decompose $2-x^2-y^2-z^2>0 ∧ 1-x^2-(y-1)^2-z^2>0 ∧ 1-z^2>0$ into cells:

import sympy as sp
import gcad

x, y, z = sp.symbols("x y z")
p1 = 2 - x**2 - y**2 - z**2
p2 = 1 - x**2 - (y-1)**2 - z**2
p3 = 1 - z**2
cells = gcad.merge(gcad.GCAD([p1 > 0, p2 > 0, p3 > 0], [x, y, z]))
for i, cell in enumerate(cells):
    print(f"Cell {i}:")
    for axis in cell:
        print(f"  {axis.cell_lo} < {axis.var} < {axis.cell_hi}")
print()
for i, cell in enumerate(cells):
    print(f"Cell {i} example point:")
    for axis in cell:
        print(f"  {axis.var} = {axis.point}")

The expected output is:

Cell 0:
  Root[#+1, 0] < x < Root[#-1, 0]
  Root[#^2-2*#+x^2, 0] < y < Root[#-1, 0]
  Root[#^2+x^2+y^2-2*y, 0] < z < Root[#^2+x^2+y^2-2*y, 1]
Cell 1:
  Root[#+1, 0] < x < Root[#-1, 0]
  Root[#-1, 0] < y < Root[#^2+x^2-2, 1]
  Root[#^2+x^2+y^2-2, 0] < z < Root[#^2+x^2+y^2-2, 1]

Cell 0 example point:
  x = -7/8
  y = 2/3
  z = 0
Cell 1 example point:
  x = -7/8
  y = 17/16
  z = 0

How to install

The latest release is available from PyPI.org/project/gcad, and can be installed via pip with

pip install gcad

If you prefer to experiment with the source code, then clone it, make sure you have a C++ compiler toolchain (GCC or Clang), install the dependencies listed in pyproject.toml, and build with make. Run tests with make test.

Note that the build system will build FLINT, MPFR, and GMP automatically (via hepware).

If the desired compiler is located in a non-standard location, you can set CC, CFLAGS, CXX, CXXFLAGS, and LDFLAGS environment variables before attempting installation to convey this information to the build system.

References

This implementation is based on the work of A. Strzeboński (S00). The root finding used here is based on the works of A. Akritas, A. Strzeboński, and others (ASV08, AAS08, ASV06, AS05). The variable ordering heuristics are from DSS04, dRE22, and PdRAEC23. The handling of rationals and polynomials is done by FLINT, MPFR, and GMP.

License

This library is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

API

Functions

GCAD (relations: Expr | list[Expr], varlist: list[Symbol]) → list[list[AxisBound]]

Generic Cylindrical Algebraic Decomposition (Algorithm 3.1 of S00). Take one or multiple multivariate polynomial inequalities, and decompose the region their conjunction defines into a set of cells. (The list of boundaries from the original algorithm is not computed here).

GPROJ (positives: list[Poly], varlist: list[Symbol]) → list[list[Poly]]

Generic projection (Algorithm 3.4 of S00). Return the list of $pr_i$, as polynomials with integer coefficients.

RSFC (positives: list[Poly], pr: list[list[Poly]], varlist: list[Symbol], max_cells: int | None) → list[list]

Recursive Solution Formula Construction (Algorithm 3.5 of S00).

find_instance (relations: Expr | list[Expr], varlist: list[Symbol]) → list[Rational] | None

Find a rational point that lies inside a region defined by the conjunction of the given polynomial relations (using GCAD), and return its coordinates, or None, if no such point exists.

greedy_mods_order (relations: Expr | list[Expr], var_groups: list[list[Symbol]]) → list[Symbol]

Variable order that greedily minimizes the "multiplication of degree sum" metric, as advocated in dRE22.

greedy_sotd_order (relations: Expr | list[Expr], var_groups: list[list[Symbol]]) → list[Symbol]

Variable order that greedily minimizes the "sum of total degrees" metric, as advocated in DSS04.

greedy_t1_order (relations: Expr | list[Expr], var_groups: list[list[Symbol]]) → list[Symbol]

Variable order that greedily minimizes the "t1" metric, as advocated in PdRAEC23.

merge (cells: list[list[AxisBound]]) → list[list[AxisBound]]

Merge adjacent cells (Remark 3.7 of S00) using an aggressive merging strategy, merging cells even if the inequalities may not hold at the boundaries between the cells.

relations_to_positives (relations: Expr | list[Expr], varlist: list[Symbol]) → list[Poly]

Turn one or more $a>b$ or $a<b$ relations into a list of positive expressions (i.e. turn $a>b$ into $a-b$, and $a<b$ into $b-a$).

Classes

AxisBound. Cell boundaries along one axis.

• var: Symbol. Variable name corresponding to this axis.
• point: Rational. A sample var value inside the cell.
• cell_lo: PolyRoot | None. Lower bounding value. None means negative infinity.
• cell_hi: PolyRoot | None. Higher bounding value. None means positive infinity.

PolyRoot. Real root of a polynomial.

• poly: Poly. Polynomial defining the root.
• var_idx: int. Index of the variable in which to find the roots.
• idx: int. Index of the real root, if roots were ordered by value.
• value_lo: Rational. Lower bound on the root's value at the sample point.
• value_hi: Rational. Upper bound on the root's value at the sample point.