# Testing whether polynomial is in algebra of bunch of given polynomials

I know that mathematica can check whether one polynomial is in the ideal of a given collection of polynomials, but I was wondering whether it can check if a polynomial is in the algebra of a given collection of polynomials. By "algebra", I mean that you can only add, subtract, and multiply polynomials in the given collection, and multiply them by real numbers.

To be completely clear: we say a collection $$\Sigma$$ of polynomials is an algebra if: (1) $$\lambda f + \eta g \in \Sigma$$ for any $$f,g \in \Sigma, \lambda,\eta \in \mathbb{R}$$ and (2) $$f,g \in \Sigma$$ implies $$fg \in \Sigma$$. We say that $$P$$ is in the algebra of $$\{P_1,\dots,P_n\}$$ if $$P$$ is in the smallest algebra containing $$P_1,\dots,P_n$$.

I think mathematica defines an ideal as I defined an algebra above except (2) has the stronger form $$fg \in \Sigma$$ if $$f \in \mathbb{R}[x], g \in \Sigma$$, and mathematica can check if a given $$P$$ is in the smallest ideal containing given $$P_1,\dots,P_n$$.

It would be nice to have this test for any number of variables, i.e., in any $$\mathbb{R}[x_1,\dots,x_n]$$.

Fundamental example: take $$n \ge 1$$ and let $$P_1 = x_1+\dots+x_n, P_2 = x_1^2+\dots+x_n^2,\dots P_n = x_1^n+\dots+x_n^n$$. Then all $$n$$ of the following symmetric functions are in the algebra generated by $$P_1,\dots,P_n$$: $$x_1+\dots+x_n$$ $$x_1x_2+\dots+x_{n-1}x_n$$ $$x_1x_2x_3+\dots+x_{n-2}x_{n-1}x_n$$ $$\dots$$ $$x_1\dots x_n$$

I'd appreciate any help.

• Comments are not for extended discussion; this conversation has been moved to chat. – Mr.Wizard Jan 19 at 1:41

The trouble with my answering this sort of question is that I haven't kept up with symbolic algebra over last 25 years or so. This problem in particular seems like an area where there's been some progress over that time, as well as pitfalls that I don't quite remember.

In any case, here's an idea that's too long for a comment. The basic idea is to introduce variables $$p_i$$ for each polynomial $$P_i(x_1,\dots,x_n)$$ and $$y$$ for $$P(x_1,\dots,x_n)$$, and to use GroebnerBasis to eliminate the $$x_i$$ from $$\{ p_i - P_i(x_1,\dots,x_n),\ y - P(x_1,\dots,x_n)\}$$ and solve for $$y$$, which, if successful, will be in terms of the $$p_i$$. If unsuccessful, then I hope that means $$P$$ is not in the subalgebra generated by the $$P_i$$, but I'm not certain of that.

Example 1: p is in the subalgebra generated by p1, p2:

