Function Repository Resource:

TriangularSets

Source Notebook

Create a triangular set decomposition for a given list of polynomials and variables

Contributed by: Daniel Lichtblau

ResourceFunction["TriangularSets"][polys,vars]

computes a triangular decomposition of polynomials polys with respect to the ordered list of variables vars.

Details and Options

A triangular decomposition is a type of characteristic set, as defined by Ritt and Wu.
A set of polynomials is triangular with respect to a set of variables if the last polynomial is univariate in the last variable, the penultimate one has only the last two variables, etc.
The solution set for a polynomial ideal with finitely many solutions can be represented by a list of triangular sets.
Every solution for the polynomial system is also a solution of at least one triangular set. Every solution of a triangular set is either a solution of the original system, or else causes a leading term of a set member to vanish (where the nth element is regarded as a polynomial in the nth variable, with coefficients in lesser variables).
A Gröbner basis in lexicographic order is often a triangular set. Specifically, it is a triangular set whenever the ideal is radical (has no multiple solutions), and in general position with respect to the variables (has no coordinates taking the same value in more than one solution).
ResourceFunction["TriangularSets"] uses iterated pseudoremainders, while also removing polynomial divisors common to the multiplier and remainder coefficients.
ResourceFunction["TriangularSets"] replaces powers of a polynomial with the polynomial itself. In so doing, it destroys multiplicity in a solution set.
ResourceFunction["TriangularSets"] takes a Modulus option to compute triangular sets over a prime field.

Examples

Basic Examples (2) 

Triangularize a system of three polynomials in three variables:

In[1]:=
polys = {x^2 + y^2 - 1, x^3 + (2 + z)*x*y + y^3 - 1, z^2 - 2};
vars = {x, y, z};
ts = ResourceFunction["TriangularSets"][polys, vars]
Out[1]=

Check that all solutions to this system (as found by NSolve) are solutions of at least one of the triangular sets:

In[2]:=
sols = NSolve[polys, WorkingPrecision -> 100];
len = Length[sols];
zposns = Position[Chop[N[ts /. sols], 10^(-5)], {0 ..}][[All, 1]];
Union[zposns] == Range[len]
Out[2]=

Scope (3) 

TriangularSets can handle more than three variables, as illustrated on a standard benchmark known as the Trinks system:

In[3]:=
polys = {45*p + 35*s - 165*b - 36, 35*p + 40*z + 25*t - 27*s, 15*w + 25*p*s + 30*z - 18*t - 165*b^2, -9*w + 15*p*t + 20*z*s, w*p + 2*z*t - 11*b^3, 99*w - 11*s*b + 3*b^2};
vars = {w, t, z, b, s, p};
ts = ResourceFunction["TriangularSets"][polys, vars]
Out[3]=

Check the size of the result in terms of number of sets and leaf count:

In[4]:=
{Length[ts], LeafCount[ts]}
Out[4]=

Show that the solutions are all solutions of at least one triangular set:

In[5]:=
sols = NSolve[polys == 0, vars, WorkingPrecision -> 100];
len = Length[sols];
zposns = Position[Chop[N[ts /. sols], 10^(-5)], {0 ..}][[All, 1]];
{len, Union[zposns] == Range[len]}
Out[5]=

Find a triangular set decomposition for a perturbation of another benchmark known as the Caprasse system:

In[6]:=
polys = {-2*x + 2*t*x*y - (100001*z)/100000 + y^2*z, 2000001/1000000 + 4*x^2 - 10*t*y + 4*t*x^2*y - 10*y^2 + 2*t*y^3 + 4*x*z - x^3*z + 4*x*y^2*z, -x + t^2*x - 2*z + 2*t*y*z, 2 - 10*t^2 - 10*t*y + 2*t^3*y + 4*x*z + 4*t^2*x*z + 4*z^2 + 4*t*y*z^2 - x*z^3};
vars = {x, y, z, t};
Timing[ts = ResourceFunction["TriangularSets"][polys, vars];]
Out[6]=

Check the results:

In[7]:=
sols = Union[NSolve[polys == 0, vars, WorkingPrecision -> 100]];
len = Length[sols];
zposns = Position[Chop[N[ts /. sols], 10^(-5)], {0 ..}];
{len, Union[zposns[[All, 1]]] == Range[len]}
Out[7]=

Decompose the Cassou–Nouges benchmark system into triangular sets:

