Algebraic substitution with PolynomialReduce
Guido Walter Pettinari
in order to perform algebraic substitutions in expressions, I usually
rely on replacementFunction, that is a very nice wrapper to
PolynomialReduce made by Daniel Lichtblau and Andrzej Kozlowski (see
links below).
However, I realized that replacementFunction is not enough when
dealing with multiple substitutions. Take this expression (which
itself is a sub-expression of a more complicated one in my code):
terms = (2 kx1^2 kx2^2)/3 - (kx2^2 ky1^2)/3 + 2 kx1 kx2 ky1 ky2 - (
kx1^2 ky2^2)/3 + (2 ky1^2 ky2^2)/3 - (kx2^2 kz1^2)/3 - (
ky2^2 kz1^2)/3 + 2 kx1 kx2 kz1 kz2 + 2 ky1 ky2 kz1 kz2 - (
kx1^2 kz2^2)/3 - (ky1^2 kz2^2)/3 + (2 kz1^2 kz2^2)/3
If I define k1 = {kx1, ky1, kz1} and k2 = {kx2, ky2, kz2}, then the
above is equal to the following combination of scalar products:
-(1/3 k1.k1 k2.k2 - (k1.k2)^2)
In order to explicitly reproduce the equality (and thus simplify the
equation), I tried to apply replacementFunction several times, with no
success. So I turned to PolynomialReduce. However, the following:
polylist = { k1.k1 k2.k2, (k1.k2)^2 };
PolynomialReduce[ terms, polylist, Join[k1, k2] ] // Last
did not return a zero remainder.
By reading old messages in the mailing list, I see that this may be
due to the ordering of the variables in the third argument, i.e.
Join[ k1,k2 ]. Using GroebnerBasis, I managed to prove that 'terms'
can be expressed in terms of the polynomials in 'polylist'. In fact:
gb = GroebnerBasis[polylist, Join[k1, k2] ];
PolynomialReduce[terms, gb, Join[k1, k2] ] // Last
yields zero.
The question is: how can I obtain the coefficients of my expression
('term') with respect to my polynomials ('polylist'), given that I
already have the coefficients of 'terms' with respect to the Groebner
basis of 'polylist'?
Thank you very much for your consideration!
Best wishes,
Guido W. Pettinari
REPLACEMENT FUNCTION LINKS:
http://groups.google.com/group/comp.soft-sys.math.mathematica/browse_thread/thread/634a54e2e844837f/ddf3803cab4f3a5f?lnk=gst&q=replacementfunction#ddf3803cab4f3a5f
Daniel Lichtblau
Yeah. Okay. What you need is a conversion matrix. Not entirely trivial
to get those.
This is discussed in a manuscript long ago abandoned, located here:
http://library.wolfram.com/infocenter/MathSource/7523/
I'll provide the needed code below, and skip minor details, like why
anyone might expect it to do anything useful. If you know, you know,
else just regard it as you would any other jinn.
(* Groebner basis of a set of vectors wrt POT order *)
moduleGroebnerBasis[polys_, vars_, cvars_, opts___] := Module[
{newpols, rels, len = Length[cvars], gb, j, k, rul},
rels = Flatten[
Table[cvars[[j]]*cvars[[k]], {j, len}, {k, j, len}]];
newpols = Join[polys, rels];
gb = GroebnerBasis[newpols, Join[cvars, vars], opts];
rul = Map[(# :> {}) &, rels];
gb = Flatten[gb /. rul];
Collect[gb, cvars]]
(* Finds a Groebner basis wrt lex order, and a conversion matrix. *)
conversionMatrix[polys_, vars_] := Module[
{aa, coords, pmat, len = Length[polys], newpolys, mgb, gb, convmat,
fvar, rvars},
coords = Array[aa, len + 1];
fvar = First[coords];
rvars = Rest[coords];
pmat = Transpose[Join[{polys}, IdentityMatrix[len]]];
newpolys = pmat.coords;
mgb = moduleGroebnerBasis[newpolys, vars, coords];
gb = mgb /. Join[{fvar -> 1}, Thread[rvars -> 0]] /. 0 :> Sequence[];
convmat = Select[mgb, ! FreeQ[#, fvar] &] /. fvar -> 0;
{gb, convmat /. Thread[rvars -> Table[UnitVector[len, j], {j, len}]]}
]
Here is your example.
terms = (2 kx1^2 kx2^2)/3 - (kx2^2 ky1^2)/3 +
2 kx1 kx2 ky1 ky2 - (kx1^2 ky2^2)/3 + (2 ky1^2 ky2^2)/
3 - (kx2^2 kz1^2)/3 - (ky2^2 kz1^2)/3 + 2 kx1 kx2 kz1 kz2 +
2 ky1 ky2 kz1 kz2 - (kx1^2 kz2^2)/3 - (ky1^2 kz2^2)/
3 + (2 kz1^2 kz2^2)/3;
k1 = {kx1, ky1, kz1};
k2 = {kx2, ky2, kz2};
polylist = {k1.k1 k2.k2, (k1.k2)^2};
vars = Join[k1, k2];
gb1 = GroebnerBasis[polylist, vars];
In[82]:= {gb, convmat} = conversionMatrix[polylist, vars]
Out[82]= {{kx2^6*ky1^4 + 3*kx2^4*ky1^4*ky2^2 + 3*kx2^2*ky1^4*ky2^4 +
ky1^4*ky2^6 +
2*kx2^6*ky1^2*kz1^2 + 4*kx2^4*ky1^2*ky2^2*kz1^2 +
2*kx2^2*ky1^2*ky2^4*kz1^2 + kx2^6*kz1^4 + kx2^4*ky2^2*kz1^4 +
4*kx2^4*ky1^3*ky2*kz1*kz2 + 8*kx2^2*ky1^3*ky2^3*kz1*kz2 +
4*ky1^3*ky2^5*kz1*kz2 + 4*kx2^4*ky1*ky2*kz1^3*kz2 +
4*kx2^2*ky1*ky2^3*kz1^3*kz2 + kx2^4*ky1^4*kz2^2 +
2*kx2^2*ky1^4*ky2^2*kz2^2 + ky1^4*ky2^4*kz2^2 +
4*kx2^4*ky1^2*kz1^2*kz2^2 +
10*kx2^2*ky1^2*ky2^2*kz1^2*kz2^2 +
6*ky1^2*ky2^4*kz1^2*kz2^2 + 3*kx2^4*kz1^4*kz2^2 +
2*kx2^2*ky2^2*kz1^4*kz2^2 + 4*kx2^2*ky1^3*ky2*kz1*kz2^3 +
4*ky1^3*ky2^3*kz1*kz2^3 + 8*kx2^2*ky1*ky2*kz1^3*kz2^3 +
4*ky1*ky2^3*kz1^3*kz2^3 + 2*kx2^2*ky1^2*kz1^2*kz2^4 +
6*ky1^2*ky2^2*kz1^2*kz2^4 + 3*kx2^2*kz1^4*kz2^4 +
ky2^2*kz1^4*kz2^4 +
4*ky1*ky2*kz1^3*kz2^5 + kz1^4*kz2^6, kx2^6*ky1^3*ky2 +
3*kx2^4*ky1^3*ky2^3 + 3*kx2^2*ky1^3*ky2^5 + ky1^3*ky2^7 +
kx2^6*ky1*ky2*kz1^2 + 2*kx1*kx2^5*ky2^2*kz1^2 +
4*kx2^4*ky1*ky2^3*kz1^2 + 2*kx1*kx2^3*ky2^4*kz1^2 +
3*kx2^2*ky1*ky2^5*kz1^2 - kx2^6*ky1^2*kz1*kz2 +
kx2^4*ky1^2*ky2^2*kz1*kz2 + 5*kx2^2*ky1^2*ky2^4*kz1*kz2 +
3*ky1^2*ky2^6*kz1*kz2 - kx2^6*kz1^3*kz2 +
2*kx2^4*ky2^2*kz1^3*kz2 +
3*kx2^2*ky2^4*kz1^3*kz2 + kx2^4*ky1^3*ky2*kz2^2 +
2*kx2^2*ky1^3*ky2^3*kz2^2 + ky1^3*ky2^5*kz2^2 +
2*kx1*kx2^5*kz1^2*kz2^2 + 2*kx2^4*ky1*ky2*kz1^2*kz2^2 +
4*kx1*kx2^3*ky2^2*kz1^2*kz2^2 + 7*kx2^2*ky1*ky2^3*kz1^2*kz2^2 +
3*ky1*ky2^5*kz1^2*kz2^2 - kx2^4*ky1^2*kz1*kz2^3 +
2*kx2^2*ky1^2*ky2^2*kz1*kz2^3 + 3*ky1^2*ky2^4*kz1*kz2^3 +
5*kx2^2*ky2^2*kz1^3*kz2^3 + ky2^4*kz1^3*kz2^3 +
2*kx1*kx2^3*kz1^2*kz2^4 + kx2^2*ky1*ky2*kz1^2*kz2^4 +
3*ky1*ky2^3*kz1^2*kz2^4 + kx2^2*kz1^3*kz2^5 +
ky2^2*kz1^3*kz2^5,
(-kx2^4)*ky1^2 + 2*kx1*kx2^3*ky1*ky2 + 2*kx1*kx2*ky1*ky2^3 +
ky1^2*ky2^4 - kx2^4*kz1^2 - kx2^2*ky2^2*kz1^2 +
2*kx1*kx2^3*kz1*kz2 +
2*kx2^2*ky1*ky2*kz1*kz2 + 2*kx1*kx2*ky2^2*kz1*kz2 +
2*ky1*ky2^3*kz1*kz2 - kx2^2*ky1^2*kz2^2 +
2*kx1*kx2*ky1*ky2*kz2^2 +
ky1^2*ky2^2*kz2^2 + ky2^2*kz1^2*kz2^2 + 2*kx1*kx2*kz1*kz2^3 +
2*ky1*ky2*kz1*kz2^3 + kz1^2*kz2^4, (-kx2^6)*ky1^3 -
3*kx2^4*ky1^3*ky2^2 - 3*kx2^2*ky1^3*ky2^4 - ky1^3*ky2^6 -
kx2^6*ky1*kz1^2 - 2*kx1*kx2^5*ky2*kz1^2 -
4*kx2^4*ky1*ky2^2*kz1^2 -
2*kx1*kx2^3*ky2^3*kz1^2 - 3*kx2^2*ky1*ky2^4*kz1^2 +
2*kx1*kx2^5*ky1*kz1*kz2 - 4*kx2^2*ky1^2*ky2^3*kz1*kz2 -
2*kx1*kx2*ky1*ky2^4*kz1*kz2 - 4*ky1^2*ky2^5*kz1*kz2 -
2*kx2^4*ky2*kz1^3*kz2 - 2*kx2^2*ky2^3*kz1^3*kz2 -
kx2^4*ky1^3*kz2^2 -
2*kx2^2*ky1^3*ky2^2*kz2^2 - ky1^3*ky2^4*kz2^2 -
4*kx1*kx2^3*ky2*kz1^2*kz2^2 - 7*kx2^2*ky1*ky2^2*kz1^2*kz2^2 -
2*kx1*kx2*ky2^3*kz1^2*kz2^2 - 5*ky1*ky2^4*kz1^2*kz2^2 +
2*kx1*kx2^3*ky1*kz1*kz2^3 - 2*kx1*kx2*ky1*ky2^2*kz1*kz2^3 -
4*ky1^2*ky2^3*kz1*kz2^3 - 4*kx2^2*ky2*kz1^3*kz2^3 -
2*ky2^3*kz1^3*kz2^3 + kx2^2*ky1*kz1^2*kz2^4 -
2*kx1*kx2*ky2*kz1^2*kz2^4 - 5*ky1*ky2^2*kz1^2*kz2^4 -
2*ky2*kz1^3*kz2^5,
(-kx2^5)*ky1^3 - 4*kx2^3*ky1^3*ky2^2 +
2*kx1*kx2^2*ky1^2*ky2^3 -
3*kx2*ky1^3*ky2^4 + 2*kx1*ky1^2*ky2^5 - kx2^5*ky1*kz1^2 -
2*kx1*kx2^4*ky2*kz1^2 - 5*kx2^3*ky1*ky2^2*kz1^2 -
2*kx1*kx2^2*ky2^3*kz1^2 - 4*kx2*ky1*ky2^4*kz1^2 +
2*kx1*kx2^4*ky1*kz1*kz2 - 2*kx2^3*ky1^2*ky2*kz1*kz2 +
6*kx1*kx2^2*ky1*ky2^2*kz1*kz2 - 2*kx2*ky1^2*ky2^3*kz1*kz2 +
4*kx1*ky1*ky2^4*kz1*kz2 - 4*kx2^3*ky2*kz1^3*kz2 -
4*kx2*ky2^3*kz1^3*kz2 - kx2^3*ky1^3*kz2^2 -
3*kx2*ky1^3*ky2^2*kz2^2 +
2*kx1*ky1^2*ky2^3*kz2^2 - 3*kx2*ky1*ky2^2*kz1^2*kz2^2 +
2*kx1*ky2^3*kz1^2*kz2^2 + 2*kx1*kx2^2*ky1*kz1*kz2^3 -
2*kx2*ky1^2*ky2*kz1*kz2^3 + 4*kx1*ky1*ky2^2*kz1*kz2^3 -
4*kx2*ky2*kz1^3*kz2^3 + kx2*ky1*kz1^2*kz2^4 +
2*kx1*ky2*kz1^2*kz2^4,
kx1*kx2^4*ky1^2 + 2*kx2^3*ky1^3*ky2 + 2*kx2*ky1^3*ky2^3 -
kx1*ky1^2*ky2^4 + kx1*kx2^4*kz1^2 + 2*kx2^3*ky1*ky2*kz1^2 +
kx1*kx2^2*ky2^2*kz1^2 + 2*kx2*ky1*ky2^3*kz1^2 +
2*kx2^3*ky1^2*kz1*kz2 -
2*kx1*kx2^2*ky1*ky2*kz1*kz2 + 2*kx2*ky1^2*ky2^2*kz1*kz2 -
2*kx1*ky1*ky2^3*kz1*kz2 + 2*kx2^3*kz1^3*kz2 +
2*kx2*ky2^2*kz1^3*kz2 +
kx1*kx2^2*ky1^2*kz2^2 + 2*kx2*ky1^3*ky2*kz2^2 -
kx1*ky1^2*ky2^2*kz2^2 +
2*kx2*ky1*ky2*kz1^2*kz2^2 - kx1*ky2^2*kz1^2*kz2^2 +
2*kx2*ky1^2*kz1*kz2^3 - 2*kx1*ky1*ky2*kz1*kz2^3 +
2*kx2*kz1^3*kz2^3 -
kx1*kz1^2*kz2^4,
kx2^2*ky1^2 - 2*kx1*kx2*ky1*ky2 + kx1^2*ky2^2 +
kx2^2*kz1^2 + ky2^2*kz1^2 - 2*kx1*kx2*kz1*kz2 -
2*ky1*ky2*kz1*kz2 +
kx1^2*kz2^2 + ky1^2*kz2^2, kx1^2*kx2^2 + 2*kx1*kx2*ky1*ky2 +
ky1^2*ky2^2 + 2*kx1*kx2*kz1*kz2 + 2*ky1*ky2*kz1*kz2 +
kz1^2*kz2^2},
{{kx2^4*ky1^2 + 2*kx1*kx2^3*ky1*ky2 + 3*kx2^2*ky1^2*ky2^2 +
kx2^4*kz1^2 +
2*kx1*kx2^3*kz1*kz2 + 6*kx2^2*ky1*ky2*kz1*kz2 +
3*kx2^2*kz1^2*kz2^2,
(-kx2^4)*ky1^2 - 2*kx1*kx2^3*ky1*ky2 - 2*kx1*kx2*ky1*ky2^3 +
ky1^2*ky2^4 - kx2^4*kz1^2 - kx2^2*ky2^2*kz1^2 -
2*kx1*kx2^3*kz1*kz2 +
2*kx2^2*ky1*ky2*kz1*kz2 - 2*kx1*kx2*ky2^2*kz1*kz2 +
2*ky1*ky2^3*kz1*kz2 - kx2^2*ky1^2*kz2^2 -
2*kx1*kx2*ky1*ky2*kz2^2 +
ky1^2*ky2^2*kz2^2 + ky2^2*kz1^2*kz2^2 - 2*kx1*kx2*kz1*kz2^3 +
2*ky1*ky2*kz1*kz2^3 + kz1^2*kz2^4},
{kx2^4*ky1*ky2 + 2*kx1*kx2^3*ky2^2 + 3*kx2^2*ky1*ky2^3 -
kx2^4*kz1*kz2 +
3*kx2^2*ky2^2*kz1*kz2, (-kx2^4)*ky1*ky2 - 2*kx1*kx2^3*ky2^2 -
2*kx1*kx2*ky2^4 + ky1*ky2^5 + kx2^4*kz1*kz2 +
2*kx2^2*ky2^2*kz1*kz2 +
ky2^4*kz1*kz2 - kx2^2*ky1*ky2*kz2^2 - 2*kx1*kx2*ky2^2*kz2^2 +
ky1*ky2^3*kz2^2 + kx2^2*kz1*kz2^3 + ky2^2*kz1*kz2^3},
{-kx2^2,
kx2^2 + ky2^2 + kz2^2}, {(-kx2^4)*ky1 - 2*kx1*kx2^3*ky2 -
3*kx2^2*ky1*ky2^2 - 2*kx2^2*ky2*kz1*kz2,
kx2^4*ky1 + 2*kx1*kx2^3*ky2 +
2*kx1*kx2*ky2^3 - ky1*ky2^4 - 2*kx2^2*ky2*kz1*kz2 -
2*ky2^3*kz1*kz2 +
kx2^2*ky1*kz2^2 + 2*kx1*kx2*ky2*kz2^2 - ky1*ky2^2*kz2^2 -
2*ky2*kz1*kz2^3}, {(-kx2^3)*ky1 - 2*kx1*kx2^2*ky2 -
4*kx2*ky1*ky2^2 -
4*kx2*ky2*kz1*kz2,
kx2^3*ky1 + 2*kx1*kx2^2*ky2 + kx2*ky1*ky2^2 +
2*kx1*ky2^3 + kx2*ky1*kz2^2 + 2*kx1*ky2*kz2^2},
{kx1*kx2^2 + 2*kx2*ky1*ky2 + 2*kx2*kz1*kz2, (-kx1)*kx2^2 -
kx1*ky2^2 -
kx1*kz2^2}, {1, -1}, {0, 1}}}
Quick check:
In[84]:= gb === gb1
Out[84]= True
As you did, we'll reduce terms with respect to the Groebner basis of
polylist. We show the vector of reducing multipliers.
In[85]:= reduction = First[PolynomialReduce[terms, gb, vars]]
Out[85]= {0, 0, 0, 0, 0, 0, -(1/3), 2/3}
To get the reduction in terms of the original polylist, we compute
reduction.mgb.
In[76]:= origreduction = Expand[reduction.convmat]
Out[76]= {-(1/3), 1}
Let's check this.
In[77]:= Expand[origreduction.polylist - terms]
Out[77]= 0
So we have now successfully rewritten 'terms' in terms of polylist.
I will remark that for production code purposes, one might wish to alter
the above to use a faster term order. If that is a compelling need then
I might be able to dredge up some old code for the occasion.
Daniel Lichtblau
Wolfram Research
Guido Walter Pettinari
In the meanwhile I found a temporary solution, that is running PolynomialReduce on all permutations of 'vars' until a zero-remainder result was found, but yours is much faster and cleaner.
I am going to apply your solution to much tougher expression, and I'll let you know if I succeed.
Thank you again.
Cheers,
Guido
On 9 Mar 2011, at 16:26, Daniel Lichtblau wrote:
> Guido Walter Pettinari wrote:
>> Dear group,
>> in order to perform algebraic substitutions in expressions, I usually
>> rely on replacementFunction, that is a very nice wrapper to
>> PolynomialReduce made by Daniel Lichtblau and Andrzej Kozlowski (see
>> links below).
>> However, I realized that replacementFunction is not enough when
>> dealing with multiple substitutions. Take this expression (which
>> itself is a sub-expression of a more complicated one in my code):
>> terms == (2 kx1^2 kx2^2)/3 - (kx2^2 ky1^2)/3 + 2 kx1 kx2 ky1 ky2 - (
>> kx1^2 ky2^2)/3 + (2 ky1^2 ky2^2)/3 - (kx2^2 kz1^2)/3 - (
>> ky2^2 kz1^2)/3 + 2 kx1 kx2 kz1 kz2 + 2 ky1 ky2 kz1 kz2 - (
>> kx1^2 kz2^2)/3 - (ky1^2 kz2^2)/3 + (2 kz1^2 kz2^2)/3
>> If I define k1 == {kx1, ky1, kz1} and k2 == {kx2, ky2, kz2}, then the
>> above is equal to the following combination of scalar products:
>> -(1/3 k1.k1 k2.k2 - (k1.k2)^2)
>> In order to explicitly reproduce the equality (and thus simplify the
>> equation), I tried to apply replacementFunction several times, with no
>> success. So I turned to PolynomialReduce. However, the following:
>> polylist == { k1.k1 k2.k2, (k1.k2)^2 };
>> PolynomialReduce[ terms, polylist, Join[k1, k2] ] // Last
>> did not return a zero remainder.
>> By reading old messages in the mailing list, I see that this may be
>> due to the ordering of the variables in the third argument, i.e.
>> Join[ k1,k2 ]. Using GroebnerBasis, I managed to prove that 'terms'
>> can be expressed in terms of the polynomials in 'polylist'. In fact:
>> gb == GroebnerBasis[polylist, Join[k1, k2] ];
>> PolynomialReduce[terms, gb, Join[k1, k2] ] // Last
>> yields zero.
>> The question is: how can I obtain the coefficients of my expression
>> ('term') with respect to my polynomials ('polylist'), given that I
>> already have the coefficients of 'terms' with respect to the Groebner
>> basis of 'polylist'?
>> Thank you very much for your consideration!
>> Best wishes,
>> Guido W. Pettinari
>> REPLACEMENT FUNCTION LINKS:
>> http://groups.google.com/group/comp.soft-sys.math.mathematica/browse_thr=
ead/thread/634a54e2e844837f/ddf3803cab4f3a5f?lnk==gst&q==replacementfunctio=
n#ddf3803cab4f3a5f
>> http://groups.google.com/group/comp.soft-sys.math.mathematica/browse_thr=
ead/thread/1397ca25b770b9af/59f63501668f7fbe?lnk==gst&q==replacementfunctio=
n#59f63501668f7fbe
>
> Yeah. Okay. What you need is a conversion matrix. Not entirely trivial to=
get those.
>
> This is discussed in a manuscript long ago abandoned, located here:
>
> http://library.wolfram.com/infocenter/MathSource/7523/
>
> I'll provide the needed code below, and skip minor details, like why anyo=
ne might expect it to do anything useful. If you know, you know, else just =
regard it as you would any other jinn.
>
> (* Groebner basis of a set of vectors wrt POT order *)
>
> moduleGroebnerBasis[polys_, vars_, cvars_, opts___] :== Module[
> {newpols, rels, len == Length[cvars], gb, j, k, rul},
> rels == Flatten[
> Table[cvars[[j]]*cvars[[k]], {j, len}, {k, j, len}]];
> newpols == Join[polys, rels];
> gb == GroebnerBasis[newpols, Join[cvars, vars], opts];
> rul == Map[(# :> {}) &, rels];
> gb == Flatten[gb /. rul];
> Collect[gb, cvars]]
>
> (* Finds a Groebner basis wrt lex order, and a conversion matrix. *)
>
> conversionMatrix[polys_, vars_] :== Module[
> {aa, coords, pmat, len == Length[polys], newpolys, mgb, gb, convmat,
> fvar, rvars},
> coords == Array[aa, len + 1];
> fvar == First[coords];
> rvars == Rest[coords];
> pmat == Transpose[Join[{polys}, IdentityMatrix[len]]];
> newpolys == pmat.coords;
> mgb == moduleGroebnerBasis[newpolys, vars, coords];
> gb == mgb /. Join[{fvar -> 1}, Thread[rvars -> 0]] /. 0 :> Sequence[];
> convmat == Select[mgb, ! FreeQ[#, fvar] &] /. fvar -> 0;
> {gb, convmat /. Thread[rvars -> Table[UnitVector[len, j], {j, len}]]}
> ]
>
> Here is your example.
>
> terms == (2 kx1^2 kx2^2)/3 - (kx2^2 ky1^2)/3 +
> 2 kx1 kx2 ky1 ky2 - (kx1^2 ky2^2)/3 + (2 ky1^2 ky2^2)/
> 3 - (kx2^2 kz1^2)/3 - (ky2^2 kz1^2)/3 + 2 kx1 kx2 kz1 kz2 +
> 2 ky1 ky2 kz1 kz2 - (kx1^2 kz2^2)/3 - (ky1^2 kz2^2)/
> 3 + (2 kz1^2 kz2^2)/3;
> k1 == {kx1, ky1, kz1};
> k2 == {kx2, ky2, kz2};
> polylist == {k1.k1 k2.k2, (k1.k2)^2};
> vars == Join[k1, k2];
>
> gb1 == GroebnerBasis[polylist, vars];
>
> In[82]:== {gb, convmat} == conversionMatrix[polylist, vars]
>
> Out[82]== {{kx2^6*ky1^4 + 3*kx2^4*ky1^4*ky2^2 + 3*kx2^2*ky1^4*ky2^4 +
> In[84]:== gb ====== gb1
> Out[84]== True
>
> As you did, we'll reduce terms with respect to the Groebner basis of poly=
list. We show the vector of reducing multipliers.
>
> In[85]:== reduction == First[PolynomialReduce[terms, gb, vars]]
> Out[85]== {0, 0, 0, 0, 0, 0, -(1/3), 2/3}
>
> To get the reduction in terms of the original polylist, we compute reduct=
ion.mgb.
>
> In[76]:== origreduction == Expand[reduction.convmat]
> Out[76]== {-(1/3), 1}
>
> Let's check this.
>
> In[77]:== Expand[origreduction.polylist - terms]
> Out[77]== 0
>
> So we have now successfully rewritten 'terms' in terms of polylist.
>
> I will remark that for production code purposes, one might wish to alter =
the above to use a faster term order. If that is a compelling need then I m=