basis = {p1 -> x1^2 - x2^2, p2 -> x1^3 x2};
p = x1^4 - x1^5 x2 - 2 x1^2 x2^2 + x1^3 x2^3 + x2^4;
vars = Variables[{p, Last /@ basis}]; (* assumes all symbols are variables *)
ans = Last@
GroebnerBasis[Join[Subtract @@@ basis, {y - p}],
Join[{y}, First /@ basis], vars] // Solve[# == 0, y] & //
ConditionalExpression[y /. #,
Length@# == 1 &&
AllTrue[#, PolynomialQ[y /. #, First /@ basis] &]] &
(*  {p1^2 - p1 p2}  *)

First@ans == p /. basis // Simplify
(*  True  *)


Example 2: p is not in the subalgebra generated by p1, p2:

basis = {p1 -> x1^2 - x2^2, p2 -> x1^3 x2};
p = x1 + x1^4 - x1^5 x2 - 2 x1^2 x2^2 + x1^3 x2^3 + x2^4; (* x1 + previous p *)
vars = Variables[{p, Last /@ basis}];
Last@GroebnerBasis[Join[Subtract @@@ basis, {y - p}],
Join[{y}, First /@ basis], vars] // Solve[# == 0, y] & //
ConditionalExpression[y /. #,
Length@# == 1 &&
AllTrue[#, PolynomialQ[y /. #, First /@ basis] &]] &
(*  Undefined  *)


Example 3: The example in my comment featuring a constant is problematic. It can be handled by having a variable k that stands for a constant, and the eliminating it.

basis = {p1 -> x^2 + 3 k, p2 -> x^2};
p = 2 k;
vars = Variables[{p, Last /@ basis}];
Last@GroebnerBasis[Join[Subtract @@@ basis, {y - p}],
Join[{y}, First /@ basis], vars] // Solve[# == 0, y] & //
ConditionalExpression[y /. #,
Length@# == 1 &&
AllTrue[#, PolynomialQ[y /. #, First /@ basis] &]] &
(*  {(2 (p1 - p2))/3}  *)


It's possible using, say, CoefficientArrays to rewrite a polynomial with a constant term a replaced by a*k.

• Hi Michael, thanks a lot for your answer. I am really appreciative you took the time. However, it seems it is returning "undefined" for the "fundamental example" with $n=2$ (cf the edited question). – mathworker21 Jan 18 at 9:58
• (1) This looks to be correct, by the way (but see next note). Might be slow due to the term ordering, and a bit more work than necessary to determine solvability. – Daniel Lichtblau Jan 19 at 19:09
• (2) But i wonder if this GB can produce something that factors as (y-poly1(p1,p2,...)*(y-poly2(p1,p2,...)? In such a case p would not be in the subalgebra, but you would deliver an affirmative result. (If it is in the subalgebra, that polynomial in y needs to be of the form y+poly(p1,...)). – Daniel Lichtblau Jan 19 at 19:11

The way to do this using a Groebner basis is to force a way to write subalgebra elements only in terms of some generating set. Well, we are given the generating set. The idea then is

(i) Create a new variable for each generator. (ii) Order variables so these new ones are "smallest" in a well-founded term ordering. I won't go into details on what that means but suffice it to say its fairly standard in computational commutative algebra and related. (iii) Compute a Groebner basis with respect to this term order. (iv) To test whether a new polynomial is in the algebra, reduce it by the Groebner basis. It is in the subalgebra iff the reduced form is comprised of the new variables alone.

I'll illustrate with an example. First we need a bit of code to create the monomial orders. A type that is well suited for the task at hand makes terms in the new variables strictly less than terms that have the original variables in them, but in other respects is close to the "degree reverse lexicographic" order.

drlMatrix[n_] :=
Prepend[Table[-KroneckerDelta[j + k - (n + 1)], {j, n - 1}, {k, n}],
Table[1, {n}]]

elimMatrix[n1_, n2_] :=
Module[{row1, rest}, row1 = Join[Table[1, {n1}], Table[0, {n2 - n1}]];
rest = drlMatrix[n2];
rest = Drop[rest, {n1}];
Prepend[rest, row1]]


Now we create a simple example using three polynomials in x and y. I'll show the ordering matrix explicitly, as well as the augmented polynomial system.

polys = {3*x^3 + x*y^2, x + y^3, x^2*y + y};
vars = Variables[polys];
avars = Array[a, Length[polys]];
orderMatrix = elimMatrix[Length[vars], Length[vars] + Length[avars]]
genPolys = avars - polys

(* Out[46]= {{1, 1, 0, 0, 0}, {1, 1, 1, 1, 1}, {0, 0, 0, -1, 0}, {0,
0, -1, 0, 0}, {0, -1, 0, 0, 0}}

(* Out[47]= {-3 x^3 - x y^2 + a[1], -x - y^3 + a[2], -y - x^2 y + a[3]} *)


Now form the GB.

gb1 = GroebnerBasis[genPolys, Join[vars, avars], MonomialOrder -> orderMatrix];


We'll use a simple example where we know in advance that the polynomial is in the subalgebra (we create it as an explicit sum of products of generating elements).

PolynomialReduce[polys[[1]]*polys[[3]] + polys[[2]]^2, gb1,
Join[vars, avars], MonomialOrder -> orderMatrix][[2]]

(* Out[49]= a[2]^2 + a[1] a[3] *)


So we see it is indeed written appropriately in terms of the new variables.

Now a symmetric polynomial example (which we'll see is not in the subalgebra).

PolynomialReduce[x^3 + y^3, gb1, Join[vars, avars],
MonomialOrder -> orderMatrix][[2]]

(* Out[50]= -((87 x)/32) + (107 y)/288 - (55 a[1])/96 +
13/144 x a[1]^2 + 1/32 y a[1]^2 + (87 a[2])/32 - (
3225047 x a[1] a[2])/3052176 - (400235 y a[1] a[2])/508696 + (
493465 a[1]^2 a[2])/4578264 - (38311 x a[1]^3 a[2])/1017392 - (
148173 y a[1]^3 a[2])/2034784 + (6187465 x a[2]^2)/2034784 + (
32739307 y a[2]^2)/18313056 + (6579931 a[1] a[2]^2)/6104352 + (
109197 x a[1]^2 a[2]^2)/2034784 + (207697 y a[1]^2 a[2]^2)/763044 + (
38311 a[1]^3 a[2]^2)/1017392 - (163849 a[2]^3)/508696 - (
1813561 x a[1] a[2]^3)/3052176 + (206459 y a[1] a[2]^3)/2034784 - (
109197 a[1]^2 a[2]^3)/2034784 - (129837 x a[2]^4)/2034784 - (
109197 y a[2]^4)/2034784 - (107 a[3])/288 - 3/32 x a[1] a[3] -
107/288 y a[1] a[3] + (24499501 x a[2] a[3])/18313056 - (
18311051 y a[2] a[3])/6104352 + (1652099 a[1] a[2] a[3])/1017392 + (
830467 x a[1]^2 a[2] a[3])/3052176 + (
100029 y a[1]^2 a[2] a[3])/2034784 - (
50434999 a[2]^2 a[3])/18313056 - (
2918511 x a[1] a[2]^2 a[3])/2034784 + (
2115501 y a[1] a[2]^2 a[3])/2034784 - (
327377 a[1]^2 a[2]^2 a[3])/6104352 - (
1655769 x a[2]^3 a[3])/2034784 - (2292923 y a[2]^3 a[3])/6104352 - (
38311 a[1] a[2]^3 a[3])/1017392 + (109197 a[2]^4 a[3])/2034784 +
55/96 x a[3]^2 - 81/32 y a[3]^2 + (18311051 a[2] a[3]^2)/6104352 + (
159645 x a[1] a[2] a[3]^2)/2034784 + (
441266 y a[1] a[2] a[3]^2)/190761 - (
18631649 x a[2]^2 a[3]^2)/6104352 - (
982131 y a[2]^2 a[3]^2)/2034784 - (
114933 a[1] a[2]^2 a[3]^2)/508696 + (
2292923 a[2]^3 a[3]^2)/6104352 + (13 a[3]^3)/16 - (
7235327 x a[2] a[3]^3)/2034784 + (984057 y a[2] a[3]^3)/2034784 - (
344799 a[1] a[2] a[3]^3)/1017392 + (
1637527 a[2]^2 a[3]^3)/2034784 + (982131 a[2] a[3]^4)/2034784 *)


I will remark that without explicit examples it is not obvious how one might show this. Basically, I did some of your work, in making up (possibly irrelevant) examples to use for illustration. This is of course to explain why I and others were pressing for concrete examples.

Stated differently, would a response that just stated steps (i)-(iv) have been at all of use? I'm doubtful. To the point where I actually think the moderators might remove it for lacking adequate detail. Would you want a response like that? Would you write a response like that? Again, on these points I am doubtful.

• Did you forget to include the definition of elimMatrix[]? – Michael E2 Jan 19 at 1:38
• Thank you for taking the time to write this answer. I am still going through it. I ended up providing an example, so not sure what the fuss is still about. Regardless, I am genuinely appreciative of your answer and will get back to you on it shortly. – mathworker21 Jan 19 at 13:30
• @MichaelE2 Yeah. I left it in my office (literally). I'll find it and put it in later. – Daniel Lichtblau Jan 19 at 15:57
• @MichaelE2 Code is repaired. – Daniel Lichtblau Jan 19 at 19:12