In[8]:=
polys = {15*b^2*c*d^2 + 6*b^2*c^3 + 21*b^2*c^2*d - 144*b*c - 8*b*c^2*e - 28*b*c*d*e - 648*b*d + b*d^2*e + 9*b^2*d^3 - 120, 30*b^2*c^3*d - 32*c*d*e^2 - 720*b*c*d - 24*b*c^3*e - 432*b*c^2 + 576*c*e - 576*d*e + 16*b*c*d^2*e + 16*d^2*e^2 + 16*c^2*e^2 + 9*b^2*c^4 + 5184 + 39*b^2*c^2*d^2 + 18*b^2*c*d^3 - 432*b*d^2 + 24*b*d^3*e - 16*b*c^2*d*e - 240*c, 216*b*c*d - 162*b*d^2 - 81*b*c^2 + 5184 + 1008*c*e - 1008*d*e + 15*b*c^2*d*e - 15*b*c^3*e - 80*c*d*e^2 + 40*d^2*e^2 + 40*c^2*e^2, 261 + 4*b*c*d - 3*b*d^2 - 4*b*c^2 + 22*c*e - 22*d*e};
vars = {e, d, c, b};
Timing[ts = TriangularSetsWithSplitting[polys, vars];]
Out[8]=

Check the results:

In[9]:=
sols = Union[NSolve[polys == 0, vars, WorkingPrecision -> 100]];
len = Length[sols];
zposns = Position[Chop[N[ts /. sols], 10^(-5)], {0 ..}];
{len, Union[zposns[[All, 1]]] == Range[len]}
Out[9]=

Options (3) 

Modulus (3) 

TriangularSets can be given a prime modulus:

In[10]:=
polys = {-6*p^3 + 4*p^3*phi^3 + 15*phi^3*s^3*p - 3*phi^3*s^5 - 12*phi^3*s*p^2 - 3*phi^3*s*p + phi^3*s^3 + 7, -9*phi^3*s^2*p - 5*phi^3*s^2 - 6*s*p^3 + 3*phi^3*s^4 + 5*phi^3*p - 4, -12*s^2*p - 6*s^2 + 3*s^4 + 4*phi^2 + 3 + 12*p};
vars = {s, p, phi};
Timing[ts2mod113 = ResourceFunction["TriangularSets"][polys, vars, Modulus -> 113];]
Out[10]=

Obtain the modular solutions to the original system:

In[11]:=
modsols = Solve[polys == 0, vars, Modulus -> 113]
Out[11]=

The fifth triangular set contains all three solutions:

In[12]:=
Mod[ts2mod113 /. modsols, 113]
Out[12]=

Possible Issues (6) 

Speed can vary considerably depending on the order of the variables:

In[13]:=
(* Evaluate this cell to get the example input *) CloudGet["https://www.wolframcloud.com/obj/5bd411b8-e7ed-4d25-aa65-569e3a291e95"]
Out[13]=

The size of result can also vary considerably:

In[14]:=
Map[{Length[#], LeafCount[#]} &, {tsr, ts}]
Out[14]=

Even for a seemingly small system, this function can be slow and can return a very large result. Demonstrate this on the Katsura-5 benchmark system:

In[15]:=
katsura5polys = {x1^2 - x1 + 2*x2^2 + 2*x3^2 + 2*x4^2 + 2*x5^2 + 2*x6^2, 2*x1*x2 + 2*x2*x3 + 2*x3*x4 + 2*x4*x5 + 2*x5*x6 - x2, 2*x1*x3 + x2^2 + 2*x2*x4 + 2*x3*x5 + 2*x4*x6 - x3, 2*x1*x4 + 2*x2*x3 + 2*x2*x5 + 2*x3*x6 - x4, 2*x1*x5 + 2*x2*x4 + 2*x2*x6 + x3^2 - x5, x1 + 2*x2 + 2*x3 + 2*x4 + 2*x5 + 2*x6 - 1};
vars = {x1, x2, x3, x4, x5, x6};
Timing[ts = ResourceFunction["TriangularSets"][katsura5polys, Reverse@vars];]
Out[15]=

Check the number of sets and the leaf count of the result:

In[16]:=
{Length[ts], LeafCount[ts]}
Out[16]=

Check that all solutions of this system are contained in the solution sets of the triangular sets:

In[17]:=
sols = Union[
   NSolve[katsura5polys == 0, vars, WorkingPrecision -> 100]];
len = Length[sols];
chopvals = Chop[N[ts /. sols], 10^(-6)];
zposns = Position[chopvals, {0 ..}];
{len, Union[zposns[[All, 1]]] == Range[len]}
Out[17]=

Version History

  • 1.0.0 – 14 June 2021

Source Metadata

Related Resources

Author Notes

This is perhaps a relatively slow implementation. It uses a direct method involving successive pseudoremaindering sequences, with some amount of factorization at intermediate steps to split into subsystems.

License Information