Example Worksheet - Maple Help

All Products Maple MapleSim

Home : Support : Online Help : Mathematics : Factorization and Solving Equations : RegularChains : Example Worksheet

Studying and solving polynomial systems with the RegularChains library

Parisa Alvandi, Changbo Chen, James Davenport, Juergen Gerhard, Mahsa Kazemi, Francois Lemaire, Liyun Li, Xin Li, John May, Marc Moreno Maza, Wei Pan, Bican Xia, Rong Xiao, and Yuzhen Xie

1. Introduction

The RegularChains library offers a variety of commands for solving polynomial systems symbolically and studying their solutions. The input systems may contain polynomial equations $p = 0$ , polynomial inequations $p \neq 0$ , and polynomial inequalities $p > 0$ or $p \geq 0$ . All of those may be nonlinear. For a given system, the coefficients of the input polynomials may be rational numbers, integers modulo a prime, or polynomials depending on parameters. The solution sets computed by the RegularChains library are described by lists of components. Geometrically, such a component can be a point (or a set of points), a curve (or a set of curves), a surface (or a set of surfaces), etc. Computationally, components are represented by a special kind of polynomial system with a triangular shape and other algebraic properties. Depending if the input system consists of equations only, has inequations but no inequalities, or possesses inequalities, the system representing a component is called a regular chain, a regular system or a regular semi-algebraic system, respectively. The word "regular" refers to the interesting algebraic properties that these systems have. The RegularChains library provides types for these different kinds of polynomial systems; this important feature will be illustrated hereafter.

Another design feature of the RegularChains library is the organization of its 146 commands into 8 modules. The top-level of the library gathers the most commonly used commands, whereas the six submodules are dedicated to special topics. The submodules ChainTools, ConstructibleSetTools, and SemiAlgebraicSetTools deal respectively with specific operations on regular chains, regular systems, and regular semi-algebraic systems. The submodule MatrixTools allows the user to handle linear systems over domains with zero-divisors. This is a fundamental tool in the theory of regular chains with many potential applications. The AlgebraicGeometryTools submodule provides facilities for studying algebraic curves, surfaces and algebraic sets of higher dimension.

The submodule ParametricSystemTools is devoted to solving systems with parameters, including real root classification and complex root classification of such systems, as demonstrated below. Finally, the submodule FastArithmeticTools provides highly optimized implementation of some fundamental operations on regular chains, similar to some commands of the top-level module and ChainTools. However, the commands of this submodule can only be used under some assumptions and their usage requires care and advanced knowledge.

>	restart; with(RegularChains); with(ChainTools); with(MatrixTools); with(ConstructibleSetTools); with(ParametricSystemTools); with(SemiAlgebraicSetTools); with(FastArithmeticTools); with(AlgebraicGeometryTools);

$[AlgebraicGeometryTools, ChainTools, ConstructibleSetTools, Display, DisplayPolynomialRing, Equations, ExtendedRegularGcd, FastArithmeticTools, Inequations, Info, Initial, Intersect, Inverse, IsRegular, LazyRealTriangularize, MainDegree, MainVariable, MatrixCombine, MatrixTools, NormalForm, ParametricSystemTools, PolynomialRing, Rank, RealTriangularize, RegularGcd, RegularizeInitial, SamplePoints, SemiAlgebraicSetTools, Separant, SparsePseudoRemainder, SuggestVariableOrder, Tail, Triangularize]$

$[Chain, ChangeOfCoordinates, ChangeOfOrder, Construct, Cut, DahanSchostTransform, Dimension, Empty, EqualSaturatedIdeals, EquiprojectableDecomposition, Extend, ExtendedNormalizedGcd, IsAlgebraic, IsEmptyChain, IsInRadical, IsInSaturate, IsIncluded, IsPrimitive, IsStronglyNormalized, IsZeroDimensional, IteratedResultant, LastSubresultant, Lift, ListConstruct, NormalizeRegularChain, NumberOfSolutions, Polynomial, Regularize, RemoveRedundantComponents, SeparateSolutions, Squarefree, SquarefreeFactorization, SubresultantChain, SubresultantOfIndex, Under, Upper]$

$[IsZeroMatrix, JacobianMatrix, LowerEchelonForm, MatrixInverse, MatrixMultiply, MatrixOverChain]$

$[Complement, ConstructibleSet, CylindricalDecompose, Difference, EmptyConstructibleSet, GeneralConstruct, Intersection, IsContained, IsEmpty, MakePairwiseDisjoint, PolynomialMapImage, PolynomialMapPreimage, Projection, QuasiComponent, RationalMapImage, RationalMapPreimage, RefiningPartition, RegularSystem, RegularSystemDifference, RepresentingChain, RepresentingInequations, RepresentingRegularSystems, SeparateZeros, Union]$

$[BelongsTo, BorderPolynomial, CoefficientsInParameters, ComplexRootClassification, ComprehensiveTriangularize, DefiningSet, DiscriminantSequence, DiscriminantSet, PreComprehensiveTriangularize, RealComprehensiveTriangularize, RealRootClassification, Specialize]$

$[BoxValues, Complement, CylindricalAlgebraicDecompose, Difference, DisplayParametricBox, DisplayQuantifierFreeFormula, EmptySemiAlgebraicSet, Intersection, IsContained, IsEmpty, IsParametricBox, LinearSolve, PartialCylindricalAlgebraicDecomposition, PositiveInequalities, Projection, QuantifierElimination, RealRootCounting, RealRootIsolate, RefineBox, RefineListBox, RemoveRedundantComponents, RepresentingBox, RepresentingChain, RepresentingQuantifierFreeFormula, RepresentingRootIndex, SignAtBox, VariableOrdering]$

$[BivariateModularTriangularize, IteratedResultantDim0, IteratedResultantDim1, NormalFormDim0, NormalizePolynomialDim0, NormalizeRegularChainDim0, RandomRegularChainDim0, RandomRegularChainDim1, ReduceCoefficientsDim0, RegularGcdBySpecializationCube, RegularizeDim0, ResultantBySpecializationCube, SubresultantChainSpecializationCube]$

$[Cylindrify, IntersectionMultiplicity, IsTransverse, LimitPoints, RationalFunctionLimit, RegularChainBranches, TangentCone, TangentPlane, TriangularizeWithMultiplicity]$

(1.1)

2. Polynomial systems and regular chains

This section introduces the user to the top-level module of the RegularChains library. The concept of a regular chain is based on a recursive and univariate vision of a polynomial. Several commands manipulating polynomials in this manner are illustrated below. The most commonly used command of the library is called Triangularize and is illustrated by several examples. It takes a polynomial system as input and returns a description of its solution set. If the input system consists of equations only, then this description is a list of regular chains. The case where the input system admits inequations can be handled by Triangularize or by the command GeneralConstruct of the ConstructibleSetTools submodule. In case of inequalities, and more generally for computing the real solutions of a polynomial system, the commands RealTriangularize, LazyRealTriangularize and SamplePoints of the top-level module can be used. Different commands are also available in the modules ParametricSystemTools and SemiAlgebraicSetTools for the same purpose of real solving but for particular types of systems: they will be presented later in this document.

2.1 A univariate vision of polynomials

Define the polynomial ring $&rationals; [x, y, z]$ with the variable ordering $x > y > z$ . The field $&rationals;$ of rational numbers is the default coefficient ring.

>	R := PolynomialRing([x, y, z]);

$R ≔ polynomial_ring$

(2.1.1)

The internal representation of this polynomial ring can be accessed as follows.

>	Display(R);

$\{\begin{array}{l} Variables & &colon; & [x, y, z] \\ Parameters & &colon; & \emptyset \\ Characteristic & &colon; & 0 \end{array}$

(2.1.2)

Consider the following polynomial.

>	p := (y+1)x^3+(z+4)x+3;

$p ≔ (y + 1) x^{3} + (z + 4) x + 3$

(2.1.3)

The main variable of $p$ as a polynomial of $R$ is given by:

$(y + 1) x^{3} + (z + 4) x + 3$

$polynomial_ring$

(2.1.4)

>	MainVariable(p, R);

$x$

(2.1.5)

The initial (or leading coefficient) of $p$ :

$(y + 1) x^{3} + (z + 4) x + 3$

(2.1.6)

>	Initial(p, R);

$y + 1$

(2.1.7)

The degree of $p$ in the main variable:

$(y + 1) x^{3} + (z + 4) x + 3$

(2.1.8)

>	MainDegree(p, R);

$3$

(2.1.9)

The rank of $p$ , which is the main variable of $p$ raised to the main degree of $p$ :

$(y + 1) x^{3} + (z + 4) x + 3$

(2.1.10)

>	Rank(p, R);

$x^{3}$

(2.1.11)

The tail of $p$ .

>	Tail(p, R);

$x z + 4 x + 3$

(2.1.12)

The separant of $p$ , which is the derivative of $p$ with regards to its main variable:

>	Separant(p, R);

$3 x^{2} y + 3 x^{2} + z + 4$

(2.1.13)

Change the ordering to $z > y > x$ :

>	R := PolynomialRing([z, y, x]); Display(R); p := expand((y+1)x^3+(z+4)x+3); MainVariable(p, R); Initial(p, R); Rank(p, R); Tail(p, R); Separant(p, R);

$R ≔ polynomial_ring$

$\{\begin{array}{l} Variables & &colon; & [z, y, x] \\ Parameters & &colon; & \emptyset \\ Characteristic & &colon; & 0 \end{array}$

$p ≔ x^{3} y + x^{3} + x z + 4 x + 3$

$z$

$x$

$z$

$x^{3} y + x^{3} + 4 x + 3$

$x$

(2.1.14)

Consider the polynomial ring ${&integers;}_{3} [z, y, x]$ with the ordering $z > y > x$ .

>	R := PolynomialRing([z, y, x], 3); Display(R); Initial(p, R); Tail(p, R);

$R ≔ polynomial_ring$

$\{\begin{array}{l} Variables & &colon; & [z, y, x] \\ Parameters & &colon; & \emptyset \\ Characteristic & &colon; & 3 \end{array}$

$x$

$x^{3} y + x^{3} + x$

(2.1.15)

2.2 Computing the subresultants of two polynomials

The command SubresultantChain in the ChainTools module is used to compute the subresultant chain of two polynomials. Let us illustrate this by an example.

Define a ring of polynomials.

>	R := PolynomialRing([y, x]);

$R ≔ polynomial_ring$

(2.2.1)

Define two polynomials of R.

>	f1 := (y^2+6)(x-1)-y(x^2+1);

$f1 ≔ (y^{2} + 6) (x - 1) - y (x^{2} + 1)$

(2.2.2)

>	f2 := (x^2+6)(y-1)-x(y^2+1);

$f2 ≔ (x^{2} + 6) (y - 1) - x (y^{2} + 1)$

(2.2.3)

Compute their subresultant chain w.r.t the variable $y$

>	src := SubresultantChain(f1, f2, y, R);

$src ≔ subresultant_chain$

(2.2.4)

The result is of type $subresultant_chain$ and encodes all the polynomials in the subresultant chain. To display the $i_{th}$ subresulant polynomial, one can call the command SubresultantOfIndex. For example, to know the subresultant of index $0,$ that is the resultant, one can do as follows

>	res := SubresultantOfIndex(0, src, R);

$res ≔ 2 x^{6} - 22 x^{5} + 102 x^{4} - 274 x^{3} + 488 x^{2} - 552 x + 288$

(2.2.5)

To display all the subresultant polynomials, one can call the command Display.

>	Display(src, R);

$\{\begin{array}{c} (- x^{2} + 5 x - 6) y - x^{3} + 6 x^{2} - 11 x + 6 \\ 2 x^{6} - 22 x^{5} + 102 x^{4} - 274 x^{3} + 488 x^{2} - 552 x + 288 \end{array}$

(2.2.6)

The encoding of the subresultants can be controlled by an optional argument. This encoding can be by values, over the monomial basis or via the Bezout Matrix.

When this optional argument is not specified, the encoding is chosen according to a heuristical algorithm.

>	src := SubresultantChain(f1, f2, y, R, representation = BezoutMatrix); op(src);

$src ≔ subresultant_chain$

$table ([SRC_MDEG = 1, SRC_POLYQ = x^{2} y - x y^{2} - x^{2} - x + 6 y - 6, representation = BezoutMatrix, SRC_MVAR = y, type = subresultant_chain, SRC_MATRIX = [\begin{array}{c} - x^{2} + 5 x - 6 & - x^{3} + 6 x^{2} - 11 x + 6 \\ - x^{3} + 6 x^{2} - 11 x + 6 & x^{4} - 5 x^{3} + 13 x^{2} - 35 x + 42 \end{array}], SRC_POLYP = - x^{2} y + x y^{2} - y^{2} + 6 x - y - 6])$

(2.2.7)

>	src := SubresultantChain(f1, f2, y, R, representation = MonomialBasis); op(src);

$src ≔ subresultant_chain$

$table ([SRC_MDEG = 1, SRC_POLYQ = x^{2} y - x y^{2} - x^{2} - x + 6 y - 6, representation = MonomialBasis, SRC_MVAR = y, subresultant_chain_vector = [2 x^{6} - 22 x^{5} + 102 x^{4} - 274 x^{3} + 488 x^{2} - 552 x + 288, - x^{3} - x^{2} y + 6 x^{2} + 5 x y - 11 x - 6 y + 6, x^{2} y - x y^{2} - x^{2} - x + 6 y - 6, - x^{2} y + x y^{2} - y^{2} + 6 x - y - 6], type = subresultant_chain, SRC_POLYP = - x^{2} y + x y^{2} - y^{2} + 6 x - y - 6])$

(2.2.8)

2.3 Solving systems of equations with regular chains

Define a ring of polynomials.

>	R := PolynomialRing([x, y, z]);

$R ≔ polynomial_ring$

(2.3.1)

Define a set of polynomials of $R$ .

>	sys := {x+y+z^2-1, x+y^2+z-1, x^2+y+z-1};

$sys ≔ \{z^{2} + x + y - 1, y^{2} + x + z - 1, x^{2} + y + z - 1\}$

(2.3.2)

Ideally, one would like to decompose the solutions of sys into a list of points. This is what the command Triangularize does using symbolic expressions. In the output decomposition, some points are grouped because they share some properties. These groups are called regular chains.

A regular chain is a system of equations and inequations, satisfying special algebraic properties. First, the equation part consists of non-constant polynomials with pairwise different main variables; thus, the equation part has a triangular shape. Second, the inequation part consists of the initials of the polynomials defining the equations; moreover, the product of those initials is regular (that is, not a zero-divisor) modulo the saturated ideal of the equation part. This technical condition can be skipped by the non-expert reader, since in most practical examples the inequations are trivial and can be ignored. This is the case in our example where all initials are equal to 1, leading to trivial inequations.

>	l := Triangularize(sys, R);

$l ≔ [regular_chain, regular_chain, regular_chain, regular_chain]$

(2.3.3)

>	map(Equations, l, R); map(Inequations, l, R); map(NumberOfSolutions, l, R);

$[[x - z, y - z, z^{2} + 2 z - 1], [x, y, z - 1], [x, y - 1, z], [x - 1, y, z]]$

$[\emptyset, \emptyset, \emptyset, \emptyset]$

$[2, 1, 1, 1]$

(2.3.4)

The result is to be interpreted as follows: The system $sys$ is equivalent to a union of four regular chains. The equations for the first such regular chain are $x - 1 = y = z = 0$ , the equations for the second regular chain are $x = y - 1 = z = 0$ , etc. None of the four regular chains have any inequations. The first three regular chains have exactly one solution, and the last one has two solutions.

Note that you can also specify inequations of the form $p \neq 0$ . For example, you can specify that $x - z$ must not vanish.

>	l := Triangularize(sys, [x-z], R); map(Equations, l, R); map(NumberOfSolutions, l, R);

$l ≔ [regular_chain, regular_chain]$

$[[x - 1, y, z], [x, y, z - 1]]$

$[1, 1]$

(2.3.5)

2.4 Solving polynomial systems with infinitely many solutions

In the previous examples, the polynomial systems have finitely many solutions. To illustrate how the Triangularize command handles systems with infinitely many solutions, consider the following parametric linear system with unknowns $x, y$ and parameters $a, b, c, d, g, h$ .

>	R := PolynomialRing([x, y, a, b, c, d, g, h]); sys := {ax+by-g, cx+dy-h};

$R ≔ polynomial_ring$

$sys ≔ \{a x + b y - g, c x + d y - h\}$

(2.4.1)

By default, the Triangularize command computes the generic solutions of the input system. For this parametric linear system, this implies its determinant is assumed to be nonzero.

>	l := Triangularize(sys, R); map(Equations, l, R); map(Inequations, l, R);

$l ≔ [regular_chain]$

$[[c x + y d - h, (a d - b c) y - a h + c g]]$

$[\{c, d a - b c\}]$

(2.4.2)

Thus, the resulting regular chain consists of the two equations $c x + d y - h = (a d - b c) y + c g - a h = 0$ and the two inequations $c \neq 0 \neq a d - b c$ , the first of which reflects the choice of $c$ as a pivot.

The other available Maple commands for solving such a system do a similar job.

>	solve(sys, {x, y});

$\{x = - \frac{b h - d g}{d a - b c}, y = \frac{a h - c g}{d a - b c}\}$

(2.4.3)

>	Groebner[Solve](sys, {x, y});

$\{[[d x a - b c x + b h - d g, y a d - y b c - a h + c g], plex (y, x), \emptyset]\}$

(2.4.4)

The Triangularize command can do more if the option "output=lazard" is used. In this case, all the solutions of the input are computed, generic or not. This implies that the case where the determinant vanishes is also considered and solved. For each regular chain below, its equations and inequations are printed.

>	l := Triangularize(sys, R, output = lazard): seq(Display(l[i], R), i = 1..nops(l));

$\{\begin{array}{c} c x + d y - h = 0 \\ (d a - b c) y - h a + c g = 0 \\ d a - b c \neq 0 \\ c \neq 0 \end{array}, \{\begin{array}{c} c x + d y - h = 0 \\ d a - b c = 0 \\ h b - d g = 0 \\ c \neq 0 \\ d \neq 0 \\ h \neq 0 \end{array}, \{\begin{array}{c} a x + b y - g = 0 \\ d y - h = 0 \\ c = 0 \\ a \neq 0 \\ d \neq 0 \end{array}, \{\begin{array}{c} d y - h = 0 \\ a = 0 \\ h b - d g = 0 \\ c = 0 \\ d \neq 0 \\ h \neq 0 \end{array}, \{\begin{array}{c} c x - h = 0 \\ h a - c g = 0 \\ b = 0 \\ d = 0 \\ c \neq 0 \\ h \neq 0 \end{array}, \{\begin{array}{c} a x + b y - g = 0 \\ c = 0 \\ d = 0 \\ h = 0 \\ a \neq 0 \end{array}, \{\begin{array}{c} c x + d y = 0 \\ d a - b c = 0 \\ g = 0 \\ h = 0 \\ c \neq 0 \\ d \neq 0 \end{array}, \{\begin{array}{c} b y - g = 0 \\ a = 0 \\ c = 0 \\ d = 0 \\ h = 0 \\ b \neq 0 \end{array}, \{\begin{array}{c} y = 0 \\ a = 0 \\ c = 0 \\ g = 0 \\ h = 0 \end{array}, \{\begin{array}{c} x = 0 \\ b = 0 \\ d = 0 \\ g = 0 \\ h = 0 \end{array}, \{\begin{array}{c} a = 0 \\ b = 0 \\ c = 0 \\ d = 0 \\ g = 0 \\ h = 0 \end{array}$

(2.4.5)

One may also want to "push" some parameters into the field of coefficients, that is, to solve over a field of rational functions. Consider a new polynomial ring over the field $&rationals; (g, h)$ of rational functions in $g$ and $h$ .

>	$R_{2} ≔ PolynomialRing ([x, y, a, b, c, d], \{g, h\}) &colon;$ $l ≔ Triangularize (sys, R_{2}, output = lazard) &colon;$ $seq (Display (l_{i}, R_{2}), i = 1 .. nops (l)) &semi;$

$\{\begin{array}{c} c x + d y - h = 0 \\ (d a - b c) y - h a + c g = 0 \\ d a - b c \neq 0 \\ c \neq 0 \end{array}, \{\begin{array}{c} c x + d y - h = 0 \\ d a - b c = 0 \\ h b - d g = 0 \\ c \neq 0 \\ d \neq 0 \end{array}, \{\begin{array}{c} a x + y b - g = 0 \\ d y - h = 0 \\ c = 0 \\ a \neq 0 \\ d \neq 0 \end{array}, \{\begin{array}{c} d y - h = 0 \\ a = 0 \\ h b - d g = 0 \\ c = 0 \\ d \neq 0 \end{array}, \{\begin{array}{c} c x - h = 0 \\ h a - c g = 0 \\ b = 0 \\ d = 0 \\ c \neq 0 \end{array}$

(2.4.6)

Finally, we solve over $ℤ_{3} (g, h)$ , that is, over the field of rational functions in $g$ and $h$ modulo 3.

>	$R_{3} := PolynomialRing ([x, y, a, b, c, d], \{g, h\}, 3) &semi;$ $l := Triangularize (sys, R_{3}, output = lazard) &semi;$ $op (Display (l, R_{3})) &semi;$

$R_{3} ≔ polynomial_ring$

$l ≔ [regular_chain, regular_chain, regular_chain, regular_chain, regular_chain]$

$\{\begin{array}{c} c x + d y + 2 h = 0 \\ (d a + 2 b c) y + 2 h a + c g = 0 \\ d a + 2 b c \neq 0 \\ c \neq 0 \end{array}, \{\begin{array}{c} c x + d y + 2 h = 0 \\ d a + 2 b c = 0 \\ h b + 2 d g = 0 \\ c \neq 0 \\ d \neq 0 \end{array}, \{\begin{array}{c} a x + y b + 2 g = 0 \\ d y + 2 h = 0 \\ c = 0 \\ a \neq 0 \\ d \neq 0 \end{array}, \{\begin{array}{c} d y + 2 h = 0 \\ a = 0 \\ h b + 2 d g = 0 \\ c = 0 \\ d \neq 0 \end{array}, \{\begin{array}{c} c x + 2 h = 0 \\ h a + 2 c g = 0 \\ b = 0 \\ d = 0 \\ c \neq 0 \end{array}$

(2.4.7)

2.5 Regular chains and polynomial gcds

The main subroutine of the Triangularize command is the RegularGcd command, which computes polynomial gcds modulo regular chains. This routine is illustrated hereafter together with basic commands on regular chains from the ChainTools submodule.

>	R := PolynomialRing([y, x]);

$R ≔ polynomial_ring$

(2.5.1)

Here's the empty regular chain and its internal representation.

>	Empty(R); op(Empty(R));

$regular_chain$

$property, polynomials, type, ModulePrint, module () option record &semi; export property, polynomials, type, ModulePrint &semi; end module$

(2.5.2)

You can construct a new regular chain by adding a (suitable) polynomial $p$ to an existing regular chain $T$ . This operation may return a list of regular chains. Indeed, checking that the initial of $p$ is regular with regards to (the saturated ideal of) $T$ may split $T$ .

On the example below, this does not happen; the regular chain $T$ is empty.

>	l := Construct(x^2+1, Empty(R), R); rc := l[1];

$l ≔ [regular_chain]$

$rc ≔ regular_chain$

(2.5.3)

The internal representation of the new regular chain $rc$ is as follows. For each polynomial in the regular chain, the data structure stores its main variable, its main degree, and its initial. These quantities are needed quite frequently and are therefore cached to avoid recomputing them each time they are needed.

The internal representation also shows a property about the regular chain: isPrime. This means that its saturated ideal is a prime ideal. The knowledge of such property helps to speed up computations.

op(rc);

$regular_chain$

(2.5.4)

Now compute the greatest common divisor of two polynomials $p_{1}$ and $p_{2}$ modulo $rc$ . The example is designed such that $y + 1$ is a gcd.

>	$p_{1} := expand (((x + 1) \cdot y + x) \cdot (y + 1)) &semi;$ $p_{2} ≔ expand (((x - 1) \cdot y + x) \cdot (y + 1)) &semi;$

$p_{1} ≔ x y^{2} + 2 x y + y^{2} + x + y$

$p_{2} ≔ x y^{2} + 2 x y - y^{2} + x - y$

(2.5.5)

In general, the output of a gcd computation modulo a regular chain is a list of "cases".

In the above example, no case splits can happen since $rc$ defines a field. However, note that the output, below, is not exactly the one expected.

>	$rg := RegularGcd (p_{1}, p_{2}, y, rc, R) &semi;$ $factor ({rg}_{1}_{1})$

$rg ≔ [[2 x y + 2 x, regular_chain]]$

$2 x (y + 1)$

(2.5.6)

However, the extraneous factor $2 x$ is a unit modulo $rc$ . Thus, the expected and the computed gcds are in fact associate, and therefore the result is correct.

Consider an example where a split is necessary during a gcd computation. To create such an example, you need to make sure that the input equations do not factor, because otherwise the Construct command will already split the system. The example below has three variables; the two variables $x$ and $y$ are used for building a regular chain with two polynomials $mx$ and $my$ .

>	R := PolynomialRing([z, y, x]);

$R ≔ polynomial_ring$

(2.5.7)

The first polynomial is irreducible over <.

>	mx := x^2-2;

$mx ≔ x^{2} - 2$

(2.5.8)

Thus, constructing a regular chain from it produces a single output. Make a copy of it for later use.

>	rc := Construct(mx, Empty(R), R)[1]; rcx := rc;

$rc ≔ regular_chain$

$rcx ≔ rc$

(2.5.9)

The Construct command has detected that (the saturated ideal of) $rc$ is a prime ideal.

op(rc);

$regular_chain$

(2.5.10)

The second polynomial is not irreducible in $&rationals; [x] [y]$ . The Construct command would easily discover this fact and would split it when constructing a regular chain involving $rc$ .

>	my := expand((y-x)*(y+x-1));

$my ≔ - x^{2} + y^{2} + x - y$

(2.5.11)

By reducing $my$ modulo $rc$ , the reduced polynomial becomes irreducible in $ℚ [x] [y]$ .

>	my := SparsePseudoRemainder(my, rc, R);

$my ≔ y^{2} + x - y - 2$

(2.5.12)

Extend $rc$ by $my$ .

>	l := Construct(my, rc, R);rc := l[1];

$l ≔ [regular_chain]$

$rc ≔ regular_chain$

(2.5.13)

The Construct command uses inexpensive criteria to detect whether the associated saturated ideal is prime or not. These criteria fail here. So the constructed regular chain $rc$ has no particular properties.

>	op(rc);Equations(rc, R);

$regular_chain$

$[y^{2} - y + x - 2, x^{2} - 2]$

(2.5.14)

Compute the gcd of $p_{1}$ and $p_{2}$ , defined below, modulo $rc$ .

>	$p_{1} ≔ (y - x) z + (x + y - 1) (z + 1) &semi;$ $p_{2} ≔ y - x + (x + y - 1) (z + 1)$

$p_{1} ≔ (y - x) z + (y + x - 1) (z + 1)$

$p_{2} ≔ y - x + (y + x - 1) (z + 1)$

(2.5.15)

The example is constructed in such a way that splitting is needed: modulo the first factor of $my$ , namely $y - x$ , the gcd is $z + 1$ , while it is constant modulo the second factor, $y + x - 1$ . Correspondingly, there are two branches.

>	$rg := RegularGcd (p_{1}, p_{2}, z, rc, R) &semi;$

$rg ≔ [[(y + x - 1) z + 2 y - 1, regular_chain], [3 y^{2} + (- 2 x - 2) y - x^{2} + 2 x, regular_chain]]$

(2.5.16)

Check the first case.

>	$g_{1}, rc1 ≔ op ({rg}_{1}) &semi;$ ${eq}_{1} := Equations (rc1, R) &semi;$ $SparsePseudoRemainder ((2 x + 1) {eq}_{1}_{1}, rcx, R)$

$g_{1}, rc1 ≔ (y + x - 1) z + 2 y - 1, regular_chain$

${eq}_{1} ≔ [(2 x - 1) y + x - 4, x^{2} - 2]$

$- 7 x + 7 y$

(2.5.17)

Note that ${rc}_{1}$ is not normalized. The SparsePseudoRemainder computation shows that it corresponds to the "modulo $y - x$ " branch.

>	$u := SparsePseudoRemainder (g_{1}, rc1, R) &semi;$ $factor (u)$

$u ≔ - 4 z x - 4 x + 9 z + 9$

$- (z + 1) (4 x - 9)$

(2.5.18)

This is the desired result, up to a multiplicative factor $- 4 x + 9$ which is a unit modulo ${rc}_{1}$ .

Check the second case. Again, the regular chain ${rc}_{2}$ is not normalized, and the SparsePseudoRemainder computation shows that it corresponds to the "modulo $y + x - 1$ " branch.

>	$g_{2}, rc2 ≔ op ({rg}_{2}) &semi;$ ${eq}_{2} ≔ Equations (rc2, R) &semi;$ $SparsePseudoRemainder ((2 x + 1) {eq}_{2}_{1}, rcx, R)$

$g_{2}, rc2 ≔ 3 y^{2} + (- 2 x - 2) y - x^{2} + 2 x, regular_chain$

${eq}_{2} ≔ [(2 x - 1) y - 3 x + 5, x^{2} - 2]$

$7 x + 7 y - 7$

(2.5.19)

Check that $g_{2}$ is a unit modulo ${rc}_{2}$ . First, reduce it modulo ${rc}_{2}$ . Then ask if the result is invertible modulo ${rc}_{2}$ . The answer proves it is.

>	$h := SparsePseudoRemainder (g_{2}, rc2, R) &semi;$ $out := Inverse (h, rc2, R) &semi;$ $ih := {out}_{1}_{1}_{1} &semi;$ $SparsePseudoRemainder (h ih, rc2, R)$

$h ≔ - 72 x + 113$

$out ≔ [[[113 + 72 x, 2401, regular_chain]], []]$

$ih ≔ 113 + 72 x$

$2401$

(2.5.20)

2.6 Solving systems of equations incrementally

The central command of the RegularChains library is the Triangularize command. The algorithm behind this command solves systems of equations incrementally. To do this, we rely on an important operation Intersect, which computes the common part of a hypersurface and the quasi-component of a regular chain. See the help page of Intersect for a more formal description of it. Let us now illustrate how to use it to solve a polynomial system incrementally.

Define a ring of polynomials.

>	vars := [x, y, z];R := PolynomialRing(vars);

$vars ≔ [x, y, z]$

$R ≔ polynomial_ring$

(2.6.1)

Define a set of equations.

>	sys := [x^2+y+z-1, x+y^2+z-1, x+y+z^2-1];

$sys ≔ [x^{2} + y + z - 1, y^{2} + x + z - 1, z^{2} + x + y - 1]$

(2.6.2)

Define the empty regular chain.

>	rc := Empty(R);

$rc ≔ regular_chain$

(2.6.3)

Solve the first equation.

>	dec := Intersect(sys[1], rc, R);map(Equations, dec, R);

$dec ≔ [regular_chain]$

$[[x^{2} + y + z - 1]]$

(2.6.4)

Solve the first and second equations.

>	dec := [seq(op(Intersect(sys[2], rc, R)), `in`(rc, dec))];map(Equations, dec, R);

$dec ≔ [regular_chain, regular_chain]$

$[[x - y, y^{2} + y + z - 1], [x + y - 1, y^{2} - y + z]]$

(2.6.5)

Solve the three equations together.

>	dec := [seq(op(Intersect(sys[3], rc, R)), `in`(rc, dec))];Display(dec, R);

$dec ≔ [regular_chain, regular_chain, regular_chain, regular_chain]$

$[\{\begin{array}{c} x - z = 0 \\ y - z = 0 \\ z^{2} + 2 z - 1 = 0 \end{array}, \{\begin{array}{c} x = 0 \\ y = 0 \\ z - 1 = 0 \end{array}, \{\begin{array}{c} x - 1 = 0 \\ y = 0 \\ z = 0 \end{array}, \{\begin{array}{c} x = 0 \\ y - 1 = 0 \\ z = 0 \end{array}]$

(2.6.6)

3. Real solutions of polynomial systems

The RegularChains library offers a variety of tools to compute the real solutions of polynomial systems. A first group of commands (namely RealTriangularize, LazyTriangularize and SamplePoints) deals with arbitrary semi-algebraic systems. That is, given any system $S$ of polynomial equations, polynomial inequations and polynomial inequalities (strict or large) these commands produce information about the real solutions of this system.

RealTriangularize returns a full description of the real solutions of $S$ : it computes simpler systems $S_{1}, .. &period;, S_{e}$ such that a point is a solution of $S$ if and only if it is a solution of one of the systems $S_{1}, .. &period;, S_{e}$ . Each of these systems

has a triangular shape and remarkable properties: for this reason it is called a regular semi-algebraic system and the set of the $S_{1}, .. &period;, S_{e}$ is called a full triangular decomposition of $S$ .

LazyRealTriangularize allows the user to compute a triangular decomposition of $S$ in an interactive manner. This feature is particularly well adapted for systems that are hard to solve. For such systems, LazyRealTriangularize returns the components of $S$ of maximum dimension together with unevaluated recursive calls, such that, when fully evaluated, these calls produce the other components of $S$ (which are generally harder to compute).

SamplePoints is even a lazier (and thus much cheaper) way of solving:

it produces at least one sample point per connected component

of the solution set of $S$ . This way of solving is often sufficient in practical problems.

A second group of commands compute the real solutions of particular

types of polynomial systems (such as systems with finitely many complex solutions or parametric systems) or provide advanced features (such as CylindricalAlgebraicDecompose, LinearSolve, Projection, Difference). All these commands except RealRootClassification and RealComprehensiveTriangularize are located in the subpackage SemiAlgebraicSetTools.

3.1 RealTriangularize: solving systems of equations, inequations and inequalities

Consider the generic equation of degree two.

>	R := PolynomialRing([x, c, b, a]); sys := [ax^2+bx+c=0];

$R ≔ polynomial_ring$

$sys ≔ [a x^{2} + b x + c = 0]$

(3.1.1)

Compute a triangular decomposition of the 4-variable hypersurface it defines.

>	dec := RealTriangularize(sys, R);

$dec ≔ [regular_semi_algebraic_system, regular_semi_algebraic_system, regular_semi_algebraic_system, regular_semi_algebraic_system]$

(3.1.2)

>	Display(dec, R);

$[\{\begin{array}{c} a x^{2} + b x + c = 0 \\ - 4 c a + b^{2} > 0 and a \neq 0 \end{array}, \{\begin{array}{c} 2 a x + b = 0 \\ 4 a c - b^{2} = 0 \\ a \neq 0 \end{array}, \{\begin{array}{c} b x + c = 0 \\ a = 0 \\ b \neq 0 \end{array}, \{\begin{array}{c} c = 0 \\ b = 0 \\ a = 0 \end{array}]$

(3.1.3)

Consider the record output format.

>	RealTriangularize(sys, R, output=record);

$\{\begin{array}{c} a x^{2} + b x + c = 0 \\ - 4 c a + b^{2} > 0 \\ a \neq 0 \end{array}, \{\begin{array}{c} 2 a x + b = 0 \\ 4 a c - b^{2} = 0 \\ a \neq 0 \end{array}, \{\begin{array}{c} b x + c = 0 \\ b \neq 0 \\ a = 0 \end{array}, \{\begin{array}{c} c = 0 \\ b = 0 \\ a = 0 \end{array}$

(3.1.4)

Next, we consider a system of equations, inequations and inequalities.

>	R := PolynomialRing([y, x, b, a]);

$R ≔ polynomial_ring$

(3.1.5)

>	sys := [x^3-3y^2x+ax+b=0, 3x^2-y^2+a=0, 1-y*x>0, y<>0];

$sys ≔ [x^{3} - 3 y^{2} x + a x + b = 0, 3 x^{2} - y^{2} + a = 0, 0 < - y x + 1, y \neq 0]$

(3.1.6)

>	out := RealTriangularize(sys, R);

$out ≔ [regular_semi_algebraic_system, regular_semi_algebraic_system]$

(3.1.7)

>	Display(out, R);

$[\{\begin{array}{c} y^{2} - 3 x^{2} - a = 0 \\ 8 x^{3} + 2 a x - b = 0 \\ - y x + 1 > 0 \\ 4 a^{3} + 27 b^{2} > 0 and 4 b^{2} a^{3} - 16 a^{4} + 27 b^{4} - 512 a^{2} - 4096 \neq 0 \end{array}, \{\begin{array}{c} x y + 1 = 0 \\ (2 a^{3} + 18 b^{2} + 32 a) x + (- a^{2} - 48) b = 0 \\ 27 b^{4} + 4 a^{3} b^{2} - 16 a^{4} - 512 a^{2} - 4096 = 0 \end{array}]$

(3.1.8)

Consider now an example which has finitely many complex solutions.

>	R := PolynomialRing([x, y, z]):

Define a set of equations.

>	sys := [x^3 + y + z -1=0, x + y^3 + z -1=0, x + y + z^3 -1=0];

$sys ≔ [x^{3} + y + z - 1 = 0, y^{3} + x + z - 1 = 0, z^{3} + x + y - 1 = 0]$

(3.1.9)

Compute the real solutions of sys.

>	dec := RealTriangularize(sys, R);

$dec ≔ [regular_semi_algebraic_system, regular_semi_algebraic_system, regular_semi_algebraic_system, regular_semi_algebraic_system, regular_semi_algebraic_system, regular_semi_algebraic_system, regular_semi_algebraic_system]$

(3.1.10)

>	Display(dec, R);

$[\{\begin{array}{c} x - z = 0 \\ y - z = 0 \\ z^{3} + 2 z - 1 = 0 \end{array}, \{\begin{array}{c} x - 1 = 0 \\ y + 1 = 0 \\ z - 1 = 0 \end{array}, \{\begin{array}{c} x = 0 \\ y = 0 \\ z - 1 = 0 \end{array}, \{\begin{array}{c} x + 1 = 0 \\ y - 1 = 0 \\ z - 1 = 0 \end{array}, \{\begin{array}{c} x - 1 = 0 \\ y = 0 \\ z = 0 \end{array}, \{\begin{array}{c} x = 0 \\ y - 1 = 0 \\ z = 0 \end{array}, \{\begin{array}{c} x - 1 = 0 \\ y - 1 = 0 \\ z + 1 = 0 \end{array}]$

(3.1.11)

Returning now to systems which have infinitely many complex solutions, consider now the intersection of the algebraic surfaces Sofa and Cylinder from the Algebraic Surface Gallery. As their names suggest, they are respectively the equations of a sofa and a (sort of) cylinder. One expects to find a real curve in their intersection, as we shall verify.

>	R := PolynomialRing([z, y, x]); Sofa := x^2 + y^3 + z^5; Cyl := x^4 + z^2 - 1;

$R ≔ polynomial_ring$

$Sofa ≔ z^{5} + y^{3} + x^{2}$

$Cyl ≔ x^{4} + z^{2} - 1$

(3.1.12)

>	RealTriangularize([Sofa, Cyl], R, output=record);

$\{\begin{array}{c} (x^{8} - 2 x^{4} + 1) z + y^{3} + x^{2} = 0 \\ y^{6} + 2 x^{2} y^{3} + x^{20} - 5 x^{16} + 10 x^{12} - 10 x^{8} + 6 x^{4} - 1 = 0 \\ x < 1 \\ x + 1 > 0 \\ x^{12} - 4 x^{8} + 5 x^{4} - 1 \neq 0 \end{array}, \{\begin{array}{c} (x^{8} - 2 x^{4} + 1) z - x^{2} = 0 \\ y^{3} + 2 x^{2} = 0 \\ x^{12} - 4 x^{8} + 5 x^{4} - 1 = 0 \end{array}, \{\begin{array}{c} (x^{8} - 2 x^{4} + 1) z + x^{2} = 0 \\ y = 0 \\ x^{12} - 4 x^{8} + 5 x^{4} - 1 = 0 \end{array}, \{\begin{array}{c} z = 0 \\ y + 1 = 0 \\ x + 1 = 0 \end{array}, \{\begin{array}{c} z = 0 \\ y + 1 = 0 \\ x - 1 = 0 \end{array}$

(3.1.13)

3.2 LazyRealTriangularize: interactive solving for systems of equations, inequations and inequalities

Consider again the generic equation of degree two. Now we solve it interactively.

>	R := PolynomialRing([x, c, b, a]); sys := [ax^2+bx+c=0];

$R ≔ polynomial_ring$

$sys ≔ [x^{2} a + x b + c = 0]$

(3.2.1)

Use LazyRealTriangularize to start the decomposition.

>	dec := LazyRealTriangularize(sys, R);

$dec ≔ \{\begin{array}{c} [[a x^{2} + b x + c = 0]] & 0 < - 4 c a + b^{2} \land a \neq 0 \\ LazyRealTriangularize ([a = 0, a x^{2} + b x + c = 0], polynomial_ring) & a = 0 \\ LazyRealTriangularize ([- 4 c a + b^{2} = 0, a x^{2} + b x + c = 0], polynomial_ring) & - 4 c a + b^{2} = 0 \\ [] & otherwise \end{array}$

(3.2.2)

Go one step further, computing components in lower dimension.

>	dec2 := value(dec);

$dec2 ≔ \{\begin{array}{c} [[x^{2} a + x b + c = 0]] & 0 < - 4 a c + b^{2} \land a \neq 0 \\ \{\begin{array}{c} [[x b + c = 0, a = 0]] & b \neq 0 \\ LazyRealTriangularize ([a = 0, b = 0, x b + c = 0, x^{2} a + x b + c = 0], polynomial_ring) & b = 0 \end{array} & a = 0 \\ \{\begin{array}{c} [[2 a x + b = 0, 4 a c - b^{2} = 0]] & a \neq 0 \\ LazyRealTriangularize ([a = 0, - 4 a c + b^{2} = 0, 4 a c - b^{2} = 0, 2 a x + b = 0, x^{2} a + x b + c = 0], polynomial_ring) & a = 0 \end{array} & - 4 a c + b^{2} = 0 \\ [] & otherwise \end{array}$

(3.2.3)

Go the last step, computing components in dimension zero.

>	value(dec2);

$\{\begin{array}{c} [[x^{2} a + x b + c = 0]] & 0 < - 4 a c + b^{2} \land a \neq 0 \\ \{\begin{array}{c} [[x b + c = 0, a = 0]] & b \neq 0 \\ [[c = 0, b = 0, a = 0]] & b = 0 \end{array} & a = 0 \\ \{\begin{array}{c} [[2 a x + b = 0, 4 a c - b^{2} = 0]] & a \neq 0 \\ [[c = 0, b = 0, a = 0]] & a = 0 \end{array} & - 4 a c + b^{2} = 0 \\ [] & otherwise \end{array}$

(3.2.4)

If one is only interested in computing the main components, the list output format does the job.

>	dec := LazyRealTriangularize(sys, R, output=list);

$dec ≔ [regular_semi_algebraic_system]$

(3.2.5)

>	Display(dec, R);

$[\{\begin{array}{c} a x^{2} + b x + c = 0 \\ - 4 c a + b^{2} > 0 and a \neq 0 \end{array}]$

(3.2.6)

Another way to conduct the computation interactively is to use the record output format.

>	dec := [LazyRealTriangularize(sys, R, output=record)];

$dec ≔ [\{\begin{array}{c} a x^{2} + b x + c = 0 \\ - 4 c a + b^{2} > 0 \\ a \neq 0 \end{array}, LazyRealTriangularize ([a = 0, a x^{2} + b x + c = 0], polynomial_ring, output = record), LazyRealTriangularize ([- 4 c a + b^{2} = 0, a x^{2} + b x + c = 0], polynomial_ring, output = record)]$

(3.2.7)

Go one step further.

>	dec2 := value(dec);

$dec2 ≔ [\{\begin{array}{c} a x^{2} + b x + c = 0 \\ - 4 c a + b^{2} > 0 \\ a \neq 0 \end{array}, \{\begin{array}{c} b x + c = 0 \\ b \neq 0 \\ a = 0 \end{array}, LazyRealTriangularize ([a = 0, b = 0, b x + c = 0, a x^{2} + b x + c = 0], polynomial_ring, output = record), \{\begin{array}{c} 2 a x + b = 0 \\ 4 a c - b^{2} = 0 \\ a \neq 0 \end{array}, LazyRealTriangularize ([a = 0, - 4 a c + b^{2} = 0, 4 a c - b^{2} = 0, 2 a x + b = 0, a x^{2} + b x + c = 0], polynomial_ring, output = record)]$

(3.2.8)

Go the last step.

>	value(dec2);

$[\{\begin{array}{c} a x^{2} + b x + c = 0 \\ - 4 a c + b^{2} > 0 \\ a \neq 0 \end{array}, \{\begin{array}{c} b x + c = 0 \\ b \neq 0 \\ a = 0 \end{array}, \{\begin{array}{c} c = 0 \\ b = 0 \\ a = 0 \end{array}, \{\begin{array}{c} 2 a x + b = 0 \\ 4 a c - b^{2} = 0 \\ a \neq 0 \end{array}, \{\begin{array}{c} c = 0 \\ b = 0 \\ a = 0 \end{array}]$

(3.2.9)

Computing a lazy triangular decomposition is usually much less expensive than computing a full one. The following is an example.

>	unassign('u'); variables := [x, u, v, w]; R := PolynomialRing(variables); sys := [ux^2+vx+1=0, vx^3+wx+u=0, wx^2+vx+u<=0];

$variables ≔ [x, u, v, w]$

$R ≔ polynomial_ring$

$sys ≔ [u x^{2} + v x + 1 = 0, v x^{3} + w x + u = 0, w x^{2} + v x + u \leq 0]$

(3.2.10)

Computing a lazy decomposition takes than a second.

>	dec := LazyRealTriangularize(sys,R,output=list);

$dec ≔ [regular_semi_algebraic_system]$

(3.2.11)

Computing a full one does not terminate within an hour.

3.3 SamplePoints: computing at least one point per connected component

Consider again the generic equation of degree two.

>	R := PolynomialRing([x, c, b, a]); F := [ax^2+bx+c=0];

$R ≔ polynomial_ring$

$F ≔ [a x^{2} + b x + c = 0]$

(3.3.1)

Compute sample points of the 4-variable hypersurface it defines.

>	S := SamplePoints(F, R, output=record);

$S ≔ \{\begin{array}{c} x = −1 \\ c = \frac{1}{2} \\ b = 0 \\ a = - \frac{1}{2} \end{array}, \{\begin{array}{c} x = 1 \\ c = \frac{1}{2} \\ b = 0 \\ a = - \frac{1}{2} \end{array}, \{\begin{array}{c} x = −1 \\ c = - \frac{1}{2} \\ b = 0 \\ a = \frac{1}{2} \end{array}, \{\begin{array}{c} x = 1 \\ c = - \frac{1}{2} \\ b = 0 \\ a = \frac{1}{2} \end{array}, \{\begin{array}{c} x = 0 \\ c = 0 \\ b = - \frac{1}{2} \\ a = 0 \end{array}, \{\begin{array}{c} x = 0 \\ c = 0 \\ b = \frac{1}{2} \\ a = 0 \end{array}, \{\begin{array}{c} x = 0 \\ c = 0 \\ b = 0 \\ a = 0 \end{array}, \{\begin{array}{c} x = 0 \\ c = 0 \\ b = 0 \\ a = - \frac{1}{2} \end{array}, \{\begin{array}{c} x = 0 \\ c = 0 \\ b = 0 \\ a = \frac{1}{2} \end{array}$

(3.3.2)

Consider the Tacnode curve. We look for sample points in the middle of the right branch.

>	R := PolynomialRing([y, x]); F := [y^4-2y^3+y^2-3x^2y+2x^4];

$R ≔ polynomial_ring$

$F ≔ [2 x^{4} + y^{4} - 3 x^{2} y - 2 y^{3} + y^{2}]$

(3.3.3)

>	with(plots): implicitplot(F, x=-2..2, y=-1..3, numpoints=1000000);

>	SamplePoints([op(F), 2*x > 1, x < 1], R, output=record);

$\{\begin{array}{c} y = [\frac{13}{32}, \frac{105}{256}] \\ x = \frac{3}{4} \end{array}, \{\begin{array}{c} y = [\frac{237}{128}, \frac{119}{64}] \\ x = \frac{3}{4} \end{array}$

(3.3.4)

3.4 Isolating the real roots of a regular chain

Isolating real roots

Define a simple regular chain.

>	R := PolynomialRing([y, x]): rc := Chain([2x^2-1, 3y^3-x], Empty(R), R):

Isolate its real roots.

>	rr := RealRootIsolate(rc, R);

$rr ≔ [box, box]$

(3.4.1.1)

Display the box values.

>	solution := map(BoxValues, rr, R);

$solution ≔ [[y = [- \frac{20727073}{33554432}, - \frac{20727063}{33554432}], x = [- \frac{46341}{65536}, - \frac{741455}{1048576}]], [y = [\frac{10363531}{16777216}, \frac{10363537}{16777216}], x = [\frac{741455}{1048576}, \frac{46341}{65536}]]]$

(3.4.1.2)

Each element of the list is called a box. Each box encodes a real root, in the sense that it contains exactly one real root of $rc$ . The RealRootIsolate command returns all real roots of a given regular chain in this way. Thus, it is guaranteed that no root is lost. In the example, there are two real roots.

>	solution[1]; solution[2];

$[y = [- \frac{20727073}{33554432}, - \frac{20727063}{33554432}], x = [- \frac{46341}{65536}, - \frac{741455}{1048576}]]$

$[y = [\frac{10363531}{16777216}, \frac{10363537}{16777216}], x = [\frac{741455}{1048576}, \frac{46341}{65536}]]$

(3.4.1.3)

Thus, the first root satisfies $\frac{1}{2} < x < 1$ and $0 < y < 1$ , and the second one satisfies $- 1 < x < - \frac{1}{2}$ and $- 1 < x < 0$ .

The BoxValues command returns a list of the form $v = i$ , where $v$ is a variable and $i$ is either a rational number or an open interval encoded by a list. The reason why singletons are returned is that rational solutions are sometimes hit, as in the following example.

>	rc := Chain([2x-1, 3y^3-x], Empty(R), R): rr := RealRootIsolate(rc, R): map(BoxValues, rr, R);

$[[y = [\frac{577053}{1048576}, \frac{288527}{524288}], x = \frac{1}{2}]]$

(3.4.1.4)

The Display command prints the isolated real roots in a pretty manner.

>	Display(rr, R);

$[\{\begin{array}{c} y = [\frac{577053}{1048576}, \frac{288527}{524288}] \\ x = \frac{1}{2} \end{array}]$

(3.4.1.5)

Choosing a precision

The isolating boxes can be as small as needed. The abserr option allows you to obtain boxes whose widths are smaller than or equal to a certain absolute precision.

>	rc := Chain([2x^2-1, 3y^3-x], Empty(R), R): rr := RealRootIsolate(rc, R, abserr = 1/10^10): map(BoxValues, rr, R);

$[[y = [- \frac{2716737851577}{4398046511104}, - \frac{10866951406269}{17592186044416}], x = [- \frac{759250125}{1073741824}, - \frac{97184015999}{137438953472}]], [y = [\frac{2716737851567}{4398046511104}, \frac{2716737851577}{4398046511104}], x = [\frac{97184015999}{137438953472}, \frac{759250125}{1073741824}]]]$

(3.4.2.1)

evalf(%);

$[[y = [−0.6177146705, −0.6177146705], x = [−0.7071067812, −0.7071067812]], [y = [0.6177146705, 0.6177146705], x = [0.7071067812, 0.7071067812]]]$

(3.4.2.2)

In fact, this involves the algebraic numbers $\sqrt{\frac{1}{2}}$ and ${(\frac{\sqrt{\frac{1}{2}}}{3})}^{\frac{1}{3}}$ .

>	$evalf (\sqrt{\frac{1}{2}}) &semi;$ $evalf ({(\frac{\sqrt{\frac{1}{2}}}{3})}^{\frac{1}{3}}) &semi;$

$0.7071067810$

$0.6177146705$

(3.4.2.3)

3.5 Isolating and counting the real zeros of a semi-algebraic system with finitely many solutions

Recall that a semi-algebraic system (SAS) is a system containing polynomial equations, polynomial (non-strict or/and strict) inequalities, and polynomial inequations, which are denoted by $F$ , $N$ , $P$ , and $H$ , respectively. For example, the system $\{x^{2} + y^{2} - 1 = 0, x y - 1 = 0, x \geq 0, y > 0\}$ is represented by

>	F := [x^2+y^2-1, 2xy-1]: N := [x]: P := [y]: H := []:

Define the polynomial ring

>	R := PolynomialRing([x, y]):

The RealRootIsolate command isolates all the real zeros of a given semi-algebraic system.

>	rr := RealRootIsolate(F, N, P, H, R); solution := Display(rr, R);

$rr ≔ [box]$

$solution ≔ [\{\begin{array}{c} x = [\frac{5}{8}, \frac{3}{4}] \\ y = [\frac{45}{64}, \frac{91}{128}] \end{array}]$

(3.5.1)

The RealRootCounting command computes the number of distinct real solutions of a given semi-algebraic system.

>	RealRootCounting(F, N, P, H, R);

$1$

(3.5.2)

Consider a less trivial example.

>	R := PolynomialRing([z, y, x, c]): F := [1-cx-xy^2-xz^2, 1-cy-yx^2-yz^2, 1-cz-zx^2-zy^2, 8c^6+378*c^3-27]: N := []: P := [c, 1-c]: H := []: RealRootCounting(F, N, P, H, R);

$4$

(3.5.3)

If the input system has infinitely many complex solutions (that is, it is positive-dimensional over the complex numbers), RealRootCounting will return a message suggesting to call RealRootClassification instead.

>	R := PolynomialRing([x, y]): RealRootCounting([x^2+y^2], [], [], [], R);

Error, (in RegularChains:-SemiAlgebraicSetTools:-RealRootCounting) system is not zero-dimensional; try RealRootClassification |lib/RegularChains/src/user.mm:9796|

3.6 Partial cylindrical algebraic decomposition

The command RealRootClassification makes use of Partial Cylindrical Algebraic Decomposition (PCAD). The calling sequence is PartialCylindricalAlgebraicDecomposition( $p$ , $lp$ , $R$ ), where $R$ is a polynomial ring, $p$ is a polynomial in $R$ , and $lp$ is a list of polynomials in $R$ representing positivity conditions. The output of the function is a list $sp$ of sample points. Each sample point is represented by a list of rational numbers, as many as there are variables in $R$ . Each inner list gives the coordinates of a sample point of a (real) open connected component of the (real) space decomposed by the equation $p = 0$ . Moreover, sample points which do not satisfy $q > 0$ for all polynomials $q \in lp$ are discarded.

>	R := PolynomialRing([x, y, t]): PartialCylindricalAlgebraicDecomposition(y, [], R); PartialCylindricalAlgebraicDecomposition(y, [x], R);PartialCylindricalAlgebraicDecomposition(xy, [], R);PartialCylindricalAlgebraicDecomposition(xy, [x], R);

$[[0, - \frac{1}{2}, 0], [0, \frac{1}{2}, 0]]$

$[[\frac{1}{2}, - \frac{1}{2}, 0], [\frac{1}{2}, \frac{1}{2}, 0]]$

$[[- \frac{1}{2}, - \frac{1}{2}, 0], [\frac{1}{2}, - \frac{1}{2}, 0], [- \frac{1}{2}, \frac{1}{2}, 0], [\frac{1}{2}, \frac{1}{2}, 0]]$

$[[\frac{1}{2}, - \frac{1}{2}, 0], [\frac{1}{2}, \frac{1}{2}, 0]]$

(3.6.1)

As a more advanced example, compute sample points for the polynomial system of equations $sys$ below which in addition make each polynomial in $sys$ positive.

>	R := PolynomialRing([x, y]): sys := [-x^2-y+1, x+y^2-1, -x^2+y^2]: $sp ≔ PartialCylindricalAlgebraicDecomposition (\prod_{j = 1}^{nops (sys)} {sys}_{j}, sys, R)$

$sp ≔ [[0, - \frac{17}{8}], [\frac{221}{512}, - \frac{287}{256}], [\frac{197}{512}, \frac{13}{16}]]$

(3.6.2)

Check the result.

>	for i to nops(sp) do print(eval(sys, [x = sp[i][1], y = sp[i][2]])) end do;

$[\frac{25}{8}, \frac{225}{64}, \frac{289}{64}]$

$[\frac{507191}{262144}, \frac{45121}{65536}, \frac{280635}{262144}]$

$[\frac{10343}{262144}, \frac{23}{512}, \frac{134247}{262144}]$

(3.6.3)

Plot these sample points.

>	with(plots): colors := [blue, green, red, brown]: curves := seq(implicitplot(sys[i], x = -3 .. 3, y = -4 .. 3, color = colors[i], numpoints = 5000), i = 1 .. 3): points := pointplot(sp, color = colors[4]): display([curves, points]);

Compute and plot sample points for the same system $sys$ but without any sign constraints.

>	$sp ≔ PartialCylindricalAlgebraicDecomposition (\prod_{j = 1}^{nops (sys)} {sys}_{j}, [], R) &colon;$ points := pointplot(sp, color = colors[4]): display([curves, points]);

3.7 Cylindrical algebraic decomposition

Cylindrical Algebraic Decomposition (CAD) is a fundamental and powerful tool for studying systems of equations, inequations and inequalities.

Our algorithm is different from the traditional algorithm of Collins. It first computes a cylindrical decomposition of the complex space, from which a CAD of the real space can be easily extracted.

Consider the hyperbola $xy - 1 = 0.$

>	R := PolynomialRing([y, x]); F := [y*x-1];

$R ≔ polynomial_ring$

$F ≔ [x y - 1]$

(3.7.1)

A cylindrical algebraic decomposition adapted to the polynomial $xy - 1$ can be computed by the command CylindricalAlgebraicDecompose as follows:

>	outcad := CylindricalAlgebraicDecompose(F, R, output=piecewise);

$outcad ≔ \{\begin{array}{c} \{\begin{array}{c} 1 & y < \frac{1}{x} \\ 1 & y = \frac{1}{x} \\ 1 & \frac{1}{x} < y \end{array} & x < 0 \\ 1 & x = 0 \\ \{\begin{array}{c} 1 & y < \frac{1}{x} \\ 1 & y = \frac{1}{x} \\ 1 & \frac{1}{x} < y \end{array} & 0 < x \end{array}$

(3.7.2)

The output CAD is described by a nested piecewise function. The outmost piecewise function is a function with three conditions $x < 0$ , $x = 0$ , and $0 < x$ .
Each of the conditions has a corresponding expression, which is again a piecewise function. The output could be read from top to bottom and from right to left.

One can see that the CAD consists of seven cells.

For example, $x < 0 and y < \frac{1}{x}$ describes one cell of the CAD, where $[regular_chain, [[- 1, - 1], [- 2, - 2]]]$

represents a sample point in this cell.

This sample point is represented by a regular chain and an isolating box such that inside this box there is one and only one root of this regular chain.

We plot the hyperbola and all the sample points of the CAD adapted to this hyperbola as follows.

>	with(plots): sp := [[-1, -2], [-1, -1], [-1, 0], [0, 0], [1, 0], [1, 1], [1, 2]]: points := pointplot(sp, color = blue): curve := implicitplot([x*y-1, x], x = -5 .. 5, y = -5 .. 5, color=[red, black]): display([curve, points]);

The piecewise format is good when the output has few cells. When many cells are present the 'output'='cadcell' and 'output'='rootof' formats are useful.

>	R := PolynomialRing([x, c, b, a]); F := [ax^2+bx+c]; cad := CylindricalAlgebraicDecompose(F, R, output = cadcell);

$R ≔ polynomial_ring$

$F ≔ [x^{2} a + x b + c]$

$cad ≔ [cad_cell, cad_cell, cad_cell, cad_cell, cad_cell, cad_cell, cad_cell, cad_cell, cad_cell, cad_cell, cad_cell, cad_cell, cad_cell, cad_cell, cad_cell, cad_cell, cad_cell, cad_cell, cad_cell, cad_cell, cad_cell, cad_cell, cad_cell, cad_cell, cad_cell, cad_cell, cad_cell]$

(3.7.3)

>	nops(cad);

$27$

(3.7.4)

The output consists of $27$ cells, which can also been displayed as so-called Taski Formulas.

>	Display(cad, R);

$[\{\begin{array}{c} x = x \\ c < \frac{b^{2}}{4 a} \\ b = b \\ a < 0 \end{array}, \{\begin{array}{c} x < - \frac{b}{2 a} \\ c = \frac{b^{2}}{4 a} \\ b = b \\ a < 0 \end{array}, \{\begin{array}{c} x = - \frac{b}{2 a} \\ c = \frac{b^{2}}{4 a} \\ b = b \\ a < 0 \end{array}, \{\begin{array}{c} - \frac{b}{2 a} < x \\ c = \frac{b^{2}}{4 a} \\ b = b \\ a < 0 \end{array}, \{\begin{array}{c} x < \frac{- b + \sqrt{- 4 c a + b^{2}}}{2 a} \\ \frac{b^{2}}{4 a} < c \\ b = b \\ a < 0 \end{array}, \{\begin{array}{c} x = \frac{- b + \sqrt{- 4 c a + b^{2}}}{2 a} \\ \frac{b^{2}}{4 a} < c \\ b = b \\ a < 0 \end{array}, \{\begin{array}{c} \frac{- b + \sqrt{- 4 c a + b^{2}}}{2 a} < x \land x < - \frac{b + \sqrt{- 4 c a + b^{2}}}{2 a} \\ \frac{b^{2}}{4 a} < c \\ b = b \\ a < 0 \end{array}, \{\begin{array}{c} x = - \frac{b + \sqrt{- 4 c a + b^{2}}}{2 a} \\ \frac{b^{2}}{4 a} < c \\ b = b \\ a < 0 \end{array}, \{\begin{array}{c} - \frac{b + \sqrt{- 4 c a + b^{2}}}{2 a} < x \\ \frac{b^{2}}{4 a} < c \\ b = b \\ a < 0 \end{array}, \{\begin{array}{c} x < - \frac{c}{b} \\ c = c \\ b < 0 \\ a = 0 \end{array}, \{\begin{array}{c} x = - \frac{c}{b} \\ c = c \\ b < 0 \\ a = 0 \end{array}, \{\begin{array}{c} - \frac{c}{b} < x \\ c = c \\ b < 0 \\ a = 0 \end{array}, \{\begin{array}{c} x = x \\ c < 0 \\ b = 0 \\ a = 0 \end{array}, \{\begin{array}{c} x = x \\ c = 0 \\ b = 0 \\ a = 0 \end{array}, \{\begin{array}{c} x = x \\ 0 < c \\ b = 0 \\ a = 0 \end{array}, \{\begin{array}{c} x < - \frac{c}{b} \\ c = c \\ 0 < b \\ a = 0 \end{array}, \{\begin{array}{c} x = - \frac{c}{b} \\ c = c \\ 0 < b \\ a = 0 \end{array}, \{\begin{array}{c} - \frac{c}{b} < x \\ c = c \\ 0 < b \\ a = 0 \end{array}, \{\begin{array}{c} x < - \frac{b + \sqrt{- 4 c a + b^{2}}}{2 a} \\ c < \frac{b^{2}}{4 a} \\ b = b \\ 0 < a \end{array}, \{\begin{array}{c} x = - \frac{b + \sqrt{- 4 c a + b^{2}}}{2 a} \\ c < \frac{b^{2}}{4 a} \\ b = b \\ 0 < a \end{array}, \{\begin{array}{c} - \frac{b + \sqrt{- 4 c a + b^{2}}}{2 a} < x \land x < \frac{- b + \sqrt{- 4 c a + b^{2}}}{2 a} \\ c < \frac{b^{2}}{4 a} \\ b = b \\ 0 < a \end{array}, \{\begin{array}{c} x = \frac{- b + \sqrt{- 4 c a + b^{2}}}{2 a} \\ c < \frac{b^{2}}{4 a} \\ b = b \\ 0 < a \end{array}, \{\begin{array}{c} \frac{- b + \sqrt{- 4 c a + b^{2}}}{2 a} < x \\ c < \frac{b^{2}}{4 a} \\ b = b \\ 0 < a \end{array}, \{\begin{array}{c} x < - \frac{b}{2 a} \\ c = \frac{b^{2}}{4 a} \\ b = b \\ 0 < a \end{array}, \{\begin{array}{c} x = - \frac{b}{2 a} \\ c = \frac{b^{2}}{4 a} \\ b = b \\ 0 < a \end{array}, \{\begin{array}{c} - \frac{b}{2 a} < x \\ c = \frac{b^{2}}{4 a} \\ b = b \\ 0 < a \end{array}, \{\begin{array}{c} x = x \\ \frac{b^{2}}{4 a} < c \\ b = b \\ 0 < a \end{array}]$

(3.7.5)

Another important feature of the CAD command is the fact that it can take a list of semi-algebraic systems as input. This list of semi-algebraic systems represents a DNF quantifier free formula, that is, a disjunction of semi-algebraic systems.

Consider the union of a circle and a line.

>	Circle := x^2+y^2-1=0; Line := x-2=0; with(plots): implicitplot([Circle, Line], x = -1 .. 3, y = -1 .. 2, color=[red, blue]);

$Circle ≔ x^{2} + y^{2} - 1 = 0$

$Line ≔ x - 2 = 0$

Compute a CAD of the union of $Circle$ and $Line &period;$

>	R := PolynomialRing([y, x]); cad := CylindricalAlgebraicDecompose([[Circle], [Line]], R, output=cadcell);

$R ≔ polynomial_ring$

$cad ≔ [cad_cell, cad_cell, cad_cell, cad_cell, cad_cell]$

(3.7.6)

>	Display(cad, R);

$[\{\begin{array}{c} y = 0 \\ x = −1 \end{array}, \{\begin{array}{c} y = - \sqrt{- x^{2} + 1} \\ −1 < x \land x < 1 \end{array}, \{\begin{array}{c} y = \sqrt{- x^{2} + 1} \\ −1 < x \land x < 1 \end{array}, \{\begin{array}{c} y = 0 \\ x = 1 \end{array}, \{\begin{array}{c} y = y \\ x = 2 \end{array}]$

(3.7.7)

Each cad_cell has a sample point associated with it.

>	sp := map(SamplePoints, cad, R): Display(sp,R);

$[\{\begin{array}{c} y = 0 \\ x = −1 \end{array}, \{\begin{array}{c} y = −1 \\ x = 0 \end{array}, \{\begin{array}{c} y = 1 \\ x = 0 \end{array}, \{\begin{array}{c} y = 0 \\ x = 1 \end{array}, \{\begin{array}{c} y = 0 \\ x = 2 \end{array}]$

(3.7.8)

Last but not least, the CAD command is back engine of SolveTools:-SemiAlgebraic.

>	SolveTools:-SemiAlgebraic([x^2+x-c>0], [c, x]);

$[[c < - \frac{1}{4}, x = x], [c = - \frac{1}{4}, x \neq - \frac{1}{2}], [- \frac{1}{4} < c, x < - \frac{1}{2} - \frac{\sqrt{4 c + 1}}{2}], [- \frac{1}{4} < c, - \frac{1}{2} + \frac{\sqrt{4 c + 1}}{2} < x]]$

(3.7.9)

>	cad := CylindricalAlgebraicDecompose([[x^2+x-c>0]], PolynomialRing([x, c]), output=rootof);

$cad ≔ [[c < - \frac{1}{4}, x = x], [c = - \frac{1}{4}, x \neq - \frac{1}{2}], [- \frac{1}{4} < c, x < - \frac{1}{2} - \frac{\sqrt{4 c + 1}}{2}], [- \frac{1}{4} < c, - \frac{1}{2} + \frac{\sqrt{4 c + 1}}{2} < x]]$

(3.7.10)

3.8 Solving a linear semi-algebraic system

The input linear semi-algebraic system:

>	S := [f00 = c, f10 = c + cx1, f01 = c + cx2, f11 = c + cx1 + cx2 + cx1x2, cx1x2 <= 0];

$S ≔ [f00 = c, f10 = c + cx1, f01 = c + cx2, f11 = c + cx1 + cx2 + cx1x2, cx1x2 \leq 0]$

(3.8.1)

Define the elimination order (descending order).

>	R := PolynomialRing([cx1x2, cx1, cx2, c, f00, f01, f10, f11]):

Call LinearSolve to eliminate the variables.

The output is a set of equivalent linear equations and inequalities sorted in ascending order according to the largest variables appearing in the constraints. It provides conditions on lower order variables such that higher order variables having solutions.

In other words, the projection of the solutions of input system $S$ onto any lower dimensional space, say the space formed by the smallest $i$ variables, are exactly the solutions of those constraints in the output which only involve the smallest $i$ smallest variables.

>	R := PolynomialRing([cx1x2, cx1, cx2, c, f00, f01, f10, f11]): result := LinearSolve(S, R);

$result ≔ [f00 \leq - f11 + f10 + f01, c = f00, cx2 = f01 - f00, cx1 = f10 - f00, cx1x2 = f11 - f10 - f01 + f00]$

(3.8.2)

3.9 Verifying the output of different real solvers

On a given input polynomial system, two solving tools may produce correct results that look fairly different. Proving that these two results, say S1 and S2, are equivalent can be a very complex task. Here’s an example.

Given a triangle with edges $a, b, c$ (denoting the respective lengths by $a, b, c$ as well) the following two conditions $C1, C2$ are both characterizing the fact that the external bisector of the angle of $a, c$ intersects with $b$ on the other side of $a$ than the triangle:

We set S1 and S2 up first. S1 is the disjunction of C1 and C2.

C1:=[a > 0 , b > 0 , c > 0 , a < b + c ,
b < a + c , c < a + b, b^2 + a^2 - c^2 <= 0 ];

C2:=[a > 0 , b > 0 , c > 0 , a < b + c ,
b < a + c , c < a + b,
c*(b^2 + a^2 - c^2)^2 < a*b^2*(2*a*c - (c^2 + a^2 - b^2))]:

S1:=[C1, C2];

S2 := [a-c<0, a > 0 , b > 0 ,
c > 0 , a < b + c , b < a + c , c < a + b];

$C1 ≔ [0 < a, 0 < b, 0 < c, a < b + c, b < a + c, c < a + b, a^{2} + b^{2} - c^{2} \leq 0]$

$S1 ≔ [[0 < a, 0 < b, 0 < c, a < b + c, b < a + c, c < a + b, a^{2} + b^{2} - c^{2} \leq 0], [0 < a, 0 < b, 0 < c, a < b + c, b < a + c, c < a + b, c {(a^{2} + b^{2} - c^{2})}^{2} < a b^{2} (- a^{2} + 2 c a + b^{2} - c^{2})]]$

$S2 ≔ [a - c < 0, 0 < a, 0 < b, 0 < c, a < b + c, b < a + c, c < a + b]$

(3.9.1)

Compute regular semi-algebraic system representations dec1 (resp. dec2) for S1 (resp. S2).

>	R := PolynomialRing([a,b,c]): dec1 := map(op, map(RealTriangularize, S1,R)); dec2:= RealTriangularize(S2, R);

$dec1 ≔ [regular_semi_algebraic_system, regular_semi_algebraic_system, regular_semi_algebraic_system]$

$dec2 ≔ [regular_semi_algebraic_system]$

(3.9.2)

Compute the differences: S1 \ S2 and S2 \ S1.

>	Difference(dec1,dec2,R);

$[]$

(3.9.3)

>	Difference(dec2,dec1,R);

$[]$

(3.9.4)

3.10 Computing the projection of a semi-algebraic set

In the following problem, we are interesting in determining sufficient and necessary conditions on the variables a1, a2, for the prescribed semi-algebraic system to have solutions in r1, r2, x1, x2, e1, e2. This question is equivalent to compute the standard projection of the corresponding semi-algebraic set on the (a1, a2)-plane.

Problem: ∃(r1, r2, x1, x2, e1, e2) ( $\land_{f \in eqs}$ f=0 $\land_{p \in pie}$ p>0 )

vars := [r1, r2, x1, x2, e1, e2, a1, a2];
eqs := [
  r1^2 - r2^2 - e1^2 + e2^2 + x1^2 - 2*x1*x2 + x2^2,
  2*r1*r2 - 2*e1*e2,
  -r1 + e1 + x1*a2 - x2*a2,
  -r2 + e2 - x1*a1 + x2*a1,
  r1^2 + r2^2 + x1^2 + x2^2 + e1^2 + e2^2 - 1];

pie := [r1, -e1];
ineqs := [];
nie := [];
params := [a1,a2];

R := PolynomialRing(vars):
proj := Projection(eqs, nie, pie, ineqs, 2, R);

$vars ≔ [r1, r2, x1, x2, e1, e2, a1, a2]$

$eqs ≔ [- {e1}^{2} + {e2}^{2} + {r1}^{2} - {r2}^{2} + {x1}^{2} - 2 x1 x2 + {x2}^{2}, - 2 e1 e2 + 2 r1 r2, x1 a2 - x2 a2 + e1 - r1, - x1 a1 + x2 a1 + e2 - r2, {e1}^{2} + {e2}^{2} + {r1}^{2} + {r2}^{2} + {x1}^{2} + {x2}^{2} - 1]$

$pie ≔ [r1, - e1]$

$ineqs ≔ []$

$nie ≔ []$

$params ≔ [a1, a2]$

$proj ≔ [regular_semi_algebraic_system, regular_semi_algebraic_system]$

(3.10.1)

Show the result:

>	Display(proj, R);

$[\{\begin{array}{c} a1 = 0 \\ \{\begin{array}{c} a2 < −1 \\ or a2 - 1 > 0 \end{array} \end{array}, {a1}^{2} + {a2}^{2} - 1 > 0 and a1 \neq 0 and a2 \neq 0]$

(3.10.2)

4. Linear algebra over towers of (field) extensions

Regular chains encode towers of transcendental and algebraic extensions of the base field. In practice, these towers themselves are not always fields and may have zero-divisors. Many standard algorithms (for example, for solving systems of linear equations) require that the coefficient ring be a field. This difficulty is overcome by means of the so-called D5 Principle. This principle is at the core of the theory of regular chains and allows us to generalize many algorithms which are a priori restricted to coefficients from a field.

The MatrixTools submodule applies the D5 Principle to algorithms for solving linear systems with coefficient rings encoded by regular chains. Below, use of this submodule is demonstrated by two examples. As discussed earlier regarding polynomial gcd computations, splitting naturally occurs when computing modulo regular chains. Sometimes, split results can be recombined as illustrated below.

4.1 Automatic case discussion and recombination (I)

Assume that there are two algebraic entities, $y$ and $z$ ; that they have the same square; and that $z$ is a4th root of $- 1$ . Suppose you need to perform algebraic computations with $y$ and $z$ .

>	R := PolynomialRing([y, z]): rc := Chain([z^4+1, y^2-z^2], Empty(R), R): Equations(rc, R);

$[y^{2} - z^{2}, z^{4} + 1]$

(4.1.1)

For example, you may want to compute the inverse of the following matrix.

>	$A ≔ [\begin{array}{c} 1 & y + z \\ 0 & y - z \end{array}] &colon;$

Clearly, the result depends on whether $y$ and $z$ are equal or not. Note that from the assumptions above, neither $y = z$ nor $y = - z$ can be deduced. When $y \neq z$ , the matrix $A$ has an inverse. When they are equal, the matrix $A$ is singular. These facts are detected automatically by the command MatrixInverse.

>	result := MatrixInverse(A, rc, R);

$result ≔ [[[[\begin{array}{c} 1 & 0 \\ 0 & \frac{z^{3}}{2} \end{array}], regular_chain]], [[noInv, [\begin{array}{c} 1 & y + z \\ 0 & y - z \end{array}], regular_chain]]]$

(4.1.2)

Check the first result. Note the use of the command MatrixMultiply in order to multiply two matrices modulo a regular chain.

>	$B, rc1 ≔ op (result [1] [1]) &semi;$ $Equations (rc1, R) &semi;$ $MatrixMultiply (B, A, rc1, R)$

$B, rc1 ≔ [\begin{array}{c} 1 & 0 \\ 0 & \frac{z^{3}}{2} \end{array}], regular_chain$

$[y + z, z^{4} + 1]$

$[\begin{array}{c} 1 & 0 \\ 0 & 1 \end{array}]$

(4.1.3)

>	${rc2}_{} ≔ result [2] [1] [3] &semi;$ $Equations (rc2, R) &semi;$

$rc2 ≔ regular_chain$

$[y - z, z^{4} + 1]$

(4.1.4)

You can see that the computation (and the regular chain) was split into two branches. The first one corresponds to $y + z = 0$ , and in that branch the matrix is invertible. The second branch corresponds to $y = z$ , and in that branch the matrix is not invertible.

Consider the matrix $M$ below and compute its inverse modulo the regular chain $rc$ . Here again, computations split. However, the matrix $M$ is invertible in both cases.

>	M := Matrix([[1, y+z], [2, y-z]]); result := MatrixInverse(M, rc, R);

$M ≔ [\begin{array}{c} 1 & y + z \\ 2 & y - z \end{array}]$

$result ≔ [[[[\begin{array}{c} 1 & 0 \\ - z^{3} & \frac{z^{3}}{2} \end{array}], regular_chain], [[\begin{array}{c} 0 & \frac{1}{2} \\ - \frac{z^{3}}{2} & \frac{z^{3}}{4} \end{array}], regular_chain]], []]$

(4.1.5)

Double check this result.

>	op(result[1][1]);

$[\begin{array}{c} 1 & 0 \\ - z^{3} & \frac{z^{3}}{2} \end{array}], regular_chain$

(4.1.6)

>	$unassign (' N') &semi; N_{1}, rc1 ≔ op ({result}_{} [1] [1]) &semi;$ $Equations (rc1, R) &semi;$ $N_{1} &semi; MatrixMultiply (N_{1}, M, rc1, R)$

$N_{1}, rc1 ≔ [\begin{array}{c} 1 & 0 \\ - z^{3} & \frac{z^{3}}{2} \end{array}], regular_chain$

$[y + z, z^{4} + 1]$

$[\begin{array}{c} 1 & 0 \\ - z^{3} & \frac{z^{3}}{2} \end{array}]$

$[\begin{array}{c} 1 & 0 \\ 0 & 1 \end{array}]$

(4.1.7)

>	$N_{2}, rc2 ≔ op ({result}_{} [1] [2]) &semi;$ $Equations (rc2, R) &semi;$ $MatrixMultiply (N_{2}, M, rc2, R)$

$N_{2}, rc2 ≔ [\begin{array}{c} 0 & \frac{1}{2} \\ - \frac{z^{3}}{2} & \frac{z^{3}}{4} \end{array}], regular_chain$

$[y - z, z^{4} + 1]$

$[\begin{array}{c} 1 & 0 \\ 0 & 1 \end{array}]$

(4.1.8)

Since the matrix $M$ is invertible in both branches given by the regular chains ${rc}_{1}$ and ${rc}_{2}$ , it is natural to ask whether there is a "generic" answer that would hold for both cases. The answer is $yes$ . Technically speaking, one can observe that the Chinese Remainder Theorem applies since the (saturated) ideals of ${rc}_{1}$ and ${rc}_{2}$ are relatively prime. The command MatrixCombine implements this recombination.

>	$combined := MatrixCombine ([rc1, rc2], R, [N_{1}, N_{2}])$

$combined ≔ [[[\begin{array}{c} \frac{y z^{3}}{2} + \frac{1}{2} & - \frac{y z^{3}}{4} + \frac{1}{4} \\ \frac{1}{4} y z^{2} - \frac{3}{4} z^{3} & - \frac{1}{8} y z^{2} + \frac{3}{8} z^{3} \end{array}], regular_chain]]$

(4.1.9)

Check that the above matrix times $M$ gives the identity matrix modulo $rc$ .

>	MatrixMultiply(combined[1][1], M, combined[1][2], R);

$[\begin{array}{c} 1 & 0 \\ 0 & 1 \end{array}]$

(4.1.10)

4.2 Automatic case discussion and recombination (II)

Can there be several cases in the output of the command MatrixCombine? In other words, is it possible that this command fails recombining several cases into one, even if the Chinese Remainder Theorem applies? Yes, this can happen, for algebraic reasons that shall be explained below. For starters, consider one of the previous polynomial systems.

>	sys := {x+y+z^2-1, x+y^2+z-1, x^2+y+z-1}: R := PolynomialRing([x, y, z]): l := Triangularize(sys, R, normalized = yes); map(Equations, l, R);

$l ≔ [regular_chain, regular_chain, regular_chain, regular_chain]$

$[[x - z, y - z, z^{2} + 2 z - 1], [x, y, z - 1], [x, y - 1, z], [x - 1, y, z]]$

(4.2.1)

Next, generate four random $2 \times 2$ matrices with polynomial entries.

>	A := [seq(Matrix([seq([seq(randpoly([x, y, z], degree = 1), j = 1 .. 2)], i = 1 .. 2)]), k = 1 .. 4)];

$A ≔ [[\begin{array}{c} - 20 x + 51 y - 87 z - 37 & 2 x + 94 y + 70 z + 73 \\ 85 x + 21 y + 51 z - 6 & 10 x - 16 y - 95 z - 16 \end{array}], [\begin{array}{c} - 69 x - 53 y - 53 z + 65 & - 45 x - 54 y + 11 z + 38 \\ - 29 x - 27 y - 95 z - 39 & 93 x + 55 y + 65 z + 23 \end{array}], [\begin{array}{c} - 71 x + 86 y + 77 z + 45 & 36 x - 15 y + 64 z + 76 \\ 99 x + 80 y + 77 z + 36 & 40 x + 6 y + 92 z + 38 \end{array}], [\begin{array}{c} 7 x + 97 y + 68 z - 58 & - 54 x - 6 y + 63 z + 47 \\ 21 x + 4 y + 43 z - 16 & - 81 x + 93 y + 93 z - 84 \end{array}]]$

(4.2.2)

Attempt the recombination of the four cases given by the regular chains in $A$ , and observe that MatrixCombine produces two cases.

>	combined := MatrixCombine(l, R, A);

$combined ≔ [[[\begin{array}{c} 182 y - 51 & 68 y - 7 \\ 111 y + 5 & 209 y - 165 \end{array}], regular_chain], [[\begin{array}{c} \frac{105}{2} z^{2} + 49 z - \frac{179}{2} & - 95 z^{2} - 24 z + 168 \\ - \frac{285}{2} z^{2} - 128 z + \frac{273}{2} & \frac{205}{2} z^{2} + 104 z - \frac{237}{2} \end{array}], regular_chain]]$

(4.2.3)

Now investigate why these two cases cannot be merged into a single one.

>	$rc1 ≔ combined [1] [2] &colon; Equations (rc1, R) &semi;$ $rc2 ≔ {combined}_{} [2] [2] &colon;$ $Equations (rc2, R)$

$[x + y - 1, y^{2} - y, z]$

$[2 x + z^{2} - 1, 2 y + z^{2} - 1, z^{3} + z^{2} - 3 z + 1]$

(4.2.4)

The two ideals generated by ${rc}_{1}$ and ${rc}_{2}$ are obviously relatively prime (no common roots in $z$ ). Thus, the Chinese Remainder Theorem applies. But, if you try to recombine ${rc}_{1}$ and ${rc}_{2}$ into a single system, this will create a polynomial in $y$ with a zero-divisor as a leading coefficient, as seen below. This is forbidden by the properties of a regular chain.

First, create two new regular chains with the polynomials from ${rc}_{1}$ and ${rc}_{2}$ that are univariate in the variable $z$ .

>	$S ≔ PolynomialRing ([z]) &colon;$ $rc3 ≔ Chain ([z], Empty (S), (S)) &colon;$ $rc4 ≔ Chain ([z^{3} + z^{2} - 3 z + 1], Empty (S), S) &colon;$

Then define two $1 \times 1$ matrices $U_{1}, U_{2}$ , and combine them with regards to ${rc}_{3}$ and ${rc}_{4} &period;$

>	$U_{1} ≔ Matrix ([[y^{2} - y]]) &semi;$ $U_{2} ≔ Matrix ([[2 y + z^{2} - 1]]) &semi;$ $combined ≔ MatrixCombine ([rc3, rc4], S, [U_{1}, U_{2}]) &semi;$ $V, rc5 ≔ op ({combined}_{1}) &colon; Equations (rc5, S)$

$U_{1} ≔ [\begin{array}{c} y^{2} - y \end{array}]$

$U_{2} ≔ [\begin{array}{c} 2 y + z^{2} - 1 \end{array}]$

$combined ≔ [[[\begin{array}{c} y^{2} z^{3} + y^{2} z^{2} - 3 y z^{3} - 3 y^{2} z - 3 y z^{2} + z^{3} + y^{2} + 9 y z + 2 z^{2} - y - 3 z \end{array}], regular_chain]]$

$[z^{4} + z^{3} - 3 z^{2} + z]$

(4.2.5)

So, the recombination is successful. However, the initial $p$ (with regards to $z$ ) of the resulting polynomial is a zero divisor: it vanishes modulo ${rc}_{4}$ and is invertible modulo ${rc}_{3}$ . Consequently, the computation of the inverse using the command Inverse splits the combined regular chain ${rc}_{5}$ .

>	$p ≔ Initial (V_{1, 1}, R) &semi;$ $inverse ≔ Inverse (p, rc5, S) &semi;$ $Equations ({inverse}_{} [1] [1] [3], S) &semi;$ $map (Equations, {inverse}_{2}, S)$

$p ≔ z^{3} + z^{2} - 3 z + 1$

$inverse ≔ [[[1, 1, regular_chain]], [regular_chain, regular_chain]]$

$[z]$

$[[z - 1], [z^{2} + 2 z - 1]]$

(4.2.6)

5. Constructible sets and rational maps

Constructible sets are the geometrical objects naturally attached to triangular decompositions as polynomial ideals are the algebraic concept underlying the computation of Gröbner bases. This relation becomes even more complex and essential in the case of polynomial systems with infinitely many solutions. Basically, constructible sets are what you get when you take algebraic varieties (defined by polynomial ideals) and add the operation of taking complements. More precisely, a constructible set is either a set of points defined by both equations and inequalities, or a finite union of such sets. Note that, in general, the complement of an algebraic variety cannot be described by a polynomial ideal, but it is a constructible set. Constructible sets have the nice property that they are closed under intersections and finite unions, like polynomial ideals, and additionally under complements.

This section presents the ConstructibleSetTools submodule of the RegularChains library. To our knowledge, this is the first general-purpose computer algebra package providing constructible set as a type and exporting a rich collection of operations for manipulating constructible sets. Besides, this module provides routines in support of solving parametric polynomial systems, and several of its commands will be demonstrated in other parts of this document. The examples of the present section illustrate the Theorem of Chevalley which states that the image of a constructible set by a rational map is a constructible set.

5.1 Some high-school examples

Constructible sets appear naturally in many elementary questions from high-school problems. A simple one is the following: for which values of x does $f (x, y) = 0$ have solutions, with $f$ as below?

>	$R ≔ PolynomialRing ([y, x]) &colon;$ $f ≔ 1 - x \cdot y &colon;$

The answer is whenever $x \neq 0$ holds. Formally speaking, what we want there is the projection onto the x-axis of the hyperbola $f = 0$ . This object is not the zero set of a system of polynomial equations, and therefore it cannot be computed directly via traditional techniques such as a Gröbner basis computation. It requires some finer and more geometrical elimination process. The function Projection implements such a process via triangular decompositions. The output is an object of type constructible_set. Its internal representation is given by a list of so called regular systems.

>	$cs := Projection (\{f\}, 1, R) &semi;$ $l ≔ RepresentingRegularSystems (cs, R) &semi;$ $rs ≔ l_{1} &colon;$

$cs ≔ constructible_set$

$l ≔ [regular_system]$

(5.1.1)

Each regular system is a pair consisting of a regular chain plus one or more inequations:

>	$rc ≔ ConstructibleSetTools [RepresentingChain] (rs, R) &semi;$ $IsEmptyChain (rc, R) &semi;$ $map (`≠`, RepresentingInequations (rs, R), 0)$

$rc ≔ regular_chain$

$true$

$[x \neq 0]$

(5.1.2)

In this example, this regular chain $rc$ is just the empty one and the inequation is just $x \neq 0$ . This output may look more complicated than the posed problem itself. So, consider another familiar example but less trivial example: for which values of a, b, c does the equation $a x^{2} + b x + c = 0$ have solutions? This can be seen as another "projection" question.

>	$R := PolynomialRing ([x, a, b, c]) &semi;$ $f := a x^{2} + b x + c &semi;$ $cs ≔ Projection (\{f\}, 3, R) &semi;$ $RepresentingRegularSystems (cs, R) &semi;$ $Display (cs, R) &semi;$

$R ≔ polynomial_ring$

$f ≔ x^{2} a + x b + c$

$cs ≔ constructible_set$

$[regular_system, regular_system, regular_system]$

$a \neq 0, \{\begin{array}{c} a = 0 \\ b \neq 0 \end{array}, \{\begin{array}{c} a = 0 \\ b = 0 \\ c = 0 \end{array}$

(5.1.3)

We obtain three regular systems:
    (1) $a \neq 0$ ,
    (2) $a = 0,$ $b \neq 0$ ,
    (3) $a = b = c = 0$ .

5.2 Polynomial map images

The following picture illustrates a polynomial map $Π &colon; {&complexes;}^{2} \to {&complexes;}^{3}$ . The task is to describe the image of this map in ${&complexes;}^{3}$ . As one can see from the picture, a point with coordinates $(x, y, z) = (0, b, 0)$ , where $b \neq 0$ , cannot be in the image of $Π$ , whereas any other point with $z = 0$ is. For this reason, the image of $Π$ is a constructible set and not the zero set of a system of polynomial equations. The command PolynomialMapImage below computes the image of $Π$ .

>	$unassign (' u',' v') &semi; S ≔ PolynomialRing ([u, v]) &colon;$ $T ≔ PolynomialRing ([x, y, z]) &colon;$

After specifying the coordinates of the source space $S$ and target space $T$ , the polynomial map is defined and its image is computed.

>	$Π ≔ [u, u v, 0] &semi;$ $cs ≔ PolynomialMapImage ([], Π, S, T) &semi;$

$Π ≔ [u, u v, 0]$

$cs ≔ constructible_set$

(5.2.1)

The constructible set $cs$ is given by a list of two regular systems. The first one is given by $z = 0 \neq x$ , and the second one by $x = y = z = 0$ .

>	$lrs := RepresentingRegularSystems (cs, T) &semi;$ $map (Display, lrs, T)$

$lrs ≔ [regular_system, regular_system]$

$[\{\begin{array}{c} z = 0 \\ x \neq 0 \end{array}, \{\begin{array}{c} x = 0 \\ y = 0 \\ z = 0 \end{array}]$

(5.2.2)

5.3 Rational map images

Images (and pre-images) of constructible sets under rational maps are also supported by the ConstructibleSetTools submodule. To demonstrate this facility, consider the implicit equation of a curve (namely the tacnode curve) given by a parametric representation involving rational functions. The parametrization of this curve can be seen as a rational map; its one-dimensional source space $S$ , its two-dimensional target space $T$ , and the image $ρ$ of an arbitrary point are defined below.

>	$S := PolynomialRing ([t]) &colon;$ $T := PolynomialRing ([x, y]) &colon;$ $ρ ≔ [\frac{t^{3} - 6 t^{2} + 9 t - 2}{2 t^{4} - 16 t^{3} + 40 t^{2} - 32 t + 9}, \frac{t^{2} - 4 t + 4}{2 t^{4} - 16 t^{3} + 40 t^{2} - 32 t + 9}] &colon;$

Here is a parametric plot of the tacnode curve.

>	$plot ([op (ρ), t = - 10 .. 10]) &semi;$

The parametric explicit representation of a curve is very useful for plotting it. However, in order to answer questions such as "does a given point in target space lie on the curve?", an implicit representation is more useful.
The image of the full one-dimensional space under the rational map $ρ$ is computed below.

>	$F := [] &colon;$ ${cs}_{1} ≔ RationalMapImage (F, ρ, S, T) &semi;$ $l ≔ RepresentingRegularSystems ({cs}_{1}, T) &semi; Display ({cs}_{1}, T) &semi;$

${cs}_{1} ≔ constructible_set$

$l ≔ [regular_system, regular_system, regular_system]$

$\{\begin{array}{c} 2 x^{4} - 3 y x^{2} + y^{4} - 2 y^{3} + y^{2} = 0 \\ 10 y x^{2} + 2 y^{3} + 2 x^{2} - y^{2} - y \neq 0 \\ 964 x^{2} y^{6} - 88 y^{8} - 480 x^{2} y^{5} + 2104 y^{7} - 6858 x^{2} y^{4} - 2316 y^{6} - 4328 y^{3} x^{2} - 943 y^{5} - 888 y^{2} x^{2} + 892 y^{4} - 72 y x^{2} + 318 y^{3} - 2 x^{2} + 32 y^{2} + y \neq 0 \\ y \neq 0 \end{array}, \{\begin{array}{c} x = 0 \\ y - 1 = 0 \end{array}, \{\begin{array}{c} x = 0 \\ y = 0 \end{array}$

(5.3.1)

The result has three components. The second and third ones correspond to the self-intersection points of the tacnode curve, that is, its singularities. The first component $l_{1}$ defines all the other points. The role of its inequations is to exclude the two self-intersection points, so that the sets described by the three components are disjoint.

>	$rc ≔ ConstructibleSetTools [RepresentingChain] (l_{1}, T) &colon;$ $eq ≔ Equations (rc, T) &semi;$ $ineq ≔ ConstructibleSetTools [RepresentingInequations] (l_{1}, T)$

$eq ≔ [2 x^{4} - 3 y x^{2} + y^{4} - 2 y^{3} + y^{2}]$

$ineq ≔ [y, 10 y x^{2} + 2 y^{3} + 2 x^{2} - y^{2} - y, 964 x^{2} y^{6} - 88 y^{8} - 480 x^{2} y^{5} + 2104 y^{7} - 6858 x^{2} y^{4} - 2316 y^{6} - 4328 y^{3} x^{2} - 943 y^{5} - 888 y^{2} x^{2} + 892 y^{4} - 72 y x^{2} + 318 y^{3} - 2 x^{2} + 32 y^{2} + y]$

(5.3.2)

The computations below show that the solution set of the equations of the first component $l_{1}$ correspond to the entire tacnode curve. Note the use of the command IsContained which decides whether a constructible set is contained in another.

>	${cs}_{2} ≔ GeneralConstruct (eq, [], T) &semi;$ $IsContained ({cs}_{1}, {cs}_{2}, T) and IsContained ({cs}_{2}, {cs}_{1}, T)$

${cs}_{2} ≔ constructible_set$

$true$

(5.3.3)

Therefore, $2 x^{4} - 3 y x^{2} + y^{2} - 2 y^{3} + y^{4} = 0$ is the implicit representation of the tacnode curve.

6. Parametric polynomial systems

The ParametricSystemTools module is devoted to solving systems with parameters, including real root classification and complex root classification of such systems, as demonstrated below. The first two examples are dedicated to the complex solutions of parametric polynomial systems. The last, but very detailed, third subsection deals with the real case.

6.1 An example from classical invariant theory

Each of the following polynomials defines an elliptic curve in the complex plane. They depend on parameters $a_{1}$ and $a_{2}$ , respectively. In invariant theory, a classical question is whether there exists a linear fractional map from the first curve to the second one.

>	$f_{1} ≔ y^{2} - x^{3} - a_{1} x - 1 &semi;$ $f_{2} ≔ y^{2} - x^{3} - a_{2} x - 1$

$f_{1} ≔ - x^{3} - a_{1} x + y^{2} - 1$

$f_{2} ≔ - x^{3} - a_{2} x + y^{2} - 1$

(6.1.1)

These two curves are plotted below for $a_{1} = 1$ and $a_{2} = 4$ .

>	$with (plots) &colon;$ $display (implicitplot (\binom{f_{1}}{} \| \binom{}{a_{1} = 1}, x = - 3 .. 3, y = - 3 .. 3), implicitplot (\binom{f_{2}}{} \| \binom{}{a_{2} = 4}, x = - 3 .. 3, y = - 3 .. 3, color = blue))$

A generic linear fractional map is applied to the coordinates $x$ and $y$ of the second curve.

>	$unassign (' A',' B',' C',' E',' F',' H') &colon;$ $f_{3} ≔ eval (f_{2}, [x = \frac{A x + B y + C}{G x + H y + K}, y = \frac{D x + E y + F}{G x + H y + K}])$

$f_{3} ≔ - \frac{{(A x + B y + C)}^{3}}{{(G x + H y + K)}^{3}} - \frac{(A x + B y + C) a_{2}}{G x + H y + K} + \frac{{(D x + E y + F)}^{2}}{{(G x + H y + K)}^{2}} - 1$

(6.1.2)

Next, one stipulates that the rational function $f_{1} - f_{3}$ must be identically zero. This yields a system of equations $sys$ . Without loss of generality, one can assume that the origin is mapped to the origin, which implies $C = F = 0$ .

>	$sys := [coeffs (numer (f_{1} - f_{3}), [x, y]), C, F] &semi;$ $R := PolynomialRing ([B, E, H, A, D, G, K, F, C, a_{1}, a_{2}]) &colon;$

$sys ≔ [- G^{3}, - 3 G^{2} H, - 3 G^{2} K, - 3 G H^{2}, - 6 G H K, - G^{3} a_{1} - 3 G K^{2}, - H^{3}, G^{3} - 3 H^{2} K, - 3 G^{2} H a_{1} - 3 H K^{2}, A G^{2} a_{2} - 3 G^{2} K a_{1} + A^{3} - D^{2} G - K^{3}, 3 G^{2} H, - 3 G H^{2} a_{1} + 3 G^{2} K, 2 A G H a_{2} + B G^{2} a_{2} - 6 G H K a_{1} + 3 A^{2} B - D^{2} H - 2 D E G, 2 A G K a_{2} + C G^{2} a_{2} - 3 G K^{2} a_{1} + 3 A^{2} C - D^{2} K - 2 D F G, 3 G H^{2}, - H^{3} a_{1} + 6 G H K, A H^{2} a_{2} + 2 B G H a_{2} - 3 H^{2} K a_{1} + 3 A B^{2} - 2 D E H - E^{2} G + 3 G K^{2}, 2 A H K a_{2} + 2 B G K a_{2} + 2 C G H a_{2} - 3 H K^{2} a_{1} + 6 A B C - 2 D E K - 2 D F H - 2 E F G, A K^{2} a_{2} + 2 C G K a_{2} - K^{3} a_{1} + 3 A C^{2} - 2 D F K - F^{2} G, H^{3}, 3 H^{2} K, B H^{2} a_{2} + B^{3} - E^{2} H + 3 H K^{2}, 2 B H K a_{2} + C H^{2} a_{2} + 3 B^{2} C - E^{2} K - 2 E F H + K^{3}, B K^{2} a_{2} + 2 C H K a_{2} + 3 B C^{2} - 2 E F K - F^{2} H, C K^{2} a_{2} + C^{3} - F^{2} K, C, F]$

(6.1.3)

To the system of equations $sys$ , inequation constraints need to be added. Indeed, the unknowns $G, H$ , and $K$ cannot vanish simultaneously. Hence, the system to be solved consists of the zeros of $sys$ , which are not common zeros of these three polynomials. This new system is easily expressed as the set-theoretical difference of two constructible sets.

>	${cs}_{1} ≔ GeneralConstruct (sys, [], R) &colon;$ ${cs}_{2} ≔ GeneralConstruct ([G, H, K], [], R) &colon;$ ${cs}_{3} ≔ Difference ({cs}_{1}, {cs}_{2}, R) &colon;$ $Display ({cs}_{3}, R)$

$\{\begin{array}{c} B = 0 \\ E - K = 0 \\ H = 0 \\ a_{2} A - K a_{1} = 0 \\ D = 0 \\ G = 0 \\ F = 0 \\ C = 0 \\ a_{1}^{2} + a_{2} a_{1} + a_{2}^{2} = 0 \\ K \neq 0 \\ a_{2} \neq 0 \end{array}, \{\begin{array}{c} B = 0 \\ E + K = 0 \\ H = 0 \\ a_{2} A - K a_{1} = 0 \\ D = 0 \\ G = 0 \\ F = 0 \\ C = 0 \\ a_{1}^{2} + a_{2} a_{1} + a_{2}^{2} = 0 \\ K \neq 0 \\ a_{2} \neq 0 \end{array}, \{\begin{array}{c} B = 0 \\ E - K = 0 \\ H = 0 \\ A - K = 0 \\ D = 0 \\ G = 0 \\ F = 0 \\ C = 0 \\ a_{1} - a_{2} = 0 \\ K \neq 0 \end{array}, \{\begin{array}{c} B = 0 \\ E + K = 0 \\ H = 0 \\ A - K = 0 \\ D = 0 \\ G = 0 \\ F = 0 \\ C = 0 \\ a_{1} - a_{2} = 0 \\ K \neq 0 \end{array}, \{\begin{array}{c} B = 0 \\ E - K = 0 \\ H = 0 \\ A^{2} + K A + K^{2} = 0 \\ D = 0 \\ G = 0 \\ F = 0 \\ C = 0 \\ a_{1} = 0 \\ a_{2} = 0 \\ K \neq 0 \end{array}, \{\begin{array}{c} B = 0 \\ E + K = 0 \\ H = 0 \\ A^{2} + K A + K^{2} = 0 \\ D = 0 \\ G = 0 \\ F = 0 \\ C = 0 \\ a_{1} = 0 \\ a_{2} = 0 \\ K \neq 0 \end{array}$

(6.1.4)

The question can now be stated in algebraic terms: for which values of the parameters $a_{1}$ and $a_{2}$ is the constructible set non-empty? Formulas solving for the unknowns are also desired. The command ComprehensiveTriangularize addresses these two requirements. The second argument, $2$ in the example below, specifies that the last two indeterminates of the polynomial ring (that is, $a_{1}$ and $a_{2}$ ) are the parameters of the system.

>	$ct := ComprehensiveTriangularize ({cs}_{3}, 2, R)$

$ct ≔ [regular_system, regular_system, regular_system, regular_system, regular_system, regular_system], [[constructible_set, [3, 6, 4, 5]], [constructible_set, [3, 4]], [constructible_set, [1, 2]]]$

(6.1.5)

The second list returned by this command is a finite partition of the projection of $cs$ onto the parameter space. Each component $P$ of this partition is a constructible set above which the geometrical properties of $cs$ (degree, dimension, etc.) are essentially constant; the corresponding component of ${cs}_{3}$ is given by the regular systems of $ct$ whose indices are associated with $P$ . Below, the three parts of this partition are joined into a single constructible set.

>	d := map2(op, 1, ct[2]): d := Union(d[1], Union(d[2], d[3], R), R): Display(d, R);

$\{\begin{array}{c} a_{1} = 0 \\ a_{2} = 0 \end{array}, \{\begin{array}{c} a_{1} - a_{2} = 0 \\ a_{2} \neq 0 \end{array}, \{\begin{array}{c} a_{1}^{2} + a_{2} a_{1} + a_{2}^{2} = 0 \\ a_{2} \neq 0 \end{array}$

(6.1.6)

Next, check if the answer to the question is the one that is well-known from invariant theory, namely that the two curves can be matched provided that $a_{1}^{3} =$ $a_{2}^{3}$ .

>	$e := GeneralConstruct ([a_{1}^{3} - a_{2}^{3}], [], R) &semi;$ $`and` (IsContained (d, e, R), IsContained (e, d, R))$

$e ≔ constructible_set$

$true$

(6.1.7)

Finally, the three regular systems (from the first list returned by ComprehensiveTriangularize) are displayed below.

>	for i to 3 do Display(ct[1][i], R); end do;

$\{\begin{array}{c} B = 0 \\ E - K = 0 \\ H = 0 \\ a_{2} A - K a_{1} = 0 \\ D = 0 \\ G = 0 \\ F = 0 \\ C = 0 \\ a_{1}^{2} + a_{2} a_{1} + a_{2}^{2} = 0 \\ K \neq 0 \\ a_{2} \neq 0 \end{array}$

$\{\begin{array}{c} B = 0 \\ E + K = 0 \\ H = 0 \\ a_{2} A - K a_{1} = 0 \\ D = 0 \\ G = 0 \\ F = 0 \\ C = 0 \\ a_{1}^{2} + a_{2} a_{1} + a_{2}^{2} = 0 \\ K \neq 0 \\ a_{2} \neq 0 \end{array}$

$\{\begin{array}{c} B = 0 \\ E - K = 0 \\ H = 0 \\ A - K = 0 \\ D = 0 \\ G = 0 \\ F = 0 \\ C = 0 \\ a_{1} - a_{2} = 0 \\ K \neq 0 \end{array}$

(6.1.8)

6.2 Counting complex solutions

In this short example, it is shown how to determine the number of distinct roots of a generic polynomial equation of degree $4$ depending on the parameters $a, b, c, d$ .

>	unassign('a', 'b', 'c', 'd'): R := PolynomialRing([x, a, b, c, d]): p := ax^4+bx^2+c*x+d;

$p ≔ a x^{4} + b x^{2} + c x + d$

(6.2.1)

The command ComplexRootClassification is precisely adapted to this purpose. It determines the number of distinct complex solutions of a parametric system (possibly a parametric constructible set) depending on parameters. The second argument, $4$ in the example below, indicates that the last four variables of the polynomial ring, $a, b, c, d$ , are the parameters of the equation.

>	cr := ComplexRootClassification([p], 4, R);

$cr ≔ [[constructible_set, 1], [constructible_set, 2], [constructible_set, 3], [constructible_set, 4], [constructible_set, \infty]]$

(6.2.2)

The result is a partition of the parameter space into five sets: those parameter values for which $p$ has exactly $1, 2, 3$ , or $4$ complex roots, respectively, and those parameter values for which $p$ has infinitely many roots (which happens only if $a = b = c = d = 0$ ).

Suppose you are interested in the case of three distinct roots. Display the equations and inequations of the regular systems defining the constructible set, giving the necessary and sufficient conditions for the polynomial $p_{}$ to have three distinct complex roots.

>	Info(cr[3][1], R);

$[[256 d^{3} a^{2} + (- 128 d^{2} b^{2} + 144 c^{2} d b - 27 c^{4}) a + 16 b^{4} d - 4 b^{3} c^{2}], [d, a, (8 d b - 9 c^{2}) a - 2 b^{3}, c, 32 d b - 9 c^{2}, (1073741824 b^{10} d^{10} - 66437775360 b^{9} c^{2} d^{9} + 672682475520 b^{8} c^{4} d^{8} - 2623461654528 c^{6} d^{7} b^{7} + 5164940132352 c^{8} d^{6} b^{6} - 5810557648896 c^{10} d^{5} b^{5} + 3961743851520 c^{12} d^{4} b^{4} - 1665238487040 c^{14} d^{3} b^{3} + 421513492032 c^{16} d^{2} b^{2} - 58887914328 c^{18} d b + 3486784401 c^{20}) a - 268435456 b^{12} d^{9} + 12683575296 b^{11} c^{2} d^{8} - 98524200960 b^{10} c^{4} d^{7} + 294298583040 b^{9} c^{6} d^{6} - 439871275008 b^{8} c^{8} d^{5} + 368924295168 b^{7} c^{10} d^{4} - 181579926528 b^{6} c^{12} d^{3} + 52038702720 b^{5} c^{14} d^{2} - 8035387920 b^{4} c^{16} d + 516560652 b^{3} c^{18}]], [[27 c^{2} a + 4 b^{3}, d], [c, b]], [[c, d], [b, a]]$

(6.2.3)

This output is to be interpreted as follows. There are three cases:

•	$c = d = 0$ and $a \neq 0 \neq b$ . In this case, the polynomial degenerates to $p = a x^{4} + {bx}^{2}$ , and the three roots are a double root at $x = 0$ and two distinct simple roots at $\pm \sqrt{- \frac{b}{a}}$ .

•

$4 b^{3} + 27 a c^{2} = 0 = d$ and $b \neq 0 \neq c$ . In this case, the polynomial degenerates to $p = a x^{4} + b x^{2} + c x$ , which has $x = 0$ as a simple root. The expression $4 b^{3} + 27 a c^{2}$ is the discriminant of the cofactor $\frac{p}{x} = a x^{3} + b x + c$ , and the condition that it must be zero implies that this polynomial has at least a double root. Finally, the condition $b \neq 0$ ensures that the latter polynomial does not have a triple root, that is, it has one double and one single root.

•

$256 d^{3} a^{2} + (- 128 d^{2} b^{2} + 144 d c^{2} b - 27 c^{4}) a + 16 d b^{4} - 4 c^{2} b^{3} = 0$ and none of the polynomials in the last list are zero, including $a \neq 0$ $, c \neq 0$ , $d \neq 0$ , $32 b d - 9 c^{2} \neq 0$ , and $2 b^{3} - (8 b d - 9 c^{2}) a \neq 0$ . The first polynomial is a factor of the discriminant of $p$ :

>	discrim(p, x);

$a (256 d^{3} a^{2} - 128 a b^{2} d^{2} + 144 a b c^{2} d - 27 a c^{4} + 16 b^{4} d - 4 b^{3} c^{2})$

(6.2.4)

The condition that it must be zero implies that $p$ has at least a double root, and the remaining inequations ensure that all the other roots are simple.

6.3 Real root classification and border polynomial

Real root classification

A semi-algebraic system (SAS) is a polynomial system containing polynomial equations $p = 0$ , inequations $p \neq 0$ , and inequalities $p > 0$ or $p \geq 0$ . If a semi-algebraic system contains parameters, it is called a parametric SAS. For a parametric SAS and a prescribed non-negative integer $n$ , the real root classification problem is to compute conditions on the parameters such that the system has exactly $n$ distinct real solutions. In the RegularChains package, a semi-algebraic system is given by four lists $F, N, P, H$ of polynomials, representing the equations, non-negativity conditions (weak inequalities), positivity conditions (strict inequalities), and equations, respectively. The RealRootClassification command takes these four lists as the first four arguments. As in the examples above, the fifth argument $, d$ , is an integer specifying that the last $d$ indeterminates of the underlying polynomial ring are the parameters of the system. The sixth argument, $n$ , specifies the desired number of distinct real solutions, and the last argument is the polynomial ring $R$ . For example,

>	R := PolynomialRing([x, a, b, c]): F := [x^2a+xb+c]: N := []:P := [x]:H := [a]: rr := RealRootClassification(F, N, P, H, 3, 2, R);

$rr ≔ [[regular_semi_algebraic_set], border_polynomial]$

(6.3.1.1)

Here, $a, b, c$ are viewed as parameters and the task is to find conditions for the general polynomial of degree $2$ to have exactly $2$ distinct positive real solutions.

The output is a list, ${rr}_{1}$ of regular semi-algebraic sets and a border polynomial object (BP), ${rr}_{2}$ ; together these describe a first-order logic formula. More precisely, the list ${rr}_{1}$ gives necessary and sufficient conditions for the input system to have the prescribed number of real solutions, provided that the condition encoded by the BP holds. The role of the border polynomial is to exclude degenerate cases, which can be handled later with further computations (more on this later).

The command below shows the contents of the BP, which is actually a list of polynomials. The logical condition encoded by the BP is that none of those polynomials should be zero.

>	Info(rr[2], R);

$[a, b, c, - 4 a c + b^{2}]$

(6.3.1.2)

The following commands extract information about the list of regular SAS and display the logical condition that these encode, using the command Display.

>	ss := rr[1][1]; Display(ss, R);

$ss ≔ regular_semi_algebraic_set$

$\{\begin{array}{c} c < 0 and a < 0 and b > 0 and 4 a c - b^{2} < 0 \\ or c > 0 and a > 0 and b < 0 and 4 a c - b^{2} < 0 \end{array}$

(6.3.1.3)

The output of the last command is to be interpreted as the logical disjunction ("or") of the conditions in the last two rows; each row represents the logical conjunction ("and") of the individual inequalities.

As a first approximation, a regular SAS is a subset of the real solutions of a regular chain. More precisely, it is encoded by three pieces of data:

•	a quantifier free formula $qf$ ,

•	a regular chain $rc$ , and

•	a list of lists of indices $l$ .

This encoding is called a parametric box. The formula $qf$ specifies the conditions on the variables which are free in $rc$ (non-algebraic), whereas $rc$ and $l$ specify the conditions on the variables which are algebraic in $rc$ .

In the example above, $rc$ and $l$ are empty (which is often the case); therefore, the conditions encoded by the parametric box are just the ones of $qf &period;$

>	rc := RepresentingChain(ss, R); Equations(rc, R); box := RepresentingBox(ss); l := RepresentingRootIndex(box);

$rc ≔ regular_chain$

$[]$

$box ≔ parametric_box$

$l ≔ []$

(6.3.1.4)

Summarizing, the answer to the particular question is: provided that the BP holds, the input SAS has exactly $2$ distinct positive real solutions if and only if $qf$ holds. Note that, in the example, the conditions of $qf$ imply those of the BP; thus the answer can be simplified further as follows:

the input SAS has exactly $2$ distinct positive real solutions if and only if $qf$ holds.

Setting infolevel to a nonzero value also provides information on how to interpret the above result.

>	infolevel[RegularChains] := 1: RealRootClassification(F, N, P, H, 3, 2, R); infolevel[RegularChains] := 0:

RealRootClassification: FINAL RESULT:
RealRootClassification: The system has given number of real solution(s) IF AND ONLY IF
RealRootClassification: [R[1] < 0 R[2] < 0 0 <= R[3] R[4] < 0]
RealRootClassification: OR
RealRootClassification: [0 < R[1] 0 < R[2] R[3] <= 0 R[4] < 0]
RealRootClassification: where
RealRootClassification: R[1] = c
RealRootClassification: R[2] = a
RealRootClassification: R[3] = b
RealRootClassification: R[4] = 4*a*c-b^2
RealRootClassification: PROVIDED THAT
RealRootClassification: c <> 0
RealRootClassification: a <> 0
RealRootClassification: 4*a*c-b^2 <> 0
RealRootClassification: 0.5e-2*seconds

$[[regular_semi_algebraic_set], border_polynomial]$

(6.3.1.5)

There are two alternative ways to obtain the answer in a simpler (with less commands) but more coarse manner. The first one, shown above, is to set the infolevel to $1$ . The second one, shown below, is to use the Info command on the regular SAS. With this second approach, one can directly see the encoding of the regular SAS, that is, the parametric box; however, the conditions defining the regular SAS are not shown explicitly.

>	Info(ss, R);

$[[[c, a, b, 4 a c - b^{2}], [[−1, −1, 1, −1], [1, 1, −1, −1]]], [], []]$

(6.3.1.6)

More about border polynomials

The border polynomial is actually a list of polynomials such that the union of the hypersurfaces given by these polynomials contains all the "abnormal" (non-generic) parameter values of the parametric input SAS. For the previous example, the border polynomial is $[$ $c, a, 4 a c - b^{2}]$ . It is intuitive that the zeroes of these polynomials represent degenerate cases: the first two indicate vanishing of the trailing and leading coefficients, respectively, and the last polynomial is the discriminant of the input polynomial. You can use the Info command to view the contents of the BP:

>	Info(rr[2], R);

$[a, b, c, - 4 a c + b^{2}]$

(6.3.2.1)

You might want to go further in analyzing the number of real solutions of the input SAS; that is, you might want to figure out what the answer is in these degenerate cases. This is very simple: it suffices to add the product of the polynomials in the BP to the list of equations when calling RealRootClassification (or to add the polynomials in the BP one by one to the equations of the system and call RealRootClassification accordingly).

In the example above, the question was for which parameter values does the input system has $2$ distinct positive real solutions. When adding the polynomial $4 a c - b^{2}$ to the input system, below, you see that the list of regular SAS in the output from RealRootClassification is empty. This just means that it is impossible that the equation $a x^{2} + b x + c = 0$ has $2$ distinct positive real solutions when its discriminant is zero.

>	rr := RealRootClassification([4ac-b^2, op(F)], N, P, H, 3, 2, R); bp := rr[2];Info(bp, R);

$rr ≔ [[], border_polynomial]$

$bp ≔ border_polynomial$

$[b, a]$

(6.3.2.2)

The above output means: provided that $b \neq 0 \neq c$ holds, the new system never admits two distinct positive real solutions.

More about regular semi-algebraic sets

A regular semi-algebraic set is essentially a subset of the real solutions of a regular chain. When this subset is finite, it can be encoded as a list of points with real coordinates. Each of these real coordinates must be an algebraic number defined by a regular chain and an isolating interval. Such an encoding is called a numerical box.

When this subset is not finite, the idea is to use a regular chain $rc$ that is not zero-dimensional. Such a regular chain will have variables which are free in $rc$ , that is, variables that do not belong to the main variables of the polynomials defining $rc$ . Those free variables can be regarded as parameters. For $rc$ to have real solutions, those free variables may need to satisfy some inequalities; this is where the quantifier-free formula $qf$ comes in. More technically, $qf$ defines an open semi-algebraic set (SAS) at every point of which $rc$ specializes well and separates well (again, viewing the free variables of $rc$ as parameters). This technical condition has the following important consequence: the number of real roots of $rc$ is constant on the open SAS; moreover, each root of $rc$ (complex or real) is a simple root on the open SAS. This in turn has the following consequence: the real solutions of $rc$ can be indexed uniformly on the open SAS. Now this is where the list of lists of indices $l$ comes in; it selects some of the real roots of $rc$ .

To illustrate the concept of a regular semi-algebraic set in a sufficiently general situation, the original question will be modified as follows: Determine under which conditions on $a, b, c, d$ the equation $x^{2} + d = 0$ has 2 distinct real solutions subject to $d$ satisfying $a d^{2} + b d + c = 0$ . Note that $d$ plays a special role among the parameters, so it is ordered before $a, b, c$ .

>	$F := [x^{2} + d, a d^{2} + b d + c] &colon;$ $N ≔ [] &colon;$ $P ≔ [] &colon;$ $H ≔ [] &colon;$ $R ≔ PolynomialRing ([x, d, a, b, c]) &colon;$ $rr ≔ RealRootClassification (F, N, P, H, 4, 2, R)$

$rr ≔ [[regular_semi_algebraic_set], border_polynomial]$

(6.3.3.1)

Recall that the output consists of a list of regular SAS and a BP object. The command RepresentingBox applied to a regular SAS returns its encoding as either a numerical box or a parametric box; the command IsParametricBox will determine which of these two cases holds.

>	$\{seq (IsParametricBox (RepresentingBox (ss, R)), ss = {rr}_{1})\}$

$\{true\}$

(6.3.3.2)

Display the BP using the Info command.

>	$Info ({rr}_{2}, R)$

$[d]$

(6.3.3.3)

Select the first parametric box for further inspection. The command Info shows the raw data defining the parametric box, whereas the command Display pretty prints the conditions encoded by this data.

>	$ss ≔ {rr}_{1}_{1} &colon; box ≔ RepresentingBox (ss, R) &semi;$ $Info (box, R) &semi;$ $Display (box, R) &semi;$

$box ≔ parametric_box$

$[[[d], [[−1]]], [a d^{2} + b d + c], [[1]]]$

$\{\begin{array}{c} c = - a d^{2} - b d \\ d < 0 \end{array}$

(6.3.3.4)

The three pieces of data defining a parametric box can also be accessed independently as shown below. Note that RepresentingChain applies directly to a regular SAS since both numerical boxes and parametric boxes have a representing regular chain. The components $qf$ and $l$ , however, are specific to parametric boxes.

(Note that objects of types regular_chain, regular_system, constructible_set, quantifier_free_formula, parametric_box, regular_semi_algebraic_set, and border_polynomial can all be printed by the Info and Display command.)

>	$box ≔ RepresentingBox (ss, R) &semi;$ $rc ≔ RepresentingChain (ss, R) &semi;$ $qf ≔ RepresentingQuantifierFreeFormula (box) &semi;$ $l ≔ RepresentingRootIndex (box) &semi;$ $Display (rc, R) &semi; Display (qf, R) &semi;$ $Display (box, R)$

$box ≔ parametric_box$

$rc ≔ regular_chain$

$qf ≔ quantifier_free_formula$

$l ≔ [[1]]$

$a d^{2} + b d + c = 0$

$d < 0$

$\{\begin{array}{c} c = - a d^{2} - b d \\ d < 0 \end{array}$

(6.3.3.5)

An advanced example of real root classification

Define a polynomial ring.

>	$unassign (' p',' q') &semi; R := PolynomialRing ([x, y, z, q, p]) &colon;$

Consider the following equations...

>	$g_{1} ≔ y^{2} + z^{2} - 2 y z p - 1 &colon;$ $g_{2} ≔ z^{2} + x^{2} - 2 z x q - 1 &colon;$ $g_{3} ≔ x^{2} + y^{2} - 2 x y q - 1 &colon;$ $g_{4} ≔ 5 q^{2} - 1 &colon;$ $F ≔ [g_{1}, g_{2}, g_{3}, g_{4}] &colon;$

...and constraints.

>	$N := [] &colon;$ $P ≔ [x, y, z, 1 - p^{2}, 1 - q^{2}, p + 1 - 2 q^{2}] &colon;$ $H := [] &colon;$

The last 2 unknowns, $p$ and $q$ , are viewed as parameters. Suppose you want to know the conditions for the system to have no real solutions.

>	$rr ≔ RealRootClassification (F, N, P, H, 2, 0, R)$

$rr ≔ [[regular_semi_algebraic_set, regular_semi_algebraic_set], border_polynomial]$

(6.3.4.1)

>	$Display ({rr}_{1}, R) &semi;$ $Display ({rr}_{2}, R)$

$[\{\begin{array}{c} q = RootOf (5 {_Z}^{2} - 1, index = {real}_{1}) \\ p > 0 and 2 p - 1 > 0 and 5 p - 3 \neq 0 and 5 p + 3 > 0 and 5 p^{2} - 1 > 0 \end{array}, \{\begin{array}{c} q = RootOf (5 {_Z}^{2} - 1, index = {real}_{2}) \\ p > 0 and 2 p - 1 > 0 and 5 p - 3 > 0 and 5 p + 3 > 0 and 5 p^{2} - 1 > 0 \end{array}]$

$[p, 2 p - 1, 5 p - 3, 5 p + 3, 2 p + 1, 5 p^{2} - 1]$

(6.3.4.2)

The output should be read as follows: provided that not one of the polynomials in the BP (displayed in the last row) vanishes, the input system has no real solutions if and only if one of the following two conditions holds:

(1) $p > 0 \land 5 p + 3 > 0 \land 2 p - 1 > 0 \land 5 p^{2} - 1 > 0 \land q$ is the first (from left to right) root of $5 q^{2} - 1$ ;

(2) $p > 0$ $\land 5 p - 3 > 0 \land 5 p + 3 > 0 \land 2 p - 1 > 0 \land 5 p^{2} - 1 > 0 \land q$ is the second root of $5 q^{2} - 1$ .

Why do border polynomials form a type?

Consider again the parametric system in one variable and three parameters, associated with a generic univariate polynomial equation of degree 2.

>	$R := PolynomialRing ([x, a, b, c]) &colon;$ $F ≔ [a x^{2} + b x + c] &colon;$ $N ≔ [] &colon;$ $P ≔ [] &colon;$ $H ≔ [] &colon;$

The border polynomial is intrinsically associated with this system, independent of the prescribed number of real solutions of interest. The BP can be computed directly by the BorderPolynomial command.

>	$bp := BorderPolynomial (F, N, P, H, 3, R) &colon;$ $Info (bp, R)$

$[a, - 4 a c + b^{2}]$

(6.3.5.1)

The reason why border polynomials form a type in the RegularChains package (and cannot just be viewed as lists of polynomials) is that under special circumstances, the border polynomial of a parametric semi-algebraic system may degenerate.

The first such case is when the parameters do not appear in the system of polynomials; then the border polynomial is $[1]$ , as in the next example.

>	$bp ≔ BorderPolynomial ([x^{2} - 1], [], [], [], 3, R) &colon;$ $eval (bp) &semi; Info (bp, R)$

$table ([type = border_polynomial, PolynomialList = [1], PolynomialRing = polynomial_ring])$

$[1]$

(6.3.5.2)

Another special circumstance is that of overdetermined or inconsistent systems; then the border polynomial is denoted as $[]$ , as in the example below.

>	$F := [a x^{2} + b x + c, a, b] &colon;$ $bp ≔ BorderPolynomial (F, [], [], [], 1, R) &colon;$ $eval (bp) &semi; Info (bp, R)$

$table ([type = border_polynomial, PolynomialList = [], PolynomialRing = polynomial_ring])$

$[]$

(6.3.5.3)

Finally, the last special case is when the input system has "generically" infinitely many complex solutions;

then the border polynomial is denoted as $[0],$ as in the example below (this is because only two of the variables, $b$ and $c$ , are considered as parameters).

>	$F := [a x^{2} + b x + c] &colon;$ $bp ≔ BorderPolynomial (F, [], [], [], 2, R) &colon;$ $eval (bp) &semi; Info (bp, R)$

$table ([type = border_polynomial, PolynomialList = [0], PolynomialRing = polynomial_ring])$

$[0]$

(6.3.5.4)

6.4 Solving parametric semi-algebraic systems with application to the study stability of biological systems

The biological system is described by the following system of differential equations.
Its right hand side encodes the equilibria:

>	ode := {diff(x(t),t) = -x(t)+s/(1+y(t)^2), diff(y(t),t)=-y(t)+s/(1+x(t)^2)}; F := [-x+s/(1+y^2), -y+s/(1+x^2)];

$ode ≔ \{\frac{&DifferentialD;}{&DifferentialD; t} x (t) = - x (t) + \frac{s}{1 + {y (t)}^{2}}, \frac{&DifferentialD;}{&DifferentialD; t} y (t) = - y (t) + \frac{s}{1 + {x (t)}^{2}}\}$

$F ≔ [- x + \frac{s}{y^{2} + 1}, - y + \frac{s}{x^{2} + 1}]$

(6.4.1)

The following two Hurwitz determinants determine the stability of hyperbolic equilibria:

>	D1 := -(diff(F[1],x)+diff(F[2],y)); #D1 is 2 D2 := diff(F[1],x)diff(F[2],y)-diff(F[1],y)diff(F[2],x);

$D1 ≔ 2$

$D2 ≔ 1 - \frac{4 s^{2} y x}{{(y^{2} + 1)}^{2} {(x^{2} + 1)}^{2}}$

(6.4.2)

The semi-algebraic system below encodes the asymptotically stable hyperbolic equilibria:

>	P := [numer(normal(F[1]))=0, numer(normal(F[2]))=0, x>0, y>0, s>0, numer(D2)>0];

$P ≔ [- x y^{2} + s - x = 0, - x^{2} y + s - y = 0, 0 < x, 0 < y, 0 < s, 0 < x^{4} y^{4} + 2 x^{4} y^{2} + 2 x^{2} y^{4} - 4 s^{2} y x + x^{4} + 4 x^{2} y^{2} + y^{4} + 2 x^{2} + 2 y^{2} + 1]$

(6.4.3)

Compute a real comprehensive triangular decomposition of P w.r.t. the parameter $s$ :

>	R := PolynomialRing([y, x, s]): ctd := RealComprehensiveTriangularize(P, 1, R);

$ctd ≔ [[[1, squarefree_semi_algebraic_system], [2, squarefree_semi_algebraic_system]], [[semi_algebraic_set, []], [semi_algebraic_set, [1]], [semi_algebraic_set, [2]]]]$

(6.4.4)

Derive the values of $s$ s.t. P has $2$ positive real solutions, that is the biological system is bistable:

>	ctd2 := RealComprehensiveTriangularize(ctd, R, 2); Display(ctd2[2][1][1],R);

$ctd2 ≔ [[[1, squarefree_semi_algebraic_system]], [[semi_algebraic_set, [1]]]]$

$[2 < s]$

(6.4.5)

The two asymptotically stable equilibria are represented by a $squarefree_semi_algebraic_system &period;$

>	ss := ctd2[1][1][2]; Display(ss, R);

$ss ≔ squarefree_semi_algebraic_system$

$\{\begin{array}{c} x y - 1 = 0 \\ x^{2} - s x + 1 = 0 \\ y > 0 \\ x > 0 \\ s^{7} x - s^{6} - 6 x s^{5} + 5 s^{4} + 8 x s^{3} - 4 s^{2} > 0 \end{array}$

(6.4.6)

7. FFT-based polynomial arithmetic

The module FastArithmeticTools supports the implementation of modular methods for computing with polynomials, algebraic extensions, and thus regular chains. This support consists of fundamental operations such as resultant, polynomial gcds, normal forms, etc.

The commands of this module work in prime characteristic and rely on asymptotically fast algorithms. Most of the underlying polynomial arithmetic is performed by C code and relies on (multi-dimensional) Fast Fourier Transform (FFT) and straight line programs (SLPs). This C code is highly optimized and implements the Truncated Fourier Transform (TFT) and an improved version of Montgomery's trick.

•	The commands IteratedResultantDim0 and IteratedResultantDim1 compute the iterated resultant of a polynomial with regards to a regular chain of dimension 0 and 1, respectively.

•	The commands NormalFormDim0 and ReduceCoefficientsDim0 compute the normal form of a polynomial with regards to a zero-dimensional regular chain.

•	The commands NormalizePolynomialDim0 and NormalizeRegularChainDim0 normalize a polynomial (with regards to a zero-dimensional regular chain) and a regular chain (with regards to itself).

•	The command RegularizeDim0 tests whether a polynomial is invertible modulo a zero-dimensional regular chain.

•	The commands RegularGcdBySpecializationCube, ResultantBySpecializationCube, and SubresultantChainSpecializationCube compute resultants and polynomial gcds modulo a regular chain using fast evaluation and interpolation.

•	The commands RandomRegularChainDim0 and RandomRegularChainDim1 compute random regular chains of given degrees.

•	Finally, the command BivariateModularTriangularize solves bivariate polynomial systems.

Most of the commands of FastArithmeticTools implements core operations on regular chains such as regularity test and polynomial gcd modulo a regular chain. However, these commands have several constraints. On top of the characteristic constraint (detailed below), the current regular chain must have dimension zero or one. There is only one exception: the command RegularGcdBySpecializationCube requires no assumption about the dimension. Note also that some commands do not take any regular chains as input (for instance, SubresultantChainSpecializationCube and ResultantBySpecializationCube).

Since FastArithmeticTools relies heavily on direct FFT and Montgomery's trick, the characteristic of the polynomial ring must be a prime number $p$ satisfying the following properties. First, it should not be greater than $962592769$ . Secondly, the number $p - 1$ should be divisible by a sufficiently large power of 2. The power $2^{20}$ is often sufficient. If this power of 2 is not large enough, then an appropriate error message is returned. Using try-catch statements is highly recommended when programming with the commands of this submodule.

7.1 The impact of fast arithmetic

The purpose of the session below is to demonstrate how resultants can be computed with the module FastArithmeticTools.

>	restart; with(RegularChains); with(FastArithmeticTools);

(7.1.1)

Set the polynomial ring.

>	p := 962592769; vars := [a, b]; R := PolynomialRing(vars, p);

$p ≔ 962592769$

$vars ≔ [a, b]$

$R ≔ polynomial_ring$

(7.1.2)

Define the input polynomials.

>	$degbound := 40 &semi;$ $f_{1} ≔ randpoly (vars, dense, degree = degbound) mod p &colon;$ $f_{2} ≔ randpoly (vars, dense, degree = degbound) mod p &colon;$

$degbound ≔ 40$

(7.1.3)

Evaluate/interpolate by subproduct tree techniques.

>	$t_{1} := time () &colon;$ $SCube ≔ SubresultantChainSpecializationCube (f [1], f [2], a, R, 0) &semi;$ $r_{1} := ResultantBySpecializationCube (f [1], f [2], a, SCube, R) &colon;$ $t_{1} ≔ time () - t_{1}$

$SCube ≔ subresultant_chain_specialization_cube$

$t_{1} ≔ 0.263$

(7.1.4)

Evaluate/interpolate by multi-dimensional TFT.

>	$t_{2} ≔ time () &colon;$ $SCube := SubresultantChainSpecializationCube (f_{1}, f_{2}, a, R, 1) &semi;$ $r_{2} ≔ ResultantBySpecializationCube (f_{1}, f_{2}, a, SCube, R) &colon;$ $t_{2} ≔ time () - t_{2}$

$SCube ≔ subresultant_chain_specialization_cube$

$t_{2} ≔ 0.166$

(7.1.5)

Evaluate/interpolate without fast arithmetic.

>	$t_{3} ≔ time () &colon;$ $r_{3} ≔ Resultant (f_{1}, f_{2}, a) mod p &colon;$ $t_{3} ≔ time () - t_{3}$

$t_{3} ≔ 12.306$

(7.1.6)

Compare the results.

>	$evalb (`and` (r_{1} = r_{2}, r_{2} = r_{3})) &semi;$ $degree (r_{1}) &semi;$ $t_{1}, t_{2}, t_{3} &semi;$

$true$

$1600$

$0.263, 0.166, 12.306$

(7.1.7)

7.2 The impact of modular methods together with fast arithmetic

The session below shows how polynomial gcds modulo regular chains can be computed with the module FastArithmeticTools. This applies also to polynomial gcds over towers of field extensions.

>	restart: with(RegularChains): with(FastArithmeticTools): with(ChainTools):

Set the polynomial ring.

>	p := 962592769: vars := [a, b, c]: R := PolynomialRing(vars, p):

Define the polynomials.

>	$degbound := 7 &semi;$ $f_{1} ≔ randpoly (vars, dense, degree = degbound) mod p &colon;$ $f_{2} ≔ randpoly (vars, dense, degree = degbound) mod p &colon;$

$degbound ≔ 7$

(7.2.1)

Evaluate/interpolate by multi-dimensional TFT.

$t_{2} := time () &colon;$ $SCube ≔ SubresultantChainSpecializationCube (f [1], f [2], a, R, 1) &colon;$ $r_{2} := ResultantBySpecializationCube (f [1], f [2], a, SCube, R) &colon;$ $rc := Chain ([r_{2}], Empty (R), R) &colon;$ $g_{2} := RegularGcdBySpecializationCube (f [1], f [2], rc, SCube, R) &colon;$ $t_{2} := time () - t_{2}$

$t_{2} ≔ 0.237$

(7.2.2)

Compute without fast arithmetic and without modular methods (evaluation / interpolation). Since the command RegularGcd tries to use modular methods whenever possible, change the characteristic of the ring to a small prime, so as to enforce the use of non-modular and non-FFT-based algorithms.

>	$t_{3} ≔ time () &colon;$ $S ≔ PolynomialRing (vars, 257) &colon;$ $r_{3} ≔ Resultant (f_{1}, f_{2}, a) mod p &colon;$ $rc ≔ Chain ([r_{3}], Empty (S), S) &colon;$ $g_{3} ≔ RegularGcd (f_{1}, f_{2}, a, rc, S) &colon;$ $t ≔ time () - t_{3}$

$t ≔ 12.734$

(7.2.3)

7.3 Accelerating the core operations of the RegularChains library

Testing whether a polynomial is invertible modulo (the saturated ideal of) a zero-dimensional regular chain is a fundamental operation. Its implementation RegularizeDim0 in the FastArithmeticTools takes advantage of FFT-based polynomial arithmetic and improves on the command Regularize. Note that this latter command implements a more general algorithm (with no assumptions on characteristic or dimension).

>	restart: with(RegularChains): with(FastArithmeticTools): with(ChainTools):

Set the polynomial ring.

>	p := 962592769: vars := [x1, x2, x3, x4]: R := PolynomialRing(vars, p):

Define a random (dense) regular chain and a polynomial.

>	N := nops(vars): d := 5; degrees := [3, 4, 5, 8]; rc := RandomRegularChainDim0(vars, degrees, p): p := `mod`(`mod`(randpoly(vars, dense, degree = d)+rand(), p), p):

$d ≔ 5$

$degrees ≔ [3, 4, 5, 8]$

(7.3.1)

Compute with the modular code.

>	$t_{1} := time () &colon;$ $RegularizeDim0 (p, rc, R) &colon;$ $t_{1} := time () - t_{1}$

$t_{1} ≔ 0.183$

(7.3.2)

Compute with the generic algorithm (without modular methods and without asymptotically fast arithmetic). Here again, change the characteristic to prevent use of the fast code.

>	$t_{2} ≔ time () &colon;$ $Regularize (p, rc, PolynomialRing (vars, 17)) &semi;$ $t_{2} ≔ time () - t_{2}$

$[[rc], []]$

$t_{2} ≔ 0.344$

(7.3.3)

7.4 Solving large bivariate systems

In this example, solve a dense bivariate and square system which has 2500 solutions.

>	restart: with(RegularChains): with(FastArithmeticTools): with(ChainTools):

Set the polynomial ring.

>	p := 469762049: vars := [x, y]: R := PolynomialRing(vars, p):

Define the polynomials and inspect their number of terms.

>	$degbound := 50 &colon;$ $f_{1} ≔ randpoly (vars, dense, degree = degbound) mod p &colon;$ $f_{2} ≔ randpoly (vars, dense, degree = degbound) mod p &colon;$ $nops (f_{1}) &semi;$ $nops (f_{2})$

$1314$

$1321$

(7.4.1)

Solve the system given by $f_{1} = f_{2} = 0$ using fast arithmetic.

>	$t := time () &colon;$ $l ≔ BivariateModularTriangularize ([f_{1}, f_{2}], R) &colon;$ $t ≔ time () - t$

$t ≔ 0.691$

(7.4.2)

Check the number of solutions and its number of terms.

>	map(NumberOfSolutions, l, R); map(nops, Equations(l[1], R));

$[2500]$

$[2404, 2501]$

(7.4.3)

8. Cylindrical algebraic decomposition and quantifier elimination

The RegularChains library provides a set of commands for computing cylindrical algebraic decomposition (CAD, see command
CylindricalAlgebraicDecompose) and doing quantifier elimination (QE, see command QuantifierElimination). The underlying algorithm first computes a cylindrical decomposition of the complex space (CCD, see command CylindricalDecompose), which is further refined into a cylindrical algebraic decomposition of the real space. Such an algorithm is different from the traditional projection-lifting algorithm introduced by George Collins and further refined by others.

>	with(ConstructibleSetTools):

8.1 Cylindrical decomposition of the complex space

A CCD is a partition of the complex space into disjoint cells such that they are cylindrically arranged,

meaning the projection of any two cells onto any lower dimensional space are either identical or disjoint,

and each cell is the zero set of a regular system.

>	R := PolynomialRing([y, x]); F := [x^2+y^2-1];

$R ≔ polynomial_ring$

$F ≔ [x^{2} + y^{2} - 1]$

(8.1.1)

>	ccd := CylindricalDecompose(F, R); Display(ccd, R);

$ccd ≔ c_c_d$

$[\{\begin{array}{c} y = 0 \\ x - 1 = 0 \end{array}, \{\begin{array}{c} x - 1 = 0 \\ y \neq 0 \end{array}, \{\begin{array}{c} y = 0 \\ x + 1 = 0 \end{array}, \{\begin{array}{c} x + 1 = 0 \\ y \neq 0 \end{array}, \{\begin{array}{c} y^{2} + x^{2} - 1 = 0 \\ x^{2} - 1 \neq 0 \end{array}, \{\begin{array}{c} y^{2} + x^{2} - 1 \neq 0 \\ x^{2} - 1 \neq 0 \end{array}]$

(8.1.2)

>	ccd := CylindricalDecompose(F, R, output=system); Display(ccd, R);

$ccd ≔ [regular_system, regular_system, regular_system, regular_system, regular_system, regular_system]$

(8.1.3)

A CCD is best described by a complex cylindrical tree (CCT).

Informally, a CCT T of <[x_1,x_2,...,x_n] is a rooted tree with each non-root node
described by a polynomial constraint p(x_1,...,x_i)=0, p(x_1,...,x_i) $\neq 0,$ or "any x_i".

Moreover, the polynomial constraints in any path of T form a regular system such that
the union of their zero sets form a CCD.

The CCT can be displayed by the piecewise option.

>	R := PolynomialRing([y, x]); F := [x^2+y^2-1]; CylindricalDecompose(F, R, output=piecewise);

$R ≔ polynomial_ring$

$F ≔ [y^{2} + x^{2} - 1]$

$\{\begin{array}{c} \{\begin{array}{c} 1 & y = 0 \\ 1 & otherwise \end{array} & x - 1 = 0 \\ \{\begin{array}{c} 1 & y = 0 \\ 1 & otherwise \end{array} & x + 1 = 0 \\ \{\begin{array}{c} 1 & y^{2} + x^{2} - 1 = 0 \\ 1 & otherwise \end{array} & otherwise \end{array}$

(8.1.4)

If a set of polynomials, say F, is passed to CylindricalDecompose, then an F-invariant CCD

is computed. By F-invariant, we mean for any polynomial f of F and any cell C in the CCD,

either C is contained in the zero set of f or C has no intersection with the zero set of f.

One can use the command IsContained and Intersection to verify if the output satisfies the

invariance property.

>	R := PolynomialRing([y, x]); f:= x^2+y^2-1; g := x*y-1; F := [f,g]; lrs := CylindricalDecompose(F, R, output=system);

$R ≔ polynomial_ring$

$f ≔ y^{2} + x^{2} - 1$

$g ≔ x y - 1$

$F ≔ [y^{2} + x^{2} - 1, x y - 1]$

$lrs ≔ [regular_system, regular_system, regular_system, regular_system, regular_system, regular_system, regular_system, regular_system, regular_system, regular_system, regular_system, regular_system, regular_system, regular_system, regular_system]$

(8.1.5)

For example, the first complex cell, namely the zero set of the first regular system, is contained

in the zero set of f.

>	rs1 := lrs[1]; Display(rs1, R); csf := Triangularize([f], [1], R, output=lazard); cs1 := ConstructibleSet([rs1], R); IsContained(cs1, csf, R);

$rs1 ≔ regular_system$

$\{\begin{array}{c} y - 1 = 0 \\ x = 0 \end{array}$

$csf ≔ constructible_set$

$cs1 ≔ constructible_set$

$true$

(8.1.6)

This cell has no intersection with the zero set of g.

>	rs1 := lrs[1]; Display(rs1, R); csg := Triangularize([g], [1], R, output=lazard); cs1 := ConstructibleSet([rs1], R); cs := Intersection(cs1, csg, R); IsEmpty(cs, R);

$rs1 ≔ regular_system$

$\{\begin{array}{c} y - 1 = 0 \\ x = 0 \end{array}$

$csg ≔ constructible_set$

$cs1 ≔ constructible_set$

$cs ≔ constructible_set$

$true$

(8.1.7)

Besides a list of polynomials, a list of constraints (equations or inequations) are also allowed

to be the input of CylindricalDecompose. In this case, instead of a sign-invariant CCD, a (smaller) truth-invariant CCD will be computed.

>	R := PolynomialRing([y, x]); F := [f=0,g=0]; cad := CylindricalDecompose(F, R, output=piecewise);

$R ≔ polynomial_ring$

$F ≔ [y^{2} + x^{2} - 1 = 0, y x - 1 = 0]$

$cad ≔ \{\begin{array}{c} \{\begin{array}{c} 1 & y x - 1 = 0 \\ 1 & otherwise \end{array} & x^{4} - x^{2} + 1 = 0 \\ 1 & otherwise \end{array}$

(8.1.8)

Several other output formats are supported, which are useful in some context.
If the input is a list of constraints, the options ''output=cct" and "output=tree" generate a partial CCT expressing exactly their complex zeros while the options "output=ccd", "output=fulltree" and "output=piecewise" will generate a complete CCT.
The options "output=tree" and "output=fulltree" represent the CCT in a nested list style;

R := PolynomialRing([y, x]);
F := [f=0,g=0];
cad := CylindricalDecompose(F, R, output=cct);
Display(cad, R);
cad := CylindricalDecompose(F, R, output='ccd');
Display(cad, R);
cad := CylindricalDecompose(F, R, output='piecewise');
cad := CylindricalDecompose(F, R, output='tree');
cad := CylindricalDecompose(F, R, output='fulltree');

$R ≔ polynomial_ring$

$F ≔ [y^{2} + x^{2} - 1 = 0, y x - 1 = 0]$

$cad ≔ c_c_t$

$[\{\begin{array}{c} x y - 1 = 0 \\ x^{4} - x^{2} + 1 = 0 \\ x \neq 0 \end{array}]$

$cad ≔ c_c_d$

$[\{\begin{array}{c} x y - 1 = 0 \\ x^{4} - x^{2} + 1 = 0 \\ x \neq 0 \end{array}, \{\begin{array}{c} x^{4} - x^{2} + 1 = 0 \\ x y - 1 \neq 0 \end{array}, x^{4} - x^{2} + 1 \neq 0]$

$cad ≔ \{\begin{array}{c} \{\begin{array}{c} 1 & y x - 1 = 0 \\ 1 & otherwise \end{array} & x^{4} - x^{2} + 1 = 0 \\ 1 & otherwise \end{array}$

$cad ≔ [1, [x^{4} - x^{2} + 1 = 0, [y x - 1 = 0]]]$

$cad ≔ [1, [x^{4} - x^{2} + 1, [y x - 1], [1]], [1, [1]]]$

(8.1.9)

8.2 Cylindrical algebraic decomposition of the real space

The command CylindricalAlgebraicDecompose is used to compute a cylindrical algebraic decomposition of the real space. It supports different inputs, like list of polynomials, list of polynomial constraints, or list of list of polynomial constraints. It also provides different output formats, such as 'output'='piecewise', 'output'='tree', output='list', output='cadcell', output='rootof', output='cad'. The default output option is output='cad'.

If the input F is a list of polynomials, an F-invariant CAD, which means that any polynomial of F is sign-invariant on any cell of the output CAD, is computed.

>	restart; with(RegularChains): with(SemiAlgebraicSetTools): R := PolynomialRing([y, x]); F := [y^2-x, x-1]; cad := CylindricalAlgebraicDecompose(F, R);

$R ≔ polynomial_ring$

$F ≔ [y^{2} - x, x - 1]$

$cad ≔ c_a_d$

(8.2.1)

By default, the output is a 'c_a_d' type. One can use Display or Info to show its cells.

>	Display(cad, R);

$[\{\begin{array}{c} y = y \\ x < 0 \end{array}, \{\begin{array}{c} y < 0 \\ x = 0 \end{array}, \{\begin{array}{c} y = 0 \\ x = 0 \end{array}, \{\begin{array}{c} 0 < y \\ x = 0 \end{array}, \{\begin{array}{c} y < - \sqrt{x} \\ 0 < x \land x < 1 \end{array}, \{\begin{array}{c} y = - \sqrt{x} \\ 0 < x \land x < 1 \end{array}, \{\begin{array}{c} - \sqrt{x} < y \land y < \sqrt{x} \\ 0 < x \land x < 1 \end{array}, \{\begin{array}{c} y = \sqrt{x} \\ 0 < x \land x < 1 \end{array}, \{\begin{array}{c} \sqrt{x} < y \\ 0 < x \land x < 1 \end{array}, \{\begin{array}{c} y < −1 \\ x = 1 \end{array}, \{\begin{array}{c} y = −1 \\ x = 1 \end{array}, \{\begin{array}{c} −1 < y \land y < 1 \\ x = 1 \end{array}, \{\begin{array}{c} y = 1 \\ x = 1 \end{array}, \{\begin{array}{c} 1 < y \\ x = 1 \end{array}, \{\begin{array}{c} y < - \sqrt{x} \\ 1 < x \end{array}, \{\begin{array}{c} y = - \sqrt{x} \\ 1 < x \end{array}, \{\begin{array}{c} - \sqrt{x} < y \land y < \sqrt{x} \\ 1 < x \end{array}, \{\begin{array}{c} y = \sqrt{x} \\ 1 < x \end{array}, \{\begin{array}{c} \sqrt{x} < y \\ 1 < x \end{array}]$

(8.2.2)

>	Info(cad, R);

$[[x < 0, y = y], [x = 0, y < 0], [x = 0, y = 0], [x = 0, 0 < y], [0 < x \land x < 1, y < - \sqrt{x}], [0 < x \land x < 1, y = - \sqrt{x}], [0 < x \land x < 1, - \sqrt{x} < y \land y < \sqrt{x}], [0 < x \land x < 1, y = \sqrt{x}], [0 < x \land x < 1, \sqrt{x} < y], [x = 1, y < −1], [x = 1, y = −1], [x = 1, −1 < y \land y < 1], [x = 1, y = 1], [x = 1, 1 < y], [1 < x, y < - \sqrt{x}], [1 < x, y = - \sqrt{x}], [1 < x, - \sqrt{x} < y \land y < \sqrt{x}], [1 < x, y = \sqrt{x}], [1 < x, \sqrt{x} < y]]$

(8.2.3)

The tree data structure of the CAD is best shown by the option 'output'='piecewise'.

>	R := PolynomialRing([y, x]); F := [y^2-x]; cad := CylindricalAlgebraicDecompose(F, R, output=piecewise);

$R ≔ polynomial_ring$

$F ≔ [y^{2} - x]$

$cad ≔ \{\begin{array}{c} 1 & x < 0 \\ \{\begin{array}{c} 1 & y < 0 \\ 1 & y = 0 \\ 1 & 0 < y \end{array} & x = 0 \\ \{\begin{array}{c} 1 & y < - \sqrt{x} \\ 1 & y = - \sqrt{x} \\ 1 & - \sqrt{x} < y \land y < \sqrt{x} \\ 1 & y = \sqrt{x} \\ 1 & \sqrt{x} < y \end{array} & 0 < x \end{array}$

(8.2.4)

The input can be a list of polynomial constraints, which incodes a conjunction formula or a semi-algebraic system. The Display command shows all the cells of the output CAD.

>	R := PolynomialRing([y, x]); F := [y^2-x=0, x-1=0]; cad := CylindricalAlgebraicDecompose(F, R); Display(cad, R);

$R ≔ polynomial_ring$

$F ≔ [y^{2} - x = 0, x - 1 = 0]$

$cad ≔ c_a_d$

$[\{\begin{array}{c} y = y \\ x < 1 \end{array}, \{\begin{array}{c} y < −1 \\ x = 1 \end{array}, \{\begin{array}{c} y = −1 \\ x = 1 \end{array}, \{\begin{array}{c} −1 < y \land y < 1 \\ x = 1 \end{array}, \{\begin{array}{c} y = 1 \\ x = 1 \end{array}, \{\begin{array}{c} 1 < y \\ x = 1 \end{array}, \{\begin{array}{c} y = y \\ 1 < x \end{array}]$

(8.2.5)

To see all the true cells, that is the cells satisfying the input constraints, 'output'='cadcell' can be used.

>	R := PolynomialRing([y, x]); cad := CylindricalAlgebraicDecompose([x^2+y^2-1=0, x*y-1/2=0], R, output=cadcell); Display(cad, R);

$R ≔ polynomial_ring$

$cad ≔ [cad_cell, cad_cell]$

$[\{\begin{array}{c} y = \frac{1}{2 x} \\ x = - \frac{\sqrt{2}}{2} \end{array}, \{\begin{array}{c} y = \frac{1}{2 x} \\ x = \frac{\sqrt{2}}{2} \end{array}]$

(8.2.6)

As can be seen, the option 'output'='cadcell' does not attempt to make complete back substitution.

The option 'output'='rootof' instead supports complete back substitution and tries to merge adjacent cells
to produce compact formula.

>	R := PolynomialRing([y, x]): cad := CylindricalAlgebraicDecompose([x^2+y^2-1=0, x*y-1/2=0], R, output=rootof);

$cad ≔ [[x = - \frac{\sqrt{2}}{2}, y = - \frac{\sqrt{2}}{2}], [x = \frac{\sqrt{2}}{2}, y = \frac{\sqrt{2}}{2}]]$

(8.2.7)

>	R := PolynomialRing([y, x]): cad := CylindricalAlgebraicDecompose([x^2+y^2-1<=0], R, output=rootof);

$cad ≔ [[−1 \leq x \land x \leq 1, - \sqrt{- x^{2} + 1} \leq y \land y \leq \sqrt{- x^{2} + 1}]]$

(8.2.8)

The input can also be a list of list of polynomial constraints, which represents a formula in disjunctive normal form or a union of semi-algebraic systems.

>	R := PolynomialRing([y, x]); F := [[y^2-x=0], [x-1=0]]; cad := CylindricalAlgebraicDecompose(F, R, output=cadcell); Display(cad, R);

$R ≔ polynomial_ring$

$F ≔ [[y^{2} - x = 0], [x - 1 = 0]]$

$cad ≔ [cad_cell, cad_cell, cad_cell, cad_cell, cad_cell, cad_cell, cad_cell, cad_cell, cad_cell, cad_cell]$

$[\{\begin{array}{c} y = 0 \\ x = 0 \end{array}, \{\begin{array}{c} y = - \sqrt{x} \\ 0 < x \land x < 1 \end{array}, \{\begin{array}{c} y = \sqrt{x} \\ 0 < x \land x < 1 \end{array}, \{\begin{array}{c} y < −1 \\ x = 1 \end{array}, \{\begin{array}{c} y = −1 \\ x = 1 \end{array}, \{\begin{array}{c} −1 < y \land y < 1 \\ x = 1 \end{array}, \{\begin{array}{c} y = 1 \\ x = 1 \end{array}, \{\begin{array}{c} 1 < y \\ x = 1 \end{array}, \{\begin{array}{c} y = - \sqrt{x} \\ 1 < x \end{array}, \{\begin{array}{c} y = \sqrt{x} \\ 1 < x \end{array}]$

(8.2.9)

Note that only cells making the input formula satisfied are shown. To see all the cells, one

can use the default option 'output'='cad' or the option 'output'='allcell'.

>	R := PolynomialRing([y, x]); F := [[y^2-x=0], [x-1=0]]; cad := CylindricalAlgebraicDecompose(F, R, output=allcell); Display(cad, R);

$R ≔ polynomial_ring$

$F ≔ [[y^{2} - x = 0], [x - 1 = 0]]$

(8.2.10)

Several other output formats are supported, including 'output'='tree', 'output'='list'.

>	R := PolynomialRing([y, x]); F := [y^2-x]; cad := CylindricalAlgebraicDecompose(F, R, output=list);

$R ≔ polynomial_ring$

$F ≔ [y^{2} - x]$

$cad ≔ [[[1, 1], [regular_chain, [[- \frac{1}{2}, - \frac{1}{2}], [0, 0]]]], [[2, 1], [regular_chain, [[0, 0], [- \frac{1}{2}, - \frac{1}{2}]]]], [[2, 2], [regular_chain, [[0, 0], [0, 0]]]], [[2, 3], [regular_chain, [[0, 0], [\frac{1}{2}, \frac{1}{2}]]]], [[3, 1], [regular_chain, [[\frac{1}{2}, \frac{1}{2}], [- \frac{3}{2}, - \frac{3}{2}]]]], [[3, 2], [regular_chain, [[\frac{1}{2}, \frac{1}{2}], [- \frac{1482911}{2097152}, - \frac{741455}{1048576}]]]], [[3, 3], [regular_chain, [[\frac{1}{2}, \frac{1}{2}], [0, 0]]]], [[3, 4], [regular_chain, [[\frac{1}{2}, \frac{1}{2}], [\frac{741455}{1048576}, \frac{46341}{65536}]]]], [[3, 5], [regular_chain, [[\frac{1}{2}, \frac{1}{2}], [\frac{3}{2}, \frac{3}{2}]]]]]$

(8.2.11)

>	R := PolynomialRing([y, x]); F := [y^2-x]; cad := CylindricalAlgebraicDecompose(F, R, output=tree);

$R ≔ polynomial_ring$

$F ≔ [y^{2} - x]$

$cad ≔ [1, [[[−1, x, 1], [regular_chain, [[- \frac{1}{2}, - \frac{1}{2}]]]], [[[1], [regular_chain, [[- \frac{1}{2}, - \frac{1}{2}], [0, 0]]]]]], [[[x, 1], [regular_chain, [[0, 0]]]], [[[−1, y, 1], [regular_chain, [[0, 0], [- \frac{1}{2}, - \frac{1}{2}]]]]], [[[y, 1], [regular_chain, [[0, 0], [0, 0]]]]], [[[1, y, 1], [regular_chain, [[0, 0], [\frac{1}{2}, \frac{1}{2}]]]]]], [[[1, x, 1], [regular_chain, [[\frac{1}{2}, \frac{1}{2}]]]], [[[−1, - y^{2} + x, 1], [regular_chain, [[\frac{1}{2}, \frac{1}{2}], [- \frac{3}{2}, - \frac{3}{2}]]]]], [[[- y^{2} + x, 1], [regular_chain, [[\frac{1}{2}, \frac{1}{2}], [- \frac{1482911}{2097152}, - \frac{741455}{1048576}]]]]], [[[- y^{2} + x, 1, - y^{2} + x, 2], [regular_chain, [[\frac{1}{2}, \frac{1}{2}], [0, 0]]]]], [[[- y^{2} + x, 2], [regular_chain, [[\frac{1}{2}, \frac{1}{2}], [\frac{741455}{1048576}, \frac{46341}{65536}]]]]], [[[1, - y^{2} + x, 2], [regular_chain, [[\frac{1}{2}, \frac{1}{2}], [\frac{3}{2}, \frac{3}{2}]]]]]]]$

(8.2.12)

The default algorithm of CylindricalAlgebraicDecompose is an incremental algorithm (CM) based on triangular decompositions proposed by Changbo Chen and Marc Moreno Maza on ASCM' 2012. An alternative algorithm (CMXY) proposed by Changbo Chen, Marc Moreno Maza, Bican Xia and Lu Yang on ISSAC' 2009, which is also the first generation of algorithm for computing CAD based on triangular decompositions, is available by the option 'method'='recursive'.

>	R := PolynomialRing([y, x]); F := [y^2-x]; cad := CylindricalAlgebraicDecompose(F, R, method=recursive, output=piecewise);

$R ≔ polynomial_ring$

$F ≔ [y^{2} - x]$

(8.2.13)

The CM algorithm can take advantage of equational constraints in a single semi-algebraic system.

That is, if the input is a list of polynomial constraints involving equations, by default, the option

'optimization'='EC' is enabled. To disable such optimization, one sets 'optimization'='false'.

>	R := PolynomialRing([y, x]); F := [x^2+y^2-1=0, y^2-x=0]; cad := CylindricalAlgebraicDecompose(F, R, optimization='EC', output=allcell): nops(cad); cad2 := CylindricalAlgebraicDecompose(F, R, optimization='false', output=allcell): nops(cad2);

$R ≔ polynomial_ring$

$F ≔ [y^{2} + x^{2} - 1 = 0, y^{2} - x = 0]$

$9$

$53$

(8.2.14)

If input is a list of list of polynomial constraints, by default an algorithm (RC-TTICAD) for computing truth-table invariant CAD is used to compute a smaller CAD. One could use the option 'optimization'='TTICAD' or 'optimization'='EC' to enable or disable it.

>	R := PolynomialRing([y, x]); F := [[(x-2)^2+(y-4)^2=0], [x^2+y^2-1=0, y^2-x=0]]; cad := CylindricalAlgebraicDecompose(F, R, optimization='TTICAD', output=allcell): nops(cad); cad2 := CylindricalAlgebraicDecompose(F, R, optimization='EC', output=allcell): nops(cad2);

$R ≔ polynomial_ring$

$F ≔ [[{(x - 2)}^{2} + {(y - 4)}^{2} = 0], [y^{2} + x^{2} - 1 = 0, y^{2} - x = 0]]$

$33$

$65$

(8.2.15)

To get a sample point of a CAD cell, the function SamplePoints can be called. Here no cost occurs since sample points are computed along the computation of the CAD and are stored in the type cad cell. A sample point is encoded by the type box, which is represented by a regular chain and an isolation cube. Such a representation allows one to easily test if the sign of a polynomial at the sample point by calling the function SignAtBox.

>	R := PolynomialRing([y, x]): cad := CylindricalAlgebraicDecompose([x^2+y^2-1=0, x*y-1/2=0],R,output=cadcell); Display(cad, R); sp := map(SamplePoints, cad, R); Display(sp, R); s := SignAtBox(x^2+y^2-2, sp[1], R);

$cad ≔ [cad_cell, cad_cell]$

$[\{\begin{array}{c} y = \frac{1}{2 x} \\ x = - \frac{\sqrt{2}}{2} \end{array}, \{\begin{array}{c} y = \frac{1}{2 x} \\ x = \frac{\sqrt{2}}{2} \end{array}]$

$sp ≔ [box, box]$

$[\{\begin{array}{c} y = [- \frac{46341}{65536}, - \frac{1482907}{2097152}] \\ x = [- \frac{46341}{65536}, - \frac{741455}{1048576}] \end{array}, \{\begin{array}{c} y = [\frac{185363}{262144}, \frac{741457}{1048576}] \\ x = [\frac{741455}{1048576}, \frac{46341}{65536}] \end{array}]$

$s ≔ −1$

(8.2.16)

8.3 Quantifier elimination

As one of the main applications of cylindrical algebraic decomposition, quantifier elimination is fully supported in the library. The function for doing quantifier elimination is QuantifierElimination.

The user interface of Quantifier Elimination relies on some features of the Logic library of Maple.

We introduce in addition the existential quantifier "&E" and the universal quantifier "&A".

Suppose we'd like to solve the following QE problem (due to Davenport and Heintz):

(∃ c) (∀ b, a) ((a = d ∧ b = c) ∨ ( a=c ∧ b = 1)) ⇒a^2 = b.

>	f := &E([c]), &A([b, a]), ((a=d) &and (b=c)) &or ((a=c) &and (b=1)) &implies (a^2=b);

$f ≔ &E ([c]), &A ([b, a]), (((a = d) &and (b = c)) &or ((a = c) &and (b = 1))) &implies (a^{2} = b)$

(8.3.1)

>	out := QuantifierElimination(f);

$out ≔ (d - 1 = 0) &or (d + 1 = 0)$

(8.3.2)

Note that in the previous example, no variable order is specified.

In such case, the function will try to find the best elimination order according to some heuristic strategy.

>	R := PolynomialRing([x, a, b, c]); f := &E([x]), ax^2+bx+c=0; out := QuantifierElimination(f, R);

$R ≔ polynomial_ring$

$f ≔ &E ([x]), a x^{2} + b x + c = 0$

$out ≔ &or (c = 0, &and (c < 0, b < 0, a = 0), &and (c < 0, b < 0, 0 < a), &and (c < 0, b < 0, 4 a c - b^{2} = 0), &and (c < 0, b < 0, a < 0, 4 a c - b^{2} < 0), &and (c < 0, b = 0, 0 < a), &and (c < 0, 0 < b, a = 0), &and (c < 0, 0 < b, 0 < a), &and (c < 0, 0 < b, 4 a c - b^{2} = 0), &and (c < 0, 0 < b, a < 0, 4 a c - b^{2} < 0), &and (0 < c, b < 0, a < 0), &and (0 < c, b < 0, a = 0), &and (0 < c, b < 0, 4 a c - b^{2} = 0), &and (0 < c, b < 0, 0 < a, 4 a c - b^{2} < 0), &and (0 < c, b = 0, a < 0), &and (0 < c, 0 < b, a < 0), &and (0 < c, 0 < b, a = 0), &and (0 < c, 0 < b, 4 a c - b^{2} = 0), &and (0 < c, 0 < b, 0 < a, 4 a c - b^{2} < 0))$

(8.3.3)

By default, QuantifierElimination returns the standard quantifier free formula, namely Tarski formula.

Extended Tarski formula is supported by the option 'output'='rootof'.

>	f := &E([x]), ax^2+bx+c=0; out := QuantifierElimination(f,'output'='rootof');

$f ≔ &E ([x]), x^{2} a + x b + c = 0$

$out ≔ &or ((a < 0) &and (\frac{b^{2}}{4 a} \leq c), (a = 0) &and (b \neq 0), &and (a = 0, b = 0, c = 0), (0 < a) &and (c \leq \frac{b^{2}}{4 a}))$

(8.3.4)

More examples on generating extended Tarski Formula.

>	f := &E([y]), y^2+x^2=2; out := QuantifierElimination(f, output=rootof);

$f ≔ &E ([y]), x^{2} + y^{2} = 2$

$out ≔ (- \sqrt{2} \leq x) &and (x \leq \sqrt{2})$

(8.3.5)

>	f := &E([y]), y^4+x^4=2; out := QuantifierElimination(f, output=rootof);

$f ≔ &E ([y]), x^{4} + y^{4} = 2$

$out ≔ (RootOf ({_Z}^{4} - 2, index = {real}_{1}) \leq x) &and (x \leq RootOf ({_Z}^{4} - 2, index = {real}_{2}))$

(8.3.6)

>	evalf(op(1, out)); evalf(op(2, out));

$−1.189207115 \leq x$

$x \leq 1.189207115$

(8.3.7)

An application for computing control Lyapunov function.

>	f1 := -x_1+u: f2 := -x_1-x_2^3: V := a_1x_1^2+a_2x_2^2; Vt := diff(V, x_1)f1 + diff(V, x_2)f2;

$V ≔ a_1 {x_1}^{2} + a_2 {x_2}^{2}$

$Vt ≔ 2 a_1 x_1 (- x_1 + u) + 2 a_2 x_2 (- {x_2}^{3} - x_1)$

(8.3.8)

>	QuantifierElimination( &A([x_1,x_2]), &E([u]), (x_1<>0) &or (x_2<>0) &implies ((V>0) &and (Vt<0)) );

$(0 < a_2) &and (0 < a_1)$

(8.3.9)

>	QuantifierElimination( &A([x_1, x_2]), &E([u]), (u=b_1x_1+b_2x_2) &and (a_1>0) &and (a_2>0) &and ((x_1<>0) &or (x_2<>0) &implies ((Vt<0))) );

$&and (0 < a_2, 0 < a_1, b_2 a_1 - a_2 = 0, b_1 < 1)$

(8.3.10)

Simplification of the output.

Without simplification:

>	ff := &E([i,j]), (0 <= i) &and (i <= n) &and (0 <= j) &and (j <= n) &and (t = n - j) &and (p = i + j); R := PolynomialRing([i,j,p,t,n]); sols := QuantifierElimination(ff, R, output=rootof, simplification=false);

$ff ≔ &E ([i, j]), (((((0 \leq i) &and (i \leq n)) &and (0 \leq j)) &and (j \leq n)) &and (t = n - j)) &and (p = i + j)$

$R ≔ polynomial_ring$

$sols ≔ &or (&and (n = 0, t = n, p = 0), &and (0 < n, t = 0, p = n), &and (0 < n, t = 0, n < p, p < 2 n), &and (0 < n, t = 0, p = 2 n), &and (0 < n, 0 < t, t < n, p = - t + n), &and (0 < n, 0 < t, t < n, - t + n < p, p < - t + 2 n), &and (0 < n, 0 < t, t < n, p = - t + 2 n), &and (0 < n, t = n, p = 0), &and (0 < n, t = n, 0 < p, p < n), &and (0 < n, t = n, p = n))$

(8.3.11)

With simplification:

>	sols := QuantifierElimination(ff, R, output=rootof, simplification=L4);

$sols ≔ &and (0 \leq n, 0 \leq t, t \leq n, n - t \leq p, p \leq - t + 2 n)$

(8.3.12)

9. Computing limits: from rational functions to topological closures

The AlgebraicGeometryTools submodule provides facilities for studying algebraic curves, surfaces and algebraic sets

of higher dimension. The commands currently available mainly focus on computing limits of "family of objects"

such as limits of family of secants in the case of tangent cone computation.

Different flavors of limit points

The command LimitPoints is part of AlgebraicGeometryTools package. The purpose of this command is to compute the limit points of a constructible set rc of dimension one. These limit points are corresponding to the values of the free variable of regular chain rc vanishing the product of the initials of the polynomials of rc.The command LimitPoints comes in two different types: 1) limit points w.r.t Zariski topology, 2) limit points w.r.t Euclidean topology. If one wants to compute the limit points of a constructible set with respect to Zariski topology, then one needs to use LimitPoints as follows:

>	restart:with(RegularChains): with(ChainTools): with(AlgebraicGeometryTools): R := PolynomialRing([x, y, z]): rc := Chain([y^5-z^4, x*z-y^2], Empty(R), R): lm := LimitPoints(rc, R): Display(lm, R);

$[\{\begin{array}{c} x = 0 \\ y = 0 \\ z = 0 \end{array}]$

(9.1.1)

Furthermore, the LimitPoints command uses an option called coefficient in order to indicate what the coefficient ring for the computations is. For the current implementation, coefficient can accept only complex and real corresponding to the limit points w.r.t Zariski and Euclidean topology, respectively. The option coefficient is set to complex by default. Therefore the following is equivalent to the computations above for lm.

>	rc := Chain([y^5-z^4, x*z-y^2], Empty(R), R): lm := LimitPoints(rc, R, coefficient = complex): Display(lm, R);

$[\{\begin{array}{c} x = 0 \\ y = 0 \\ z = 0 \end{array}]$

(9.1.2)

When the coefficient is set as complex, then the output of LimitPoints is called the complex limit points or for short only limit points of the regular chain rc. When the coefficient is set as real then this output is called the real limit points of the regular chain rc. The following shows how to compute the real limit points of the regular chain rc.

>	rc := Chain([y^5-z^4, x*z-y^2], Empty(R), R): lm := LimitPoints(rc, R, coefficient = real): Display(lm, R);

$[\{\begin{array}{c} x = 0 \\ y = 0 \\ z = 0 \end{array}]$

(9.1.3)

The command LimitPoints also has another option in order to indicate for which value of the free variable of regular chain rc one wants to compute the limit points of rc. In order to indicate this option, one new argument should be passed to LimitPoints which is a list of a single univariate polynomial in free variable of rc whose zeros are the corresponding values one wants to compute the limit points of rc at.

The following line shows how LimitPoints returns all the non-trivial limit points of rc.

>	rc:= Chain([y^5-z^4(z+1)^5,xz*(z+1)^2-y^2],Empty(R), R): lm := LimitPoints(rc, R); Display(lm, R);

$lm ≔ [regular_chain, regular_chain, regular_chain]$

$[\{\begin{array}{c} x = 0 \\ y = 0 \\ z = 0 \end{array}, \{\begin{array}{c} x + 1 = 0 \\ y = 0 \\ z + 1 = 0 \end{array}, \{\begin{array}{c} x^{4} - x^{3} + x^{2} - x + 1 = 0 \\ y = 0 \\ z + 1 = 0 \end{array}]$

(9.1.4)

While the following computes the limit points only related to some values of the free variable of the regular chain rc.

>	lm := LimitPoints(rc, R, [z]):Display(lm, R);

$[\{\begin{array}{c} x = 0 \\ y = 0 \\ z = 0 \end{array}]$

(9.1.5)

Real limits vs complex limits

While for some cases real limit points might be equal to the complex limit points, this is not the case all the time. For the following example the real limit points and complex ones are all equal.

>	R := PolynomialRing([y, x, z]): rc := Chain([z^5+x^4-2x^3+x^2, yz^4+x^3-x^2], Empty(R), R): lm := LimitPoints(rc, R, coefficient = complex): Display(lm, R);

$[\{\begin{array}{c} y = 0 \\ x = 0 \\ z = 0 \end{array}]$

(9.2.1)

>	m := LimitPoints(rc, R, coefficient = real): Display(lm, R);

$[\{\begin{array}{c} y = 0 \\ x = 0 \\ z = 0 \end{array}]$

(9.2.2)

But the next example shows a case for which the real limit points are different from the complex ones, even if the complex limit points all have rational coefficients.

>	R := PolynomialRing([x, y, z]): rc:=Chain([y^(3)-2* y^(3)+y^(2)+z^(5), z^(4)* x+y^(3)-y^(2)], Empty(R),R): lm := LimitPoints(rc, R, coefficient = complex): Display(lm, R);

$[\{\begin{array}{c} x = 0 \\ y = 0 \\ z = 0 \end{array}, \{\begin{array}{c} x = 0 \\ y - 1 = 0 \\ z = 0 \end{array}]$

(9.2.3)

>	lm := LimitPoints(rc, R, coefficient = real): Display(lm, R);

$[\{\begin{array}{c} x = 0 \\ y - 1 = 0 \\ z = 0 \end{array}]$

(9.2.4)

Different output formats of LimitPoints

The command LimitPoints has two different output formats, which can be controlled by the output option. A user can indicate output option to be set as chain, which is the default of LimitPoints command, or rootof, each of which gives a different representation of the limit points of rc.

>	R := PolynomialRing([x, y, z]): rc := Chain([y^4-(z^2-2)^4, (z^2-2)*x-y^2], Empty(R), R): lm := LimitPoints(rc, R, [z^2-2], output = rootof, coefficient = real): Display(lm, R);

$[[z = RootOf ({_Z}^{2} - 2), y = 0, x = 0]]$

(9.3.1)

>	lm := LimitPoints(rc, R, [z^2-2], output = chain, coefficient = real): Display(lm, R);

$[\{\begin{array}{c} x = 0 \\ y = 0 \\ z^{2} - 2 = 0 \end{array}]$

(9.3.2)

Computing the branches of the one-dimensional constructible sets

The command RegularChainBranches is part of AlgebraicGeometryTools package. The goal of this command is to compute a parametrization of the branches of a constructible set of dimension one represented by a regular chain. These parametrizations are actually the Puiseux parametrization of a constructible set when the free variable approaches zero. The following example show how to use this command without any extra option:

>	with(AlgebraicGeometryTools): R := PolynomialRing([x, y, z]): rc := Chain([-z^2+y, x*z-y^2], Empty(R), R): br := RegularChainBranches(rc, R, [z]);

$br ≔ [[z =_T, y = {_T}^{2}, x = {_T}^{3}]]$

(9.4.1)

>	rc := Chain([y^2z+y+1, (z+2)zx^2+(y+1)(x+1)], Empty(R), R): RegularChainBranches(rc, R, [z]); RegularChainBranches(rc, R, [z+2]);

$[[z =_T0, y = -_T0 - 1, x = - \frac{1}{432} {_T0}^{5} + \frac{11}{432} {_T0}^{4} + \frac{5}{432} {_T0}^{3} - \frac{5}{216} {_T0}^{2} + \frac{1}{12}_T0 - \frac{1}{2}], [z =_T0, y = -_T0 - 1, x = \frac{(_T0 - 2) ({_T0}^{2} + 4) ({_T0}^{2} - 9_T0 - 54)}{432}]]$

$[[z =_T1 - 2, y = \frac{1}{216} {_T1}^{4} + \frac{11}{216} {_T1}^{3} - \frac{5}{216} {_T1}^{2} - \frac{1}{12}_T1 - \frac{1}{2}, x = \frac{9}{8} {_T1}^{3} + 2 {_T1}^{2} + 4_T1 - 1], [z =_T1 - 2, y = - \frac{1}{216} {_T1}^{4} - \frac{11}{216} {_T1}^{3} + \frac{4}{27} {_T1}^{2} + \frac{1}{3}_T1 + 1, x = \frac{3}{8} {_T1}^{3} + \frac{1}{2} {_T1}^{2} +_T1 - 1]]$

(9.4.2)

Similar to the LimitPoints command, RegularChainBranches also works in two different modes: 1) complex, and 2) real. To switch between two different modes, user can use coefficient option which by default is set to be complex. The option coefficient indicates in which coefficient ring the corresponding parametrization of the branches will be computed. If it is set as real then it means that all the coefficients of the parametrizations are real numbers, complex, otherwise.

>	rc:=Chain([y^(3)-2* y^(3)+y^(2)+z^(5), z^(4)* x+y^(3)-y^(2)], Empty(R),R): br := RegularChainBranches(rc, R, [z], coefficient = complex);

$br ≔ [[z = {_T2}^{2}, y = - \frac{{_T2}^{5} ({_T2}^{5} + 2 RootOf ({_Z}^{2} + 1))}{2}, x = \frac{{_T2}^{2} ({_T2}^{20} + 6 {_T2}^{15} RootOf ({_Z}^{2} + 1) - 10 {_T2}^{10} - 8)}{8}], [z = {_T2}^{2}, y = - \frac{{_T2}^{5} ({_T2}^{5} - 2 RootOf ({_Z}^{2} + 1))}{2}, x = \frac{{_T2}^{2} ({_T2}^{20} - 6 {_T2}^{15} RootOf ({_Z}^{2} + 1) - 10 {_T2}^{10} - 8)}{8}], [z =_T2, y = {_T2}^{5} + 1, x = - {_T2}^{11} - 2 {_T2}^{6} -_T2]]$

(9.4.3)

>	br := RegularChainBranches(rc, R, [z], coefficient = real);

$br ≔ [[z =_T3, y = {_T3}^{5} + 1, x = - {_T3}^{11} - 2 {_T3}^{6} -_T3]]$

(9.4.4)

In the parametrization of the branches of the constructible set represented as a regular chain, a new (parametric) variable is introduced. This variable by default is chosen to be _T, possibly followed by a number, but also can be controlled by user. The following examples shows how this can be done:

>	rc := Chain([-z^2+y, x*z-y^2], Empty(R), R): br := RegularChainBranches(rc, R, [z^2+1], coefficient = complex);

$br ≔ [[z =_T4 + RootOf ({_Z}^{2} + 1), y = 2_T4 RootOf ({_Z}^{2} + 1) + {_T4}^{2} - 1, x = (- {_T4}^{8} + 3 {_T4}^{6} + 3 {_T4}^{2} - 1) RootOf ({_Z}^{2} + 1) + 3 {_T4}^{7} - {_T4}^{5} + {_T4}^{3} - 3_T4]]$

(9.4.5)

>	rc := Chain([-z^2+y, x*z-y^2], Empty(R), R): br := RegularChainBranches(rc, R, [z^2+1], coefficient = complex, W);

$br ≔ [[z = W + RootOf ({_Z}^{2} + 1), y = 2 W RootOf ({_Z}^{2} + 1) + W^{2} - 1, x = (- W^{8} + 3 W^{6} + 3 W^{2} - 1) RootOf ({_Z}^{2} + 1) + 3 W^{7} - W^{5} + W^{3} - 3 W]]$

(9.4.6)

Also it is possible to determine how many terms can be computed in the parametrization of the branches of the regular chain rc. This feature can be used as follows:

First create a list of accuracies whose elements is equal to the number of polynomials in rc.

>	L:= [1,8];

$L ≔ [1, 8]$

(9.4.7)

Then, we have:

>	R := PolynomialRing([x, y, z]): rc := Chain([y^2z+y+1, (z+2)zx^2+(y+1)(x+1)], Empty(R), R): RegularChainBranches(rc, R, [z],accuracies = L);

$[[z =_T5, y = -_T5 - 1, x = \frac{(_T5 - 2) ({_T5}^{4} + 2 {_T5}^{3} + 4 {_T5}^{2} + 8_T5 + 16) (143 {_T5}^{8} - 396 {_T5}^{7} + 1134 {_T5}^{6} - 3402 {_T5}^{5} + 10935 {_T5}^{4} - 39366 {_T5}^{3} + 177147 {_T5}^{2} - 1594323_T5 - 9565938) ({_T5}^{4} - 2 {_T5}^{3} + 4 {_T5}^{2} - 8_T5 + 16)}{4897760256}], [z =_T5, y = -_T5 - 1, x = - \frac{(_T5 - 2) ({_T5}^{4} + 2 {_T5}^{3} + 4 {_T5}^{2} + 8_T5 + 16) (143 {_T5}^{8} - 396 {_T5}^{7} + 1134 {_T5}^{6} - 3402 {_T5}^{5} + 10935 {_T5}^{4} - 39366 {_T5}^{3} + 177147 {_T5}^{2} - 1594323_T5 - 4782969) ({_T5}^{4} - 2 {_T5}^{3} + 4 {_T5}^{2} - 8_T5 + 16)}{4897760256}]]$

(9.4.8)

If now we change the accuracies, we have the following:

>	L:= [1,4]: RegularChainBranches(rc, R, [z],accuracies = L);

$[[z =_T6, y = -_T6 - 1, x = - \frac{5}{139968} {_T6}^{9} + \frac{7}{34992} {_T6}^{8} - \frac{137}{139968} {_T6}^{7} + \frac{1003}{139968} {_T6}^{6} + \frac{181}{139968} {_T6}^{5} - \frac{181}{69984} {_T6}^{4} + \frac{29}{3888} {_T6}^{3} - \frac{5}{216} {_T6}^{2} + \frac{1}{12}_T6 - \frac{1}{2}], [z =_T6, y = -_T6 - 1, x = \frac{5}{139968} {_T6}^{9} - \frac{7}{34992} {_T6}^{8} + \frac{137}{139968} {_T6}^{7} - \frac{1003}{139968} {_T6}^{6} - \frac{37}{2187} {_T6}^{5} + \frac{74}{2187} {_T6}^{4} - \frac{17}{243} {_T6}^{3} + \frac{4}{27} {_T6}^{2} - \frac{1}{3}_T6 + 1]]$

(9.4.9)

It is worth mentioning that when the list L is not presented, it means a default value is computed for the list L. This value is called list of accuracies and it meets some requirements. These requirements can be found at Computing the Limit Points of the Quasi-component of a Regular Chain in Dimension One.

Computing limits of fractional multivariate polynomials (when origin is the singular point of denominator)

RationalFunctionLimit is a command of AlgebraicGeometryTools that is used to compute the limit of the fraction of multivariate polynomials at a point where this point is an isolated zero of the denominator. RationalFunctionLimit accepts two arguments: 1) the fraction of multivariate polynomials q, 2) a list of all the variables equal to some values, indicating the point we aim at computing the limit of q at. The following examples demonstrate how this command works:

>	restart: with(RegularChains): with(AlgebraicGeometryTools): RationalFunctionLimit(x^2yz^2/(x^4+y^4+z^4), [x = 0, y = 0, z = 0]);

$0$

(9.5.1)

>	RationalFunctionLimit((w*z+x^2+y^2)/(w^2+x^2+y^2+z^2), [x = 0, y = 0, z = 0, w = 0]);

$undefined$

(9.5.2)

>	RationalFunctionLimit(x^6/(w^6+l^2+t^2+x^2+y^2+z^2), [x = 0, y = 0, z = 0, w = 0, t = 0, l = 0]);

$0$

(9.5.3)

>	RationalFunctionLimit(xyz*w/(w^4+x^4+y^4+z^4), [x = 0, y = 0, z = 0, w = 0]);

$undefined$

(9.5.4)

>	RationalFunctionLimit((xy+xz-y*z)/(x^2+y^2+z^2), [x = 0, y = 0, z = 0]);

$undefined$

(9.5.5)

Computing tangent cone of algebraic curves

The command TangentCone is integrated into AlgebraicGeometryTools package. The goal of this command is to compute the tangent cone of a polynomial system of dimension one at some points on the curve. These points can do not necessarily have rational coefficients. The following shows how to use this command.

First introduce the polynomial ring:

>	R := PolynomialRing([x, y, z]):

Second, introduce the points using a zero-dimensional regular chain:

>	with(ChainTools):

>	rc := Chain([z-1, y, x], Empty(R), R):

Then, introduce the polynomial system of dimension one, and compute the tangent cone of this system at the point represented by rc:

>	Gs := [x^2+y^2+z^2-1, x^2-y^2-z*(z-1)]: tc := TangentCone(rc, Gs, R);

$tc ≔ [[[_z - 1, 3 {_x}^{2} - {_y}^{2}], regular_chain]]$

(9.6.1)

There are two different representations for the output of TangentCone. Since the tangent cone of the polynomial systems of dimension one are lines, those lines can be represented both with their equations or slopes. The default is set as equations. Each format of the output can be used as follows:

>	tc := TangentCone(rc, Gs, R, output = equations);

$tc ≔ [[[_z - 1, 3 {_x}^{2} - {_y}^{2}], regular_chain]]$

(9.6.2)

>	tc := TangentCone(rc, Gs, R, output = slopes);

$tc ≔ [[[z, y - 1, 3 x^{2} - 1], regular_chain], [[z, y^{2} - 3, x - 1], regular_chain]]$

(9.6.3)

In order to indicate the equations or the slopes of the tangent cone, some new variables are used in order to distinguish the coordinate system of the equations and slopes from the original coordinate system. The default variables used for equation mode are the original names of variables where underscore is added as prefix to the variables. For slopes percentage added as prefix to original variables is used as default variables. However, the default for these variables can be controlled by user as follows:

>	tc := TangentCone(rc, Gs, R, [t, s, w], output = slopes);

$tc ≔ [[[w, s - 1, 3 t^{2} - 1], regular_chain], [[w, s^{2} - 3, t - 1], regular_chain]]$

(9.6.4)

Computing the tangent plane of the multivariate polynomials

TangentPlane is another command of AlgebraicGeometryTools package. Its goal is to compute the tangent plane of an algebraic surface represented by a polynomial at some point on the surface.

>	rc := Chain([z, y-1, x], Empty(R), R): tp := TangentPlane(rc, x*y, R);

$tp ≔ [[_x, rc]]$

(9.7.1)

Computing the intersection multiplicity of a polynomial system of dimension zero

The command TriangularizeWithMultiplicity is part of AlgebraicGeometryTools package. This command targets the intersection multiplicity computations of a polynomial system of dimension zero.

This command first attempts to solve this system using the Triangularize command of the RegularChains Library without taking multiplicity into consideration. Then for each point represented by the zero-dimensional regular chains in the triangular decomposition of the input system, it computes the intersection multiplicity of the original system at this related point.

>	R := PolynomialRing([z, y, x]): F := [x^2+y+z-1, y^2+x+z-1, z^2+x+y-1]: dec := TriangularizeWithMultiplicity(F, R): Display(dec, R);

$[[1, \{\begin{array}{c} z - x = 0 \\ y - x = 0 \\ x^{2} + 2 x - 1 = 0 \end{array}], [2, \{\begin{array}{c} z = 0 \\ y = 0 \\ x - 1 = 0 \end{array}], [2, \{\begin{array}{c} z = 0 \\ y - 1 = 0 \\ x = 0 \end{array}], [2, \{\begin{array}{c} z - 1 = 0 \\ y = 0 \\ x = 0 \end{array}]]$

(9.8.1)

Return to Index for Example Worksheets

Download Help Document

Maple

Maple Add-Ons

MapleSim

MapleSim Add-Ons

Systems Engineering

Consulting Services

Maple T.A. and Möbius

Education

Industries

Automotive and Aerospace

Robotics

Machine Design & Industrial Automation

Other

Application Areas

Product Pricing

Purchasing

Institutional Student Licensing

Maplesoft Elite Maintenance (EMP)

Support

Product Training

Online Product Help

Webinars & Events

Publications

Content Hubs

Examples & Applications

Community

About Maplesoft

Media Center

User Community

Contact

Online Help

All Products Maple MapleSim

Maple

Powerful math software that is easy to use

Maple Add-Ons

MapleSim

Advanced System Level Modeling

MapleSim Add-Ons

Systems Engineering

Consulting Services

Maple T.A. and Möbius

Education

Industries

Automotive and Aerospace

Robotics

Machine Design & Industrial Automation

Other

Application Areas

Product Pricing

Purchasing

Institutional Student Licensing

Maplesoft Elite Maintenance (EMP)

Support

Product Training

Online Product Help

Webinars & Events

Publications

Content Hubs

Examples & Applications

Community

About Maplesoft

Media Center

User Community

Contact

Online Help

All Products Maple MapleSim