% Setup better fontification for Mexin and Mexout
% Putting this into LaTeX-mode-hook or in the file variables does not work.
% (setq font-lock-syntactic-keywords '(t ("\\\\begin{\\(Mexin\\|Mexout\\)}\\(\n\\)" (2 '(7))) ("\\(\n\\)\\\\end{\\(Mexin\\|Mexout\\)}"   (1 '(7))) ("\\\\begin{verbatim}\\(\n\\)" (1 '(7))) ("\\(\n\\)\\\\end{verbatim}"   (1 '(7)))))


Comments from Haiman:

This problem looks more combinatorial than algebraic.

Remark: if you take the derivation operators $P_k := \sum x^k+1 \p_k$,
you can see them as vector fields; furthermore, they generate all the
S_n-invariant vector fields (those are of the form \sum f_i \p_i, with
the f_i polynomials) as a left module over symmetric functions.
Exercise: prove it! proof one: look geometrically at the chamber of
decreasing vectors in \R^n; proof two: use some elimination.

\begin{Mexin}
Pref::typeCheck(Always):
S:=examples::SymmetricFunctions():
export(experimental::Steenrod):
\end{Mexin}

Brutal conversion from the symmetric functions to Steenrod

\begin{Mexin}
toe :=
proc(p, e: DOM_DOMAIN, X=FAIL: Type::Union(Type::ListOf(Type::Indeterminate), DOM_FAIL))
    option escape;
begin
    p := S::p(p):
    if X=FAIL then
	X := e::DXVariables;
    end_if;
    _plus(_mult(_plus(e((1+q*t.x)*x^k) $ x in X)
		$k in op(term,2)) * op(term,1)
	  $ term in poly2list(p));
end_proc:
X:=[x1,x2,x3]:
Y:=[y1,y2]:
e:=OREtx(X.Y):

toe(S::p([1]), e,X);
\end{Mexin} %$

Here is the usual relations between the symmetric functions:
\begin{Mexin}
X:=[x1,x2,x3]: Y:=[y1,y2]: expand(_plus((-1)^i*S::e([i])(Y)*S::h([k-i])(X.Y) $i=0..k)-S::h([k])(X)) $k=1..5;
\end{Mexin}

Now, we try to compute some h_n(X) in the steenrod algebra by applying
brutally the corresponding relations in the symmetric functions
\begin{Mexin}
hnrelation := 
proc(X,Y,degree)
    local e;
begin
    e:=OREtx(X.Y);
//    print(pretty(toe(S::h([degree-i]),e,X.Y)),pretty(toe(S::e([i]),e,Y))) $i=0..degree;
    _plus((-1)^i*toe(S::h([degree-i]),e,X.Y)*toe(S::e([i]),e,Y) $i=0..degree);
end_proc:
\end{Mexin}
It does not work: hnrelation(X,Y,d) depends on Y for d=3, Y=[y1,y2,y3]

\subsection{Trying to reverse engineer the  Gröbner basis computed by maple}

An element in the right ideal generated by p1 and p2 whose leading
term in K[x1,...,xn] is x_{n-1}^2;

Ouch, only $c$ below works!

\begin{Mexin}
n:=4:

e:=OREtx([x.i $i=1..n]):
(hold(p).k:=qPPk(e,k)) $k=1..n:

a:=mapcoeffs(-p1^2 + p2*(1+q+q*tx1) + 2*p1*(_plus((1+q*tx.i)*x.i$i=2..n)), factor):
b:=mapcoeffs(p2*(1+q+q*tx1) + p1*(_plus((1+q*tx.i)*x.i$i=2..n)-(1+q*tx1)*x1), factor):
c:=mapcoeffs(p2*(1+q+q*tx1) + p1*p1 - 2*p1*(1+q*tx1)*x1, factor):
pretty(c);

prog::test(iszero(a-b)):
prog::test(iszero(a-c)):
prog::test(lterm(expr(a), [x.i $i=1..n], LexOrder), 
	   x2^2):
pretty(a);
\end{Mexin} %$
This is precisely the element of the gröbner basis computed with maple.
So, one should be able to obtain as well f(gbu[3]) and so on! Why doesn't
the following work?

Post comment: be careful! multiplying on the right by symmetric
functions is not enough!

\begin{Mexin}
//testargs(FALSE):Pref::typeCheck(None):
S:=examples::SymmetricFunctions():
export(experimental::Steenrod):
bla :=
proc(X,Y,d)
begin
    A:=OREtx(X.Y):
    [g,ind]:=genericHomogeneousMultiple([op(map(X.Y,A)),qPk(A,k) $ k=1..d
					 ],
					d,
					hold(a),
					vec->(_plus(op(vec,1..d))>0)):
    elements := [hold(p).k $ k=1..d, op(X),op(Y)];
    sys := [coeff(g-qPk(A,d,Y))];
    [g,
     _plus(de*_mult(op(zip(elements, [op(de,2..nops(de))], _power)))
	   $ de in ind),
     solve(sys,
	   ind)]
end_proc:
bla([x1,x2],[y1],3);
\end{Mexin} %$

A function that computes the element of degree d in the right ideal
generated by the pk(X+Y) which annihilates all the variables in X:
\begin{itemize}
\item The solution seems to be fairly unique.
\item Trying to eliminate any of the used generator fails in degree 3
\item Similar computations in OREtx fails or are not well defined
\end{itemize}
%Don't work with MuPAD 2.0.0 ?
\begin{Mexin}
testargs(FALSE):Pref::typeCheck(None):
S:=examples::SymmetricFunctions():
export(experimental::Steenrod):
computeAnnihilator :=
proc(X,Y,d)
begin
    A:=ORExs(X.Y):
    [g,ind]:=genericHomogeneousMultiple([qPPk(A,k)    $ k=1..d,
					 if X<>X.Y then
					     qPPk(A,k,X)  $ k=1..d-1;
					 else
					     null()
					 end_if,
					 S::e([k])(X) $ k=1..min(d-1,nops(X))],
					d,
					hold(a),
					vec->(_plus(op(vec,1..d))>0)):
    elements := [hold(p).k $ k=1..d,
		 if X<>X.Y then
		     hold(px).k  $ k=1..d-1;
		 else
		     null()
		 end_if,
		 hold(e).k $ k=1..min(d-1,nops(X))];
    [g,
     _plus(de*_mult(op(zip(elements, [op(de,2..nops(de))], _power)))
	   $ de in ind),
     solve([coeff(select(expand(expr(g)),t->_or(has(t,v) $v in X)),X.Y)],
	   ind)]
end_proc:
\end{Mexin}
\begin{Mexin}
[g, formalg, res] := computeAnnihilator([x1], [y1,y2,y3], 2);
mapcoeffs(poly(normal(subs(formalg, op(res,1), a[2,0,0,0]=1,a[1,0,1,0]=0)*(2*q+sx1)), [p1,p2,px1,e1]), factor);
subs(expr(g), op(res,1));
Otx:=OREtx([x1,y1,y2,y3]):Oxs:=ORExs([x1,y1,y2,y3]):
pretty(expr(Otx(Oxs(normal(g2)*(2*q+sx1)))))

[g, formalg, res] := computeAnnihilator([x1,x2], [y1,y2,y3], 3);
mapcoeffs(poly(normal(subs(formalg, op(res,1), a[2,0,0,0]=1,a[1,0,1,0]=0)*(2*q+sx1)), [p1,p2,px1,e1]), factor):
Otx:=OREtx([x1,x2,y1,y2,y3]):Oxs:=ORExs([x1,x2,y1,y2,y3]):
g:=normal(subs(expr(g), op(res,1),a[3,0$6]=1)):
g:=g*denom(g):
pretty(expr(Otx(Oxs(normal(g)))))

(res := computeAnnihilator([x1,x2], [y1,y2,y3,y4], 3):1;)

(res := computeAnnihilator([x1,x2,x3], [y1], 4); nops(res));
write(Text, "Annihilator4.mu", res):

\end{Mexin} %$


\subsection{blabla}

\begin{Mexin}
reset():
testargs(FALSE):Pref::typeCheck(None):
export(experimental::Steenrod):
generic:=genericWeylElementWithSymmetries:
computeAnnihilator :=
proc(X,Y,d,order)
    option remember;
begin
    // Build a generic element of degree d in the right ideal
    // generated by the conjugated steenrod operators qPPk
    //                  ********** FALSE!!!!
    // in ORExt(X.Y);
    // The order of the right multipliers does not exceed order-1.
    // The element is symmetric in X and Y
    userinfo(0, "Building equations");
    W := ORExt(X.Y):
    makeUnknowns := proc(t,f) option remember; begin f[op(t)] end;
    Fs := [ _plus(generic(W, d-i,o,[nops(X), nops(Y)], makeUnknowns, hold(f).i)
                  $ o=0..order-1)
            $ i=1..d ];
    F := _plus(qPPk(W,i) * Fs[i] $ i=1..d);
    unknowns := sort(valuesOfRememberTable(makeUnknowns),
                     (x,y)->if op(x,1)<>op(y,1)
			    then bool(expr2text(op(x,1))>expr2text(op(y,1)))
                            else combinat::words::degLexLess(op(y,3),op(x,3));
                            end_if);
    // Get the coefficients of all terms containing a variable in X
    // If Y is non empty, ask for the coefficient of y1^n to be 1
    equations := {coeff(t) $ t in select(F::dom::monomials(F),
                                         t -> has(expr(t),X)),
                  if nops(Y)>0 then
                      coeff(F, polylib::lsupport(W(Y[1]^d))) - 1
                  else null()
		  end_if
                 };

    userinfo(0, "Solving system with ". nops(equations)." equations and "
		 . nops(unknowns) . " unknowns"):
    linsolve(equations, unknowns);
end_proc:
\end{Mexin} %$

First experiment with Y the empty set (that is we are looking for a
non-trivial relation in the symmetric ideal generated by the qPPk):
\begin{Mexin}
n:=2: m:=0: computeAnnihilator([x.i $i=1..n], [y.i $i=1..m], n+1, 1) ;
\end{Mexin}

Oh yes, there is one!

\begin{remark}
  The qPPk satisfy the newton identities. That is, for any k, and any
  alphabet of $n$ variables:
  \begin{displaymath}
    \sum_{i=0}^n (-1)^i qPP_{k+n-i} e_i = 0 .
  \end{displaymath}
\end{remark}
\begin{proof}
  That's exactly the same trick as for the classical newton identities:
  \begin{displaymath}
    \sum_{i=0}^n (-1)^i qPP_{n+k-i} e_i =
    \sum_{i=0}^n (-1)^i \sum_x (1+q x\p_x) x^{k+n-i} e_i =
    \sum_x (1+q x\p_x) x^k (\sum_{i=0}^n (-1)^i x^{n-i} e_i) =
    \sum_x (1+q x\p_x) x^k 0 = 0
  \end{displaymath}
\end{proof}

Checking this remark by computations in ORExt and in the Weyl algebra:
\begin{Mexin}
export(experimental::Steenrod):
vars := [x,y,z,t]: n:=nops(vars): k:=3:
S:=examples::SymmetricFunctions():
W := ORExt(vars):
normal(_plus((-1)^i*qPPk(W,n+k-i)*W(S::e([i])(vars)) $i=0..n));
\end{Mexin} %$
\begin{Mexin}
export(experimental::Steenrod):
vars := [x,y,z,t]: n:=nops(vars):
S:=examples::SymmetricFunctions():
W := Dom::WeylAlgebra_Table(vars, Dom::ExpressionField()):
normal(_plus((-1)^i*qPPk(W,n+k-i)*W(S::e([i])(vars)) $i=0..n)) $ k=0..5;
\end{Mexin}

So what, what's the point ... The qPPk are *NOT* the conjugates of the
qPk by the euler operator!!! And there does not seem to be any
reasonable newton identities between the qPk's themselves.

variables:=[x1,x2]:
W:=Dom::WeylAlgebra_Table(variables, Dom::ExpressionField()):
for k from 0 to 5 do 
    hold(e).k:=W(S::e::genericElement(p->hold(a)[op(p)],k)(variables)) 
    hold(P).k :=_plus(W(v)^k * W(v*D(v)) $ v in variables);
end_for:
solve([coeff(P3*e1+P2*e2+P1*e3)])


There seems to be a non trivial
\begin{Mexin}
n:=2: m:=2: computeAnnihilator([x.i $i=1..n], [y.i $i=1..m], n+1, 2) ;
\end{Mexin}


\section{Finding a quasi elementary function}

We look for a non trivial zero-element of degree n+1 in the ideal
generated by <steen> in k[x1..xn].

\begin{Mexin}
e:=OREtx([x]):
(hold(p).k:=qPPk(e,k)) $k=1..n+2:
solve([coeff(expr(a*p1*p1+b*p2),[x1])], [b,a]) ;
\end{Mexin} % $
\begin{Mexout}
piecewise({[b = - a + a q - a q tx1]} if q tx1 + 1 <> 0,

   {[]} if q tx1 + 1 = 0)
\end{Mexout}
Variation using ORExs 
\begin{Mexin}
A:=ORExs([x]):
(hold(p).k:=qPPk(A,k)) $k=1..3:
[g,ind]:=genericHomogeneousMultiple([p1,p2],2):
map(solve([coeff(expr(g),[x])], ind), map, map, factor)
\end{Mexin} % $

In other words, in one variable, we have the following relation:
\begin{Mexin}
p1^2 - p2 * (q+sx) = 0
\end{Mexin}

Some variations in degree 3 in two variables:
\begin{Mexin}
A:=ORExs([x,y]):
(hold(p).k:=qPPk(A,k)) $k=1..3:
map(solve([coeff(expr(normal(p3*a+p2*p1*b+p1^3*c+p2*(x+y)*d+p1^2*(x+y)*e)),[x])], [a,b,c,d,e,f,g]),map, map, factor)
\end{Mexin}
\begin{Mexin}
A:=ORExs([x,y]):
(hold(p).k:=qPk(A,k)) $k=1..3:
map(solve([coeff(expr(normal(p3*a+p2*p1*b+p1^3*c+p2*(x+y)*d+p1^2*(x+y)*e+p1*(x^2+y^2)*f+p1*x*y*g)),[x])], [a,b,c,d,e,f,g]),map, map, factor)
\end{Mexin}

\subsection{Studying the usual identity}

If b=e or h, and B is the antipode basis, then they satisfy an identity
such as:

Bk(Y) = \sum_{i=0..k} -1^(i) Bi(X+Y) * bi(X)

In the classical case, this yields a relation modulo the ideal <Sym+>
generated by the symmetric functions between Bk(Y) and
bk(X). Furthermore, taking b=e, and k>|X| yields that bk(X)=0, so
Bk(Y) is in the ideal.

Can we find a proper antipode in Steenrod ?
Does this identity still hold ?
Then, Bk(Y)=bk(X) mod <Steen+>
The last step as is won't work, since there is no elementary element.

Problem: is Bk(Y)+bk(Y) zero ? is it nonzero as soon as k>|X| ?

base :=i->if i<>0 then S::p([i]) else nops; end:
abase:=i->if i<>0 then S::antipode(base(i)) else nops; end;
X:=[x1,x2,x3]: Y:=[y1,y2]:
print(expand(_plus((-1)^i*base(k-i)(X.Y)*abase(i)(Y) $i=0..k))) $k=3..3;

base:=i->S::h([i]):X:=[x1,x2,x3]: Y:=[y1,y2]: print(expand(_plus(base(k-i)(X.Y)*S::antipode(base(i))(Y) $i=0..k))) $k=3..3;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Trying to use indeterminates for P1(Y), P2(Y), ...


//testargs(FALSE):Pref::typeCheck(None):
S:=examples::SymmetricFunctions():
export(experimental::Steenrod):
bla :=
proc(X,d)
begin
    A:=ORExs(X):
    (degree(A(hold(PY).k)) := k) $k=1..d;
    [g,ind]:=genericHomogeneousMultiple
    (revert([qPPk(A,k)+A(hold(PY).k) $ k=1..d]).
     [op(map(X,A)), A(hold(PY).k) $k=1..1],
     d,
     hold(a),
     vec->(_plus(op(vec,1..d))>0)):
    elements := revert([hold(p).k $ k=1..d]).
		 [op(X), hold(PY).k $k=1..1];
    sys := [coeff(select(expand(expr(g)),t->_or(has(t,v) $v in X)),
		  X.[hold(PY).k $k=1..d])];
    [g,
     _plus(de*_mult(op(zip(elements, [op(de,2..nops(de))], _power)))
	   $ de in ind),
     solve(sys,
	   ind)]
end_proc:


[g,gformal,res] := bla([x],2);
collect(normal(subs(expr(g), op(res,1), a[2,0,0]=2*q+sx)), [PY1, PY2], factor)

q*qPPk(3,Y) is in the right ideal generated by the qPPk(d, X+Y) ? It
seems to be the only one ?
\begin{Mexin}
[g,gformal,res] := bla([x1,x2],3);
res2:=op(solve([op(op(res,1)), a[1,0,0,2,0,0,0]=0, a[1,0,0,0,2,0,0]=0,a[2,0,0,1,0,0,0]=0,a[2,0,0,0,1,0,0]=0, a[1,1,0,0,0,0,0]=-1,a[2,0,0,0,0,1,0]=0], ind),1);
factor(subs(expr(g), op(res,1)));
factor(subs(expr(g), res2));
collect(expand(subs(gformal, res2)), [p1,p2,p3]);

res2:=op(solve([op(op(res,1)), a[0,0,1,2,0,0,0]=0, a[0,0,1,0,2,0,0]=0,a[0,0,2,1,0,0,0]=0,a[0,0,2,0,1,0,0]=0, a[0,0,2,0,0,1,0]=0,a[0,1,1,0,0,0,0]=1], ind),1)
\end{Mexin}

[g,gformal,res] := bla([x],2);
collect(normal(subs(expr(g), op(res,1), a[2,0,0]=2*q+sx)), [PY1, PY2], factor)


X:=[x1,x2]:
d:=3:


There seems to be no X-anihilator in the right ideal generated by qSteen,
when one multiplies only by elements of X ?
\begin{Mexin}
bla :=
proc(X,Y,d)
begin
    A:=ORExs(X.Y):
    [g,ind]:=genericHomogeneousMultiple
    (revert([qPPk(A,k) $ k=1..d]).
     [op(map(X,A)), qPPk(A,k,X) $k=1..d-1],
     d,
     hold(a),
     vec->(_plus(op(vec,1..d))>0)):
    elements := revert([hold(p).k $ k=1..d]).
		 [op(X), hold(pX).k $k=1..d-1];
    sys := [coeff(select(expand(expr(g)),t->_or(has(t,v) $v in X)),
		  X.Y)];
    [g,
     _plus(de*_mult(op(zip(elements, [op(de,2..nops(de))], _power)))
	   $ de in ind),
     solve(sys,
	   ind)]
end_proc:
[g,gformal,res]:=bla([x1,x2],[y1,y2,y3,y4],3);
factor(subs(expr(g), op(res,1)));
\end{Mexin}

\begin{Mexin}
bla :=
proc(X,Y,d)
begin
    A:=ORExs(X.Y):
    [g,ind]:=genericHomogeneousMultiple
    (revert([qPPk(A,k) $ k=1..d]).
     [op(map(X.Y,A))],
     d,
     hold(a),
     vec->(_plus(op(vec,1..d))>0)):
    elements := revert([hold(p).k $ k=1..d]).
		 [op(X.Y)];
    sys := {coeff(select(expand(expr(g)),t->_or(has(t,v) $v in X)),
		  X.Y)};
    [g,
     _plus(de*_mult(op(zip(elements, [op(de,2..nops(de))], _power)))
	   $ de in ind),
     solve(sys,
	   ind)]
end_proc:
[g,gformal,res]:=bla([x1,x2],[y1,y2,y3,y4],3);
factor(subs(expr(g), op(res,1)));
\end{Mexin}


\section{Looking for higher derivations and elementary/homogeneous functions}

Generic higher derivation of degree $d$ in the Stenrod algebra:
  \begin{Mexin}
genericHigherDerivation(S::p, 1)
  \end{Mexin}
  \begin{Mexout}
                           a[1] p[1]
  \end{Mexout}
  \begin{Mexin}
genericHigherDerivation(S::p, 2)
  \end{Mexin}
  \begin{Mexout}
                 /     2 \
                 | a[1]  |
                 | ----- | p[1, 1] + a[2] p[2]
                 \   2   /
  \end{Mexout}
  \begin{Mexin}
genericHigherDerivation(S::p, 3)
  \end{Mexin}
  \begin{Mexout}
    /     3 \
    | a[1]  |
    | ----- | p[1, 1, 1] + a[3] p[3] + (a[1] a[2]) p[2, 1]
    \   6   /
  \end{Mexout}
  \begin{Mexin}
genericHigherDerivation(S::p, 4)
  \end{Mexin}
  \begin{Mexout}
/     4 \                 /     2      \
| a[1]  |                 | a[1]  a[2] |
| ----- | p[1, 1, 1, 1] + | ---------- | p[2, 1, 1] +
\   24  /                 \      2     /

   /     2 \
   | a[2]  |
   | ----- | p[2, 2] + (a[1] a[3]) p[3, 1] + a[4] p[4]
   \   2   /
  \end{Mexout}

That's easy to reverse engineer: the result can be expressed as
$$\sum_\lambda
  \prod_i a_i^{m_i(\lambda)} \frac{p_\lambda}{\prod_i m_i(\lambda)!},$$
where $m_i(\lambda)$ is the number of parts of size $i$ in $\lambda$.

In other words, the generating series of the higher derivation is
obtained from the one of the primitive elements by exponentiation:
$$\sum SQ_d t^d = \dots \exp(a_d p_d t^d) \dots \exp(a_2 p_2 t^2) \exp(a_1 p_1t)$$

That's not quite a surprise: for a series of operators, being an
higher derivation is equivalent to having a group-like generating
series. Furthermore, the exponential of a lie like element is group
like; the product of two group-like elements is group-like.

Let us check it:
\begin{Mexin}
reset():S:=experimental::SteenrodAlgebra():
p:= S::p::domainWrapper:
(hold(p).i := S::p([i])) $ i=1..10:
\end{Mexin} %$

\begin{Mexin}
bool(poly(expr(truncate(_mult(ncexp(p[k]*a[k]*t^k)
                              $ k in revert([$1..ORDER-1])))), [t]) = 
     poly(hold(p)[] +
          _plus(expr(genericHigherDerivation(S::p, k)) * t^k
                $ k in revert([$1..ORDER-1])), [t]))
\end{Mexin}

Furthermore, if we want the higher derivations to be pure wedges, or
in other words of homogeneous order (as the real steenrod squares),
one can just take $a_1=1, a_2=-1, a_3=1, \dots$. Then, the generating
series gives the Steenrod squares $SQ_d$ of total degree $d$, degree
$2d$ and order $d$.

\begin{Mexin}
W := Dom::WeylAlgebra_Table([x.i $i=1..5], Dom::ExpressionField());
e4 := mapcoeffs(S::p::toWeylOriginal(W, genericHigherDerivation(S::p,4)), subs, a[1]=1, a[2]=-1, a[3]=1, a[4]=-1);
\end{Mexin} %$

The image on $\pi$ of this operator should be the elementary symmetric
function (check and debug the code!):

\begin{Mexin}
W::onPoly(e4, x1*x2*x3*x4*x5) ;
\end{Mexin}

What does its antipode look like ?
\begin{Mexin}
h4:=mapcoeffs(S::p::toWeylOriginal(W, antipode(genericHigherDerivation(S::p,4))), subs, a[1]=1, a[2]=-1, a[3]=1, a[4]=-1) ;
W::onPoly(h4, x1*x2*x3*x4*x5) ;
\end{Mexin}

Conversion original Steenrod -> ORExt:
\begin{Mexin}
f:=operators::makeLinear(p->_mult( _plus(W(x^k)*W(t.x) $ x in [x1,x2,x3,x4]) $ k in p));
\end{Mexin}

Let us see how this works in the case of symmetric functions. Here we
look for the proper coefficients so that we get the generating series
of the homogeneous functions:
\begin{Mexin}
s:=truncate(_mult(ncexp(Sym::p([k])*a[k]*t^k) $ k in revert([$1..ORDER-1]))):
mapcoeffs(s, subs, a[k]=-(-1)^k $k=1..ORDER -1)
solve([coeff(Sym::h(s)-Sym::h([])-_plus(Sym::h([k]) $k=1..ORDER-1))], [a[k] $k=1..ORDER-1])
\end{Mexin} %$
This is: $a[k] := 1/(k*t^k)$.

Same thing for the elementary symmetric functions:
\begin{Mexin}
solve([coeff(Sym::e(s)-Sym::e([])-_plus(Sym::e([k]) $k=1..ORDER-1))], [a[k] $k=1..ORDER-1])
\end{Mexin}
This is: $a[k] := (-1)^(k+1)/(k*t^k)$.

Let us turn back to the $q$-Steenrod algebra. By a proper choice of the
a[i]'s, one might be able to ensure that the pure polynomial part
consists of the homogeneous/elementary functions:

normal(mapcoeffs(S::p::toWeyl(W, genericHigherDerivation(S::p,2)), subs, a[1]=1, a[2]=(1-q)/2, a[3]=1, a[4]=-1))


The expansion of exp( ... ):
\begin{Mexin}
poly(expr(ncexp(_plus(p[k]*a[k]*t^k $ k=1..ORDER-1))), [t]):
\end{Mexin} %$


Alternatively, one can use $\exp(\sum_d a_d p_d t^d)$ to see what
happens. The resulting generating series will be group like. Will the
coefficients still have a nice look?

Conversion _mult(exp) -> exp(_plus)

\begin{Mexin}
poly(expr(nclog(truncate(_mult(ncexp(p[k]*a[k]*t^k) $ k in revert([$1..ORDER-1]))))), [t]):
\end{Mexin}

The reverse conversion _mult(exp) -> exp(_plus) seems easier to
reverse engineer (all coefficients are positive, ...)
\begin{Mexin}
eqns := map({coeff(ncexp(_plus(p[k]*a[k]*t^k $ k=1..ORDER-1)) - 
                   truncate(_mult(ncexp(p[k]*b[k]*t^k)
                                  $ k in revert([$1..ORDER-1]))))}, 
	    coeff, [t]):
vars := [b[k] $ k=1..ORDER-1]:
res := op(solve(eqns, vars), 1):
print(eq) $ eq in res;
\end{Mexin} %$

The power of $q$ in a term of degree d in the a[i]'s is d-1.  For
terms of the form a[i]*a[j], the coefficient is $(j-i)/2$. Actually
one probably just need to look at the commutations when renormalizing
the expansion of ncexp(plus(...)).

The antipode equation works as expected:
\begin{Mexin}
d:=5:
normal(_plus(genericHigherDerivation(S::p, k) *
             operators::antipode(genericHigherDerivation(S::p, d-k)) $k=0..d))
\end{Mexin} %$

\section{Questions for Florent}

How much freedom do we have for the coproduct ?

What are the homogeneous functions characterized aside from the shape
of the coproduct ?

Any formula that would take into account the dual elements?

Can we find an operator dx1 such that D_k dx1 belongs to the left
ideal generated in the Weyl algebra by the Steenrod algebra?

What's the exponential relation between the p's and the h's (or the
e's) in NCSF?

Existe-t'il un coproduit dans Weyl (ou ORExs, ...) tel que les
polynômes soient un module de Hopf ? Peut-être $\Delta(x):=x\tensor 1$
and $\Delta \p:= \p \tensor 1 + 1 \tensor \p$ ?

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Looking for Steenrod-conjugated operator}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Let $t_i:=x_i\p_i$ be the Euler operator in the variable $x_i$, and
consider the "Ore-algebra" $\oxt:=\K[x_i][[t_i]]$ generated by the
$x_i's$ and $t_i$'s (we include infinite series in the $t_i$'s; maybe
we should only include those that converge for all non-negative
integers).

\begin{remark}[Computation rules $\K[x_i][[t_i]]$]
  Let $f$ be a univariate series in $\K[[t_i]]$, $g$ be a multivariate
  series in $\K[[t_i]]$ and $x^d$ a monomial in $\K[x_i]$. Then, we
  have the following derivation and commutation rules:
  \begin{align*}
    t_i    . x_i         &= x_i\\
    t_i      x_i         &= x_i (t_i+1)\\
    t_i    . x_i^k       &= k x_i^k\\
    t_i      x_i^k       &= x_i^k (t_i+k)\\
    f(t_i) . x_i^k       &= f(k) x_i^k\\
    f(t_i)   x_i^k	 &= x_i^k f(t_i+k)\\
    g(t_1,\dots,t_n).x^d &= g(d_1,\dots, d_n) x^d \\
    g(t_1,\dots,t_n) x^d &= x^d g(t_1+d_1,\dots, t_n+d_n)\\
  \end{align*}
  Some examples:
  \begin{align*}
    \exp(t_i) x_i^k      &= \exp(k) x_i^k \exp(t_i)\\
    (t_1+\dots_t_n) x^d  &= (d_1+\dots+d_n) x^d\\
    (1+qt_i) x^d         &= x^d (1+q(t_i+d_i))
  \end{align*}
\end{remark}

\subsection{Steenrod-conjugated operators}

Let $f$ and $g$ be two operators in $\oxt$. We say that $f$ \emph{has a
  right conjugate w.r.t. $g$} if $fg$ belongs to the right ideal $g
\oxt$. Furthermore, $f$ is Steenrod-conjugated if $f\qsteen\subset
\qsteen \oxt$; an equivalent condition is that, $f\P_{1,q}$ and
$f\P_{2,q}$ both belong to $\qsteen \oxt$.

\begin{remark}
  The Steenrod-conjugated operators form an algebra which contains
  $\qsteen\oxt$. Furthermore, they stabilizes the hits.
\end{remark}

\begin{problem}
  Is the Steenrod-conjugated algebra generated by $n$ homogeneous
  elements of degree $1$? Does this algebra have the same dimension as
  a polynomial algebra?
\end{problem}

An operator which commutes with any Steenrod operator (or equivalently
with both $\P_{1,q}$ and $\P_{2,q}$) is Steenrod-conjugated. But there
are no such operators.

\subsection{Finding homogeneous elements a and b of degree 1 in the Weyl
  algebra on 1 variable such that a*p1 = p1*b and a*p2 = p2*b}

\begin{Mexin}
n := 1:
S:=experimental::SteenrodAlgebra():
W:=Dom::WeylAlgebra_Table([x.i $i=1..n], Dom::ExpressionField()):
pW := k->S::p::toWeyl(W,S::p([k])):
order := ORDER:
generic := (a, order)-> W(x1)* (_plus(a.d * W(x1)^d*W(D(x1))^d $d=0..order+2)):
a:=generic(hold(a),order): 
b:=generic(hold(b),order):
truncateOrder :=
proc(p)
begin
    _plus(t $ t in select(p::dom::monomials(p), t->p::dom::order(t) < ORDER));
end:
eqns := [coeff(truncateOrder(a*p-p*b,order)) $ p in [pW(1), pW(2)]]:
vars := [(hold(a).i, hold(b).i) $ i in revert([$0..order+2])]:
solve(eqns, vars)
\end{Mexin}

Result: it's unique, up to the choice of $a_0$ and $b_0$. It's
completely unique if one further wants the result to specialize at q=0
to $a=b=x_1$.

\begin{displaymath}
a := 1 - q \sum_{k=1}^\infty \frac{(-1)^k}{k!}                            x^k\p^k\\
b := 1 +   \sum_{k=1}^\infty \frac{(-q)^{k+1}}{(2q+1)(3q+1) \dots (kq+1)} x^k\p^k
\end{displaymath}

\subsection{}

Let $A$ be the set of operators in $\otx$ that have conjugates w.r.t.
$\P_{1,q}$ and $\P_{2,q}$. It's an algebra. Key point: does there
exists operators in $A$ of degree $1$ which are not symmetric ?

\subsubsection{Degree $0$}

\begin{proposition}
  Take $f\in \K[[t_i]]$. Then, $f$ has a conjugate w.r.t. $\P_{l,q}$
  iff $f(t_1,\dots,t_i+l,\dots,t_n)=f(t_1,\dots,t_j+l,\dots,t_n)$ for
  all $i\ne j,t_1,\dots,t_i$. In other words, $f$ is periodic w.r.t.
  all vectors $(-k,0\dots 0,k,0,\dots,0)$.
  
  As a consequence, tfae:
  \begin{enumerate}
  \item $f$ is in $A$;
  \item $f$ has conjugates w.r.t. all $\P_{l,q}$;
  \item $f$ has a conjugate w.r.t. $\P_{1,k}$;
  \item $f$ is periodic w.r.t. all vectors $(-1,0\dots
    0,1,0,\dots,0)$.
  \end{enumerate}
  
  In particular, $A\cap\K[t_i]=\K[t_1+\dots+t_n]$: the polynomial
  operators in $A$ are polynomials in the symmetric Euler operator.
\end{proposition}

\begin{proof}
  Take a polynomial $p$ in $A$ and expand it on the variables
  $t_1+\dots+t_n, t_2-t_1,\dots,t_n-t_1$. It's periodic and thus
  constant w.r.t. all variables $t_i-t_1$. Hence, $p$ is a polynomial
  in $t_1+\dots+t_n$.
\end{proof}

This yields plenty of operators of degree $0$ which stabilize the
hits (and whose adjoints stabilize the harmonics).

\begin{proposition}
  Take $f\in\K[[t_i]]$, and assume $n>1$. Then, there exists $g$ in
  $\K[[x_i,t_i]]$ such that $x_1 f P_{k,q}=P_{k,q} g$ if, and only if,
  \begin{displaymath}
    (1+qx_1) f(x_1+l,\dots,x_n) = (1+q(x_1+1)) f(x_1,\dots,x_i+l,\dots,x_n)
  \end{displaymath}
  
  Hence there is no operator in $A \cap x_1\K[[t_i]]$, at least among
  the series that converge (compute $f(-2,2)$ from $f(0,0)$ by using
  on one side the $(-1,1)$ "period" twice and on the other side the
  $(-2,2)$ "period").
\end{proposition}


\begin{proposition}
  With just one variable, any function of the form $x_1f(t_1)$ has a
  conjugate with respect to any the positive Steenrod powersum:
  \begin{displaymath}
    x_1 f(t_1) P_{k,q} = P_{k,q} x_1 (1- \frac{q}{1+q(t_1+1)})
    f(t_1+l)
  \end{displaymath}
\end{proposition}

\begin{remark}
  In $1$ variable and degree $\geq1$, $P_1$ generate all of
  $\k[x_i][[t_i]]$ as a right ideal.
  
  In $2$ variables and degree $\geq3$, $P_1,P_2,P_3$ generate all of
  $\k[x_i][[t_i]]$ as a right ideal.
\end{remark}


\subsubsection{Looking at the right ideal generated by the $P_{q,k}$}

reset():
export(experimental::Steenrod):
S:=experimental::SteenrodAlgebra():
W:=ORExt([x1,x2,x3]):
(P.i:=W(S::p([i]))) $i=1..3:
toP := f->mapcoeffs(poly(expr(f), [x1,x2,x3]), factor):
[rank,unrank] := combinat::rankers::fromList(combinat::integerVectors(3,3)):
M:=_concat(toVector(toP(x), rank) $ x in homogeneousRightIdealBasis([P1,P2,P3],[W(x1),W(x2),W(x3)],3)):
subs(M, tx.i=1/q*(sx.i-1) $i=1..3) // For a more compact output


Looking for an element of the form x1 f() such that x1 f() P_2 belongs
to the right ideal generated by P_1 and P_2.

\begin{Mexin}
reset():
n:=3:
variables := [x.i $i=1..n]:
export(experimental::Steenrod):
S:=experimental::SteenrodAlgebra():
W:=ORExt(variables):
(P.i:=W(S::p([i]))) $i=1..3:
toP := f->poly(expr(f), variables):
makeUnknowns := proc(x,y) option remember; begin hold(a)[op(x),op(y)] end:
f := toP(homogeneousRightIdealGeneric([P1,P2],[W(x) $ x in variables],3,makeUnknowns)):
g := poly(_plus(x^2 * variables[1]*hold(g).x $ x in variables), variables):
unknowns := [hold(g).x $x in variables].map([op(extop(makeUnknowns,5))], op, 2):
equations := map([coeff(f-g)], expand):
res:=table(op(linsolve(equations, unknowns)));
// For n=3, this choice of a[1,0,...,2] gives a fraction-free solution
// that is symmetric in tx2, tx3, ...
res2:=map(subs(res,
               a[2, 1, 0$n-1] = 0, a[1,1,1]=0, // for the degenerate case n=2
               a[1, 0$n-1, 2] = 2*(1+q*tx.n)/(1+q*tx.1)*_mult((1+q*(t.x+1)) $x in variables)),
          factor):
print(res2[hold(g).x]) $ x in variables;
\end{Mexin}

It seems that we do not have much choice. Such a f should satisfy:
\begin{displaymath}
  f(t1+2) = a g1 + b (1 + q(tx1+1))\\
  f(ti+2) = a gi + b (1 + qtxi),
  f(t1+1) = c (1 + q(tx1+1))\\
  f(ti+1) = c (1 + qtxi),
\end{displaymath}
for some choice of $a$, $b$, and $c$ in $\K[[t_i]]$.

Ouch, for $n=4$, there is just the periodic solution!

Looking for an element of degree 1 g=\sum x g_x such that g P_2
belongs to the right ideal generated by P_1 and P_2. We write
g(x,y)=g_x(y+2)


First try: we only multiply $P_1$ and $P_2$ by polynomials in the
$x_i$'s that are symmetric w.r.t. $2,\dots,n$.  Note: an element in
this right ideal which is symmetric (say w.r.t.  $x_2,\dots,x_n$) is
not necessarily in the linear span of the $P_i*e_{3-i}(x_1,\dots,x_n)$
w.r.t. coefficients in $\K[[t_i]]$.  Indeed the action on the $x_i$'s
also change the coefficients! Note for example that we won't get
$g=P_1$ this way!

\begin{Mexin}
reset():
n:=4:
var1 := x1;
others := [ x.i $ i=2..n]:
variables := [var1].others:
export(experimental::Steenrod):
S:=experimental::SteenrodAlgebra():
elem := (examples::SymmetricFunctions())::e::domainWrapper:
W:=ORExt(variables):
(P.i:=W(S::p([i]))) $i=1..3:
toP := f->poly(expr(f), variables):
makeUnknowns := proc(x,y) option remember; begin hold(a)[op(x),op(y)] end:
fgen := toP(homogeneousRightIdealGeneric([P1,P2],[W(var1),W(elem[d](others)) $d=1..min(n-1,3)],3,makeUnknowns)):
ggen := poly(_plus(x*y^2*(1+q*t.y)*g(x,y) $ x in variables $ y in variables ), variables):
unknowns := [g(x,y) $x in variables $y in variables].map([op(extop(makeUnknowns,5))], op, 2):
equations := map([coeff(fgen-ggen)], expand):
res:=table(op(linsolve(equations, unknowns))):
res2:=map(subs(res,
               a[2, 1, 0$n-1] = 0, a[1,1,1]=0, // for the degenerate case n=2
               a[1, 0$n-1, 2] = 2*(1+q*tx.n)/(1+q*tx.1)*_mult((1+q*(t.x+1)) $x in variables)),
          factor):
print(res2[hold(g).x]) $ x in variables;
\end{Mexin} %$
Hmm, what to do with this ? 



\begin{Mexin}
reset():
n:=4:
var1 := x1;
others := [ x.i $ i=2..n]:
variables := [var1].others:
export(experimental::Steenrod):
S:=experimental::SteenrodAlgebra():
elem := (examples::SymmetricFunctions())::e::domainWrapper:
W:=ORExt(variables):
(P.i:=W(S::p([i]))) $i=1..3:
toP := f->poly(expr(f), variables):
makeUnknowns := proc(x,y) option remember; begin hold(a)[op(x),op(y)] end:
fgen := toP(homogeneousRightIdealGeneric([P1,P2],[W(v) $ v in variables],3,makeUnknowns)):
ggen := poly(_plus(x*y^2*(1+q*t.y)*g(x,y) $ x in variables $ y in variables ), variables):
unknowns := [g(x,y) $x in variables $y in variables].map([op(extop(makeUnknowns,5))], op, 2):
equations := map([coeff(fgen-ggen)], expand):
res:=table(op(linsolve(equations, unknowns))):
\end{Mexin} %$

// get rid of many denominators
res2 := map(res, factor@_mult, (1+q*tx1)*(1+q*tx2)*(1+q*tx3)*(1+q*tx4)):
//Conditions imposed by symmetries + b.i = coeff of P1.xi^2;
res3:=map(res2, factor@subs,
          a[1,2,0,0,0]=b1,a[1,0,2,0,0]=b2,a[1,0,0,2,0]=b3,a[1,0,0,0,2]=b4,
          a[2,1,0,0,0]=c1,a[2,0,1,0,0]=c2,a[2,0,0,1,0]=c3,a[2,0,0,0,1]=c4,
          a[1,0,0,1,1]=(2*f2-f3-f4)/(1+q*tx1)/(1+q*tx2),a[1,0,1,0,1]=(2*f3-f2-f4)/(1+q*tx1)/(1+q*tx3))
\end{Mexin}

\subsubsection{Families of polynomials that behave nicely w.r.t. shifting}

Let $Q_k(t):=(1+qt)(1+q(t+1))\dots(1+q(t+k-1))$.
\begin{Mexin}
Q:=
proc(k:Type::NonNegInt, x=hold(t))
begin
    _mult((1+q*(x+i)) $ i = 0..k-1)/k! 
end_proc:
Q(3);
\end{Mexin} %$

The $Q_k$ satisfy the following recurrence relation (kind of divided
difference formula):
\begin{displaymath}
  Q_k(t+1)-Q_k(t) = q \frac{1+q(t+k-1)}{1+qt} Q_{k-1}(t)
\end{displaymath}
\begin{Mexin}
prog::test(iszero(normal((Q(k,t+1)-Q(k,t)) -q*(1+q*(t+k-1))/(1+q*t)*Q(k-1,t)))) $ k=1..5
\end{Mexin} %$

The product is given by:
\begin{displaymath}
Q_k(T)Q_l(t) = \sum_{j=\max(k,l)}{k+l} \binom{j}{j-k,j-l} Q_{k+l}(t) (-q)^{k+l-j}
\end{displaymath}
\begin{Mexin}
f:=(k,l)->expand(_plus(j!/(j-k)!/(j-l)!/(k+l-j)! *Q(j)*(-q)^(k+l-j) $j=max(k,l)..k+l)):
prog::test(iszero(simplify(f(k,l)-Q(k)*Q(l))))  $k=0..5 $l=0..5;
\end{Mexin} %$


\subsection{Finding a non-symmetric degree 1 Steenrod conjugated
  operator in the Weyl algebra}

We look for a Steenrod conjugated of the form $x_1$ times a linear
combination $f$ of $x_i^d t_i^d$.

Notes: the order of the terms in the Gi's do not decrease in the
commutations; the order of the terms in F can decrease by at most 2 in
the commutations.
\begin{Mexin}
reset():
vars := [x,y,z,t]:
W:=Dom::WeylAlgebra_Table(vars, Dom::ExpressionField()):
export(experimental::Steenrod):
S:=experimental::SteenrodAlgebra():
(P.i:=S::p::toWeyl(W, S::p([i]))) $i=1..3:
alias(generic=experimental::Steenrod::genericWeylElementWithSymmetries):
prettyres := truc -> subs(map(map(subs(truc, f=eval@(()->hold(f)[args()])), map, Dom::DistributedPolynomial(subs(unknownsF, f=eval@(()->hold(f)[args()])), Dom::ExpressionField())), map, expr@mapcoeffs, c->factor(c)*x), x=1):
\end{Mexin} %$

Attempt 1: F is of the form x*P, with P fully symmetric wrt t1,t2,t3
\begin{Mexin}
reset(): export(experimental::Steenrod):
vars := [x,y,z]:
W:=Dom::WeylAlgebra_Table(vars, Dom::ExpressionField()):
order := 4:
makeUnknownsF := proc(m) option remember; begin hold(f)(op(op(m,2))) end:
F :=W(x)*_plus(generic(W, 0,o,[nops(vars)],makeUnknownsF)
	       $ o=0..order+1):
unknownsF := sort(valuesOfRememberTable(makeUnknownsF),
                  (x,y)->combinat::words::degLexLess([op(y)],[op(x)])):
res := solveInSteenrodRightIdeal(W, F, unknownsF, order)
\end{Mexin} %$
Conclusion: yeah, there is one solution: it kills any Weyl operator
with all terms of positive degree in the x's:
\begin{displaymath}
  f=\prod_i \sum_d \frac{(-1)^d}{d!} x_i^d \p_i^d
\end{displaymath}

Attempt 2: F is of the form x*P, with P a series in the t^\mu which
is symmetric w.r.t. t2,t3,...
We use the exponential generating function for P:
P = \sum_\mu (-1)^\mu /\mu! f(\mu) t^\mu
\begin{Mexin}
reset(): export(experimental::Steenrod):
vars := [x,y,z]:
W:=Dom::WeylAlgebra_Table(vars, Dom::ExpressionField()):
alias(generic=experimental::Steenrod::genericWeylElementWithSymmetries):
setuserinfo(solveInSteenrodRightIdeal,1): setuserinfo(linalg,2):
order := 3:
makeUnknownsF := proc(l) option remember; begin hold(f)(op(l)) end:
makeCoeffsF   :=
proc(m)
    local k, n;
begin
    m:=op(m,2);
    n := _plus(op(m));
    (-1)^n/_mult(k! $ k in m)*makeUnknownsF(m);
end_proc:
F :=W(x)*_plus(generic(W, 0,o,[1,nops(vars)-1],makeCoeffsF)
	       $ o=0..order+1):
unknownsF := sort(valuesOfRememberTable(makeUnknownsF),
                  (x,y)->combinat::words::degLexLess([op(y)],[op(x)])):
res := solveInSteenrodRightIdeal(W, F, unknownsF, order, Divisible):
res2 := map(select(res, eq->bool(_plus(op(op(eq,1)))<order)),
            eq->lhs(eq)=factor(rhs(eq)));
print(x) $ x in res2;
\end{Mexin} %$

Attempt 2 bis on a specific specialization of f. This specialization
of f does not work for 4 variables or more.
\begin{Mexin}
reset(): export(experimental::Steenrod):
vars := [x,y,z]:
W:=Dom::WeylAlgebra_Table(vars, Dom::ExpressionField()):
alias(generic=experimental::Steenrod::genericWeylElementWithSymmetries):
setuserinfo(solveInSteenrodRightIdeal,1): setuserinfo(linalg,2):
order := 4:
makeUnknownsF := proc(l) option remember; begin hold(f)(op(l)) end:
makeCoeffsF   :=
proc(m)
    local k, n;
begin
    m:=op(m,2);
    n := _plus(op(m));
    (-1)^n*(2*op(m,1)*q+n-1)/_mult(k! $ k in m)*(n-2)/2*makeUnknownsF(m);
end_proc:
F :=W(x)*_plus(generic(W, 0,o,[1,nops(vars)-1],makeCoeffsF)
	       $ o=0..order+1):
unknownsF := sort(valuesOfRememberTable(makeUnknownsF),
                  (x,y)->combinat::words::degLexLess([op(y)],[op(x)])):
res := solveInSteenrodRightIdeal(W, F, unknownsF, order, Divisible):
res2 := map(res, map, factor):
print(x) $ x in res2;
\end{Mexin} %$

In 3 variables, that seemed to be a correct solution:
\begin{conjecture}
  The following operator is a degree 1 Steenrod conjugated:
  \begin{displaymath}
    f := x_1 \frac12
    \sum_{\mu} \frac{(-1)^\mu}{\mu!} (2q\mu_1 + |\mu|-1)(|\mu|-2)} t^\mu,
  \end{displaymath}
  where the sum runs through all degree vectors $\mu$.
\end{conjecture}
\begin{Mexin}
  SteenrodConjugated(W,3,1)
\end{Mexin} %$

Testing that this element indeed lies in the right ideal:
\begin{Mexin}
reset(): export(experimental::Steenrod):
vars := [x,y,z]:
W:=Dom::WeylAlgebra_Table(vars, Dom::ExpressionField()):
for order in [$1..10] do
    print(order, bool(solveInSteenrodRightIdeal(W, SteenrodConjugated(W,order+2), [], order)=[]))
end
\end{Mexin} %$
\begin{Mexout}
Info: Building equations
Info: Building unknowns
Info: Solving system with 8 equations and 7 unknowns
Info: assumptions: [] [4*(q + 1)^2, -3*q*(3*q + 2)]

                                  1, TRUE
Info: Building equations
Info: Building unknowns
Info: Solving system with 29 equations and 24 unknowns
Info: assumptions: [] [-q, -3*(q + 1)*(3*q + 2), 2*q*(4*q + 3), q^2*(4*q +\
 3), -(q + 2), -(q + 1)*(9*q + 8*q^2 + 2), -2*(q + 1)*(5*q + 5*q^2 + 1), (\
q + 1)*(2*q + 1), -(q + 1), -(2*q + 1)]

                                  2, TRUE
Info: Building equations
Info: Building unknowns
Info: Solving system with 71 equations and 58 unknowns
Info: assumptions: [] [(5*q + 4)*(2*q + 1)*(q + 1), 2*(2*q + 1)^2, 8*q + 7\
*q^2 + 2, -q*(5*q + 3)*(5*q + 2), -(7*q + 7*q^2 + 2), -(q + 1)*(5*q + 3), \
-(3*q + 2)*(2*q + 1), (q + 1)^2, q*(q + 1)*(5*q + 3), 2*(q + 1)*(2*q + 1),\
 -(21*q + 23*q^2 + 4), 3*q^2*(3*q + 2), -2*q*(4*q + 3), 4*(q + 1), -(5*q +\
 3), -(4*q + 3), 3*(q + 1), q^2*(4*q + 3), 2*(2*q + 1), -2*(2*q + 1), q*(5\
*q + 3), -2*q*(2*q + 1), -(2*q + 1), -(q + 1)]

                                  3, TRUE
Info: Building equations
Info: Building unknowns
Info: Solving system with 139 equations and 112 unknowns
Info: assumptions: [] [-3*(2*q + 1)^2, -q*(q + 1), -q*(5*q + 3)*(2*q + 1),\
 (4*q + 3)*(2*q + 1), (2*q + 1)*(5*q + 2), -q*(102*q + 337*q^2 + 497*q^3 +\
 268*q^4 + 12), 3*(q + 1)*(2*q + 1), -q*(5*q + 3)*(5*q + 2), -(7*q + 8*q^2\
 + 2), q*(q + 1), -(7*q + 7*q^2 + 2), -(3*q + 2)*(2*q + 1), 3*q^4, q*(q + \
1)*(5*q + 3), (q + 1)^2, 2*(q + 1)*(2*q + 1), -3*q^3*(4*q + 1)*(4*q + 3)*(\
q + 1), -3*q*(2*q + 1)*(3*q + 1)^2, 3*(q + 1)*(5*q + 2), -3*q, 2*q, -q, -2\
*q*(4*q + 3), -3*(2*q + 1), -(7*q + 2), -(5*q + 3), -4*(q + 1), -(4*q + 3)\
, q^2*(5*q + 3), -(5*q + 2), -3*(q + 1), 2*(2*q + 1), -2*(2*q + 1), -6*q*(\
2*q + 1)*(3*q + 1), 3*q + 2, -(3*q + 2), -2*q*(3*q + 1), 4*q + 1, -(4*q + \
1), 3*q + 1, -q*(5*q + 2), -(3*q + 1), 2*q + 1, -(2*q + 1), -(q + 1)]

                                  4, TRUE
Info: Building equations
Info: Building unknowns
Info: Solving system with 241 equations and 193 unknowns
Info: assumptions: [] [-q*(2*q + 1)*(7*q + 3), -(11*q + 13*q^2 + 2), (2*q \
+ 1)*(3*q + 1)*(q + 1)*(78*q + 89*q^2 + 16), (2*q + 1)*(5*q + 2), q^3*(13*\
q + 5)*(7*q + 3)*(4*q + 3), -2*(5*q + 2)*(3*q + 1), -(7*q + 9*q^2 + 2), -q\
*(q + 1)*(7*q + 3), -q*(5*q + 3)*(5*q + 2), -(7*q + 8*q^2 + 2), 4*q + q^2 \
+ 2, -(7*q + 7*q^2 + 2), (q + 1)*(5*q + 3), q*(q + 1), -q*(7*q + 3)*(7*q +\
 2)*(4*q + 1), -2*(q + 1)*(3*q + 1), -q*(7*q + 3)*(7*q + 2), (q + 1)^2, -3\
*(q + 1)^2, 2*(q + 1)*(2*q + 1), q, q^2*(q + 1)*(21*q + 29*q^2 + 4), 2*q^2\
, (q + 2)*(2*q + 1), (2*q + 1)^2, -4*(3*q + 1), -2*(7*q + 3)*(3*q + 1), 6*\
q*(2*q + 1), 3*(q + 1)*(5*q + 2), -3*q^2*(3*q + 2), q*(5*q + 3)*(51*q + 79\
*q^2 + 8), -3*q, q + 1, -q, -3*q^2*(3*q + 1), -(9*q + 2), q*(7*q + 3)*(3*q\
 + 1), -(7*q + 3), -3*(2*q + 1), -(7*q + 2), 4*(q + 1), -(5*q + 3), -2*(3*\
q + 1), 3*(q + 1), -q^2*(4*q + 3), -3*q*(5*q + 3)*(2*q + 1), -6*q*(2*q + 1\
)*(3*q + 1), 3*q*(2*q + 1), -2*(2*q + 1), -(7*q + 3)*(5*q + 2), -(5*q + 1)\
, 2*(2*q + 1), -3*(q + 1), 4*q + 1, 3*q + 2, 2*(5*q + 2)*(7*q + 3)*(q + 1)\
, -(3*q + 2), -(4*q + 1), -q*(5*q + 2), -2*(q + 1), 3*q + 1, -(3*q + 1), 2\
*q + 1, -q*(5*q + 3)*(2*q + 1)*(5*q + 2), -(2*q + 1), -(q + 1)]

                                  5, TRUE
Info: Building equations
Info: Building unknowns
Info: Solving system with 382 equations and 304 unknowns
Info: assumptions: [] [-q*(2*q + 1)*(7*q + 3), (5*q + 2)*(2*q + 1)*(5*q + \
3)*(208*q + 971*q^2 + 2318*q^3 + 2893*q^4 + 1548*q^5 + 18), (5*q + 1)*(4*q\
 + 1)*(q + 2), 11*q + 13*q^2 + 2, -q*(q + 1), q*(7*q + 3)*(2*q + 1)^2, -q*\
(3*q + 1)*(8*q + 3)*(q + 1), q^3*(5*q + 3)*(5*q + 2), -(7*q + 2)*(11*q + 1\
3*q^2 + 2), q*(8*q + 3)*(7*q + 2)*(q + 1), -(q + 1)*(8*q + 3), (q + 1)*(8*\
q + 3), (2*q + 1)*(5*q + 2), q^2, -q*(3*q + 2)*(11*q + 27*q^2 + 1), 2*(4*q\
 + 1)*(2*q + 1)*(3*q + 1), q^4*(4*q + 1)*(25*q + 37*q^2 + 4), -2*(5*q + 2)\
*(3*q + 1), -(7*q + 10*q^2 + 2), q^2*(8*q + 3)*(4*q + 1), 2*(2*q + 1)^2, -\
(8*q + 3)*(4*q + 1), (q + 1)*(7*q + 2), 8*q^2*(q + 2)*(2*q + 1)^3, -(7*q +\
 7*q^2 + 2), q*(q + 1), (3*q + 1)^2, q*(7*q + 3)*(7*q + 2), (q + 1)^2, -(2\
3*q + 78*q^2 + 84*q^3 + 2), -q*(7*q + 3)*(q + 1)^2, 3*q^2, (2*q + 1)*(3*q \
+ 1), 2*(q + 1)*(2*q + 1), q*(2*q + 1)*(13*q + 4), -2*q*(8*q + 3)*(3*q + 1\
), (q + 1)*(5*q + 1), -2*(q + 1)*(2*q + 1), q, -q*(3*q + 2)*(7*q + 3), -q*\
(6*q + 1)*(7*q + 3)*(7*q + 2), -(2*q + 3)*(q + 1), (q + 1)*(4*q + 1), -(6*\
q + q^2 + 2), -3*q^2*(2*q + 1)*(7*q + 2), -(q + 1)*(3*q + 1), -6*(q + 1)*(\
2*q + 1)^2, (q + 1)*(3*q + 1), -6*q, -q*(8*q + 3)*(5*q + 1)*(4*q + 1), (q \
+ 1)*(2*q + 1), 3*q*(2*q + 1)*(4*q + 1), -3*q, -(7*q + 2)*(q + 1)*(15*q + \
23*q^2 + 2), -2*q, -(7*q + 3)*(3*q + 1), -4*(2*q + 1), q + 1, -q, -(3*q + \
2)*(11*q + 4), -(8*q + 3), -(7*q + 3), -2*(4*q + 1), -3*(2*q + 1), -(26*q \
+ 44*q^2 + 3), -(5*q + 4), 7*q + 2, -(7*q + 2), -(5*q + 3), q*(8*q + 3), 2\
*(3*q + 1), q^2*(- 228*q - 646*q^2 - 125*q^3 + 1957*q^4 + 2276*q^5 + 336*q\
^6 - 24), -(5*q + 2), -2*q*(8*q + 3)*(4*q + 1), -(6*q + 1), 2*q*(4*q + 1),\
 q*(3*q + 1)*(q + 1)*(67*q + 188*q^2 + 353*q^3 + 571*q^4 + 468*q^5 + 9), -\
3*q*(5*q + 3)*(2*q + 1), -3*q*(2*q + 1), 3*q*(2*q + 1), (4*q + 1)*(15*q + \
15*q^2 + 4), -(7*q + 3)*(5*q + 2), -2*(4*q + 1)*(6*q + 1), q*(5*q + 2)*(4*\
q + 1), -(5*q + 1), 5*q + 1, -(3*q + 2), 3*q + 2, 2*(3*q + 1)^2, 4*q + 1, \
-(4*q + 1), -q*(4*q + 3), q*(4*q + 3), -q*(2*q + 1)*(8*q + 3), -2*(q + 1),\
 3*q + 1, -(3*q + 1), -3*q*(q + 1), 2*q + 1, -(2*q + 1), -2*(2*q + 1)*(8*q\
 + 3), -(11*q + 16*q^2 + 2), -(q + 1)]

                                  6, TRUE
Info: Building equations
Info: Building unknowns
Info: Solving system with 570 equations and 452 unknowns
Info: assumptions: [] [3*q^2*(2*q + 1)*(7*q + 2)*(56*q + 151*q^2 + 137*q^3\
 + 7), -2*q*(q + 1), -(2*q + 1)*(7*q + 2), -(13*q + 2)*(2*q + 1)*(9*q + 2)\
, -3*q*(2*q + 1)^2, -q*(3*q + 1)*(8*q + 3)*(q + 1), -3*(q + 1)*(3*q + 1), \
3*q*(2*q + 1)*(8*q + 3), -3*(3*q + 1)*(4*q + 1), q*(5*q + 3)*(2*q + 1), 2*\
q*(2*q + 1)*(7*q^2 - 1), (2*q + 1)*(5*q + 2), -3*q*(3*q + 1)*(q + 1)^2, -(\
3*q + 2)*(3*q + 1), -(7*q + 11*q^2 + 2), q*(q + 1)*(8*q + 3), -4*q*(4*q + \
1)*(2*q + 1)*(14*q + 5)*(7*q + 7*q^2 + 2)*(88534*q + 1976356*q^2 + 2707529\
0*q^3 + 255176557*q^4 + 1754410815*q^5 + 9088336148*q^6 + 36074264387*q^7 \
+ 110372144165*q^8 + 259538454514*q^9 + 463653314768*q^10 + 615708465954*q\
^11 + 586144599348*q^12 + 376354919916*q^13 + 145422695520*q^14 + 25527830\
400*q^15 + 1828), -(7*q + 10*q^2 + 2), -2*(5*q + 2)*(3*q + 1), -3*q*(3*q +\
 1)*(9*q + 2), -3*q*(4*q + 1)*(3*q + 1)*(q + 1), 2*(q + 1)*(14*q + 23*q^2 \
+ 2), -(7*q + 9*q^2 + 2), -q*(q + 1)*(7*q + 3), -2*(q + 1)*(3*q + 2)*(2*q \
+ 1)^2, q*(8*q + 3)*(q + 1)^2, -3*(7*q + 2), -6*q*(4*q + 1), 3*q*(3*q + 1)\
^3, q*(q + 1), -q*(11*q + 4)*(2*q + 1)*(7*q + 3), (3*q + 1)^2, -3*q^5, q*(\
5*q + 3)*(q + 1)*(- 111682*q - 1928013*q^2 - 19080566*q^3 - 119473998*q^4 \
- 486912666*q^5 - 1242460887*q^6 - 1559749664*q^7 + 1088688132*q^8 + 82815\
61892*q^9 + 15215926464*q^10 + 13630910256*q^11 + 5086126080*q^12 - 2808),\
 (2*q + 1)*(13*q + 4)*(4*q + 1)*(q + 1), (2*q + 1)*(4*q + 1), -6*(5*q + 1)\
*(7*q + 2)*(3*q + 1)*(6*q + 1)*(4*q + 1)^3, -(2*q + 1)*(4*q + 1), -q*(7*q \
+ 3)*(7*q + 2), -(q + 1)*(5*q + 2), (q + 1)^2, -6*q*(2*q + 1)*(4*q + 1), (\
q + 1)*(5*q + 2), q^3*(8*q + 9*q^2 + 2)*(6338*q + 98061*q^2 + 826620*q^3 +\
 4174980*q^4 + 12773772*q^5 + 21397473*q^6 + 8581492*q^7 - 35880732*q^8 - \
66058632*q^9 - 37013760*q^10 + 168)*(3*q + 1)^2, q*(3*q + 2)*(2*q + 1), -(\
4*q + 1)*(q + 1)*(17*q + 28*q^2 + 2), -2*q*(8*q + 3)*(2*q + 1)*(7*q + 2), \
2*(q + 1)*(2*q + 1), -(q + 1)*(5*q + 1), (q + 1)*(5*q + 1), q, -q*(2*q + 1\
)*(4*q + 1)*(7*q + 3)*(- 90*q - 463*q^2 - 933*q^3 - 542*q^4 + 72*q^5 - 6),\
 -q^2, (q + 1)*(3*q + 2), -3*q*(9*q + 2)*(3*q + 1)*(6*q + 1), -(2*q + 1)^2\
, 3*q*(3*q + 1)^2, q^2*(1492442*q + 38430021*q^2 + 610248632*q^3 + 6706812\
144*q^4 + 54219469253*q^5 + 334189265957*q^6 + 1604800326046*q^7 + 6077751\
854994*q^8 + 18245289248137*q^9 + 43364302654908*q^10 + 80999119024966*q^1\
1 + 117194161595108*q^12 + 128259569036712*q^13 + 102250546488408*q^14 + 5\
5806168221712*q^15 + 18590017000320*q^16 + 2851848000000*q^17 + 26840), (q\
 + 1)*(3*q + 1), -(q + 1)*(3*q + 1), (7*q + 2)*(4*q + 1), -(- 8*q + 11*q^2\
 + 50*q^3 - 2), -q*(7*q + 2)*(4*q + 1)*(2466*q + 46325*q^2 + 531334*q^3 + \
4188902*q^4 + 24148247*q^5 + 105188336*q^6 + 351015402*q^7 + 897342874*q^8\
 + 1738672553*q^9 + 2500799431*q^10 + 2580426502*q^11 + 1801691824*q^12 + \
759965760*q^13 + 145234944*q^14 + 60), -2*(8*q + 3)*(4*q + 1)^2, -5*q, (q \
+ 1)*(2*q + 1), -(q + 1)*(3*q + 1)^2, -(5*q + 1)*(7*q + 2)*(3*q + 1)*(q + \
2), -4*q, q*(7*q + 3)*(7*q + 1), 2*(21*q + 79*q^2 + 83*q^3 + 1)*(3*q + 1)^\
2, -3*q^2*(4*q + 1), -3*q, 2*q, q^2*(4*q + 1)*(q + 1)*(45210*q + 1219708*q\
^2 + 20506716*q^3 + 231913669*q^4 + 1845477564*q^5 + 10648619465*q^6 + 454\
01328711*q^7 + 144324575426*q^8 + 341771180961*q^9 + 595535359292*q^10 + 7\
42408306818*q^11 + 627009698376*q^12 + 321378192000*q^13 + 75410092800*q^1\
4 + 844), -(7*q + 3)*(3*q + 1), -2*q, -q*(2*q + 1)*(1244*q + 16852*q^2 + 1\
31357*q^3 + 653327*q^4 + 2163509*q^5 + 4801399*q^6 + 6943746*q^7 + 6035350\
*q^8 + 2558064*q^9 + 216000*q^10 + 40), -2*(11*q + 26*q^2 + 1), -(2*q + 1)\
*(- q + 20*q^2 - 2), -3*(3*q + 1), -2*(5*q + 1), q + 1, -q, q*(7*q + 3)*(3\
*q + 1), -(17*q + 17*q^2 + 4), 2*q*(q + 1)*(7*q + 3), -(8*q + 3), q*(3*q +\
 1), 2*q^2*(66506*q + 1746292*q^2 + 29142030*q^3 + 336551906*q^4 + 2809681\
546*q^5 + 17421758451*q^6 + 81625273075*q^7 + 291689113538*q^8 + 796574295\
795*q^9 + 1652283173137*q^10 + 2561137039190*q^11 + 2875964132886*q^12 + 2\
210852956608*q^13 + 1040834940480*q^14 + 226197619200*q^15 + 1260), -(7*q \
+ 3), -2*(4*q + 1), q*(5*q + 2)*(2*q + 1)^2, -3*(2*q + 1), 3*(2*q + 1), -(\
7*q + 2), -(3*q + 1)*(7*q + 2), -(5*q + 3), (q + 1)*(41*q + 316*q^2 + 1192\
*q^3 + 2284*q^4 + 1955*q^5 + 411*q^6 + 2), (q + 1)*(87*q + 202*q^2 + 149*q\
^3 + 12), q*(8*q + 3), 7*q + 1, -2*(3*q + 1), 2*(3*q + 1), q^2*(5*q + 3), \
5*q + 2, -2*q*(8*q + 3)*(4*q + 1), -3*q*(2*q + 1)*(3*q + 1), -2*q*(4*q + 1\
), -(6*q + 1), 6*q + 1, -(5*q + 2), 2*(3*q + 2)*(5*q + 2), 3*q*(2*q + 1), \
-3*(q + 1), -(2*q + 1)*(8*q + 3), -3*q*(3*q + 1)*(5*q + 1), 5*q + 1, -2*(2\
*q + 1), -(5*q + 1), q*(5*q + 3), 3*q + 2, 2*(3*q + 1)^2, 4*q + 1, (3*q + \
1)*(12*q + 19*q^2 + 2), -(4*q + 1), -q*(4*q + 3), -q*(7*q + 3)*(5*q + 2), \
-2*(q + 1), 2*(q + 1), -(3*q + 1), 36*q^3*(5*q + 2)*(2*q + 1)*(4*q + 1)*(3\
0*q + 139*q^2 + 198*q^3 + 2), 3*q^2*(7*q + 2)*(500*q + 3799*q^2 + 14898*q^\
3 + 33049*q^4 + 41944*q^5 + 27212*q^6 + 5328*q^7 + 26)*(2*q + 1)^2*(4*q + \
1)^2, -2*q*(2*q + 1), -q^2*(8*q + 3)*(5*q + 1), 2*q + 1, -(2*q + 1), -3*q^\
3*(3*q + 2)*(5*q + 2)*(2*q + 1)*(3*q + 1)*(371*q + 1020*q^2 + 819*q^3 + 40\
)*(38*q + 275*q^2 + 985*q^3 + 1832*q^4 + 1440*q^5 + 2), -(q + 1)]

                                  7, TRUE
Info: Building equations
Info: Building unknowns
Info: Solving system with 810 equations and 640 unknowns
\end{Mexout}

Testing that this element stabilizes the harmonics. The result is only
ok for 3 variables or less.
\begin{Mexin}
  vars := [x,y,z,t]:
  W:=Dom::WeylAlgebra_Table(vars, Dom::ExpressionField()):
  export(experimental::Steenrod):
  S:=experimental::SteenrodAlgebra():
  (P.i:=S::p::toWeyl(W, S::p([i]))) $i=1..3:
  P := Dom::DistributedPolynomial(vars, Dom::ExpressionField()):
  vdm := P(_mult((vars[i]-vars[j]) $j=1..i-1 $ i=1..nops(vars))):
  order := degree(vdm)+1:
  monoms := _concat(combinat::integerVectors(d, nops(vars),
                                             Outer=[$0..nops(vars)-1])
                                             $ d=0..degree(vdm)):
  D1 := W::adjoint(P1):
  D2 := W::adjoint(P2):
  harmonics := []:
  for m in monoms do
      p := W::onPoly(W::adjoint(SteenrodConjugated(W,order,m)), vdm);
      if not iszero(D1(p)) or not iszero(D2(p)) then
          error("result is not harmonic");
      end_if;
      harmonics := append(harmonics,p);
  end_for;
  solve([coeff(mapcoeffs(_plus(hold(a).i*harmonics[i] $i=1..nops(harmonics))
                          normal))],
         [ hold(a).i $ i=1..nops(harmonics) ]);
\end{Mexin} %$

Computation for 4 variables; damned, we just have one degree of
freedom:

res2 := table(op(select(res, eq->bool(_plus(op(op(eq,1)))<order+1))));
res3:=map(res2, Pq@subs, f(0,0,0,0)=1, f(0,1,0,0)=a);

The solution can be described as follow:

mapcoeffs(F, subs, f= (d) ->  (d*(a-1)*q + _plus(args()) *(a-1) +1));

Problem: it kills almost everything, whatever choice of a we do:

F2(P(1)), F2(P(x)), F2(P(x^2)), F2(P(x*y))

Attempt 3: we do not impose that F is divisible by x
\begin{Mexin}
reset(): export(experimental::Steenrod):debugOff():
vars := [x,y,z,t]: order := 5: q:=3:
W:=Dom::WeylAlgebra_Table(vars, Dom::ExpressionField()):
P:=Dom::DistributedPolynomial(vars, Dom::ExpressionField()):
generic:=experimental::Steenrod::genericWeylElementWithSymmetries:
setuserinfo(experimental::Steenrod,1): setuserinfo(linalg,2):
makeUnknownsF :=
proc(l)
    option remember;
    local l2, i;
begin
    l2 := op(l,2);
    for i from nops(l2) downto 1 do if l2[i]>0 then break; end; end;
    l2 := [op(l2, 1..i)];
    if op(l,[1,1])>0 then
        f1(op(l2))
    else
        f2(op(l2))
    end;
end_proc:
makeCoeffsF   :=
proc(m)
    local k, n;
begin
    m2:=op(m,2);
    n := _plus(op(m2));
    1/_mult(k! $ k in m2)*makeUnknownsF(m);
end_proc:
F :=_plus(generic(W, 1,o,[1,nops(vars)-1],makeCoeffsF)
          $ o=0..order+1):
unknownsF := sort(valuesOfRememberTable(makeUnknownsF),
                  (x,y)->if op(x,0)<>op(y,0) then
			     bool(op(x,0)=f2)
			 else
			     combinat::words::degInvLexLess([op(y)],[op(x)]);
			 end_if):
time((res := solveInSteenrodRightIdeal(W, F, unknownsF, order, Q=q))):
res2 := map(table(op(select(res, eq->bool(_plus(op(op(eq,1)))<order+1)))),
            factor);

F2:=mapcoeffs(truncateOrder(F, order), subs, res2):
print(k, factor(F2(P(x^k)))) $k=0..order; %$
p:=P(x): factor(F2(p))/factor(P1(p));
\end{Mexin}

Conclusion: without restriction on $f_\mu$, we can always find a
Steenrod conjugated of the form $F:=x_1 \sum_\mu \frac1{\mu!} f_\mu +
x_2\dots$ (probably simply take $\P_1 / (1+qt_1) \sum_ \frac1{\mu!}
f_\mu$). The image of $F$ on a monomial $x^m$ of degree at least $2$
is then the following multiple of $\P_1(x^m)$:

\begin{displaymath}
F2(x^m)=\frac{1}{1+q\max(m_i)}
(\sum_{k_1=1}^{m_1} \dots \sum_{k_n=1}^{m_n}
\binom{m_1}{k_1}\dots\binom{m_1}{k_m} f_{k_1,\dots,k_n}) \P_1(x^m)
\end{displaymath}

Checking that Sq1 is a solution:
\begin{Mexin}
myf1 := ()->0:myf1():=1: myf1(1):=-q:
myf2 := ()->0:myf2():=1: myf2(0,1):=-q:
map(res2, normal@subs, f1=myf1, f2=myf2);
\end{Mexin}

\subsubsubsection{Expressing the non-symmetric degree 1 Steenrod
  conjugated using the wedge product}

\begin{Mexin}
fPoly :=
proc(P: DOM_DOMAIN, order: Type::NonNegInt)
begin
    P([[fCoeff(m), m]
       $ m in combinat::integerVectors(d, nops(polylib::variables(P)))
       $ d = 0..order-1]);
end_proc:
Exp:=
proc(P: DOM_DOMAIN, order: Type::NonNegInt)
begin
    P([[(-1)^_plus(op(m))/_mult(k! $k in m), m]
       $ m in combinat::integerVectors(d, nops(polylib::variables(P)))
       $ d = 0..order-1]);
end_proc:
\end{Mexin} %$

One variable:
\begin{Mexin}
export(experimental::Steenrod):
order := 7:
P := Dom::DistributedPolynomial([t], Dom::ExpressionField()):
p:=Exp(P,order):
solve([coeff(truncate((a+b*t+c*t^2)*p, order)-fPoly(P,order))],
      [a,b,c])
bool(P(_plus((-1)^k/k!*     t^k $k=0..order-1)) = truncate(p,         order));
bool(P(_plus((-1)^k/k!*k   *t^k $k=0..order-1)) = truncate(-p*t,      order));
bool(P(_plus((-1)^k/k!*k^2 *t^k $k=0..order-1)) = truncate(p*t*(t-1), order));
bool(fPoly(P, order) = truncate( ((q+1/2) * t^2 + (q+1)*t +1) * p, order));
\end{Mexin} %$

Several variables:
\begin{Mexin}
order := 7: vars := [t1,t2,t3]:
P := Dom::DistributedPolynomial(vars, Dom::ExpressionField()):
S:=examples::SymmetricFunctions():
p:=Exp(P,order):
bool(truncate(((1+S::e([1])+1/2*S::e([1,1])) (vars) + q*t1*(1+S::e([1]))(vars))*p, order) = fPoly(P,order))
\end{Mexin}

\begin{conjecture}
  Set $t_i:=x_i \p_i$ and $e:=t_1+\dots+t_n$. Then, for the wedge
  product, we have:
  \begin{displaymath}
    f = x_1 (1+e+\frac12 e^2 + q t_1 (1+e)) \exp(-e)
  \end{displaymath}
\end{conjecture}

\begin{Mexin}
ExpWeyl :=
proc(W: DOM_DOMAIN, order: Type::NonNegInt)
begin
    _plus((-1)^_plus(op(m))/_mult(k! $k in m) *
          W::term(_concat(zip([0,0 $ W::dimen-1],m,_plus),m))
	  $ m in combinat::integerVectors(d, W::dimen)
	  $ d = 0..order-1)
end_proc:
ExpWeylInv :=
proc(W: DOM_DOMAIN, order: Type::NonNegInt)
begin
    _plus(1/_mult(k! $k in m) *
          W::term(_concat(zip([0,0 $ W::dimen-1],m,_plus),m))
	  $ m in combinat::integerVectors(d, W::dimen)
	  $ d = 0..order-1)
end_proc:
\end{Mexin}

The expression of f looks very similar to ExpWeyl.
Looking for the inverse of ExpWeyl. That's stupid: it cannot exists,
since ExpWeyl is a zero diviser (it almost kills everything).
\begin{Mexin}
order := 4:
W:=Dom::WeylAlgebra_Table([x], Dom::ExpressionField()):
makeUnknownsF := proc(l) option remember; begin hold(f)(op(l,[2,1])) end:
F := _plus(generic(W, 0,o,[1],makeUnknownsF) $ o=0..order-1);
unknowns := sort(valuesOfRememberTable(makeUnknownsF),
                 (x,y)->combinat::words::degLexLess([op(y)],[op(x)])):
equations := [ coeff(truncateOrder(F*ExpWeyl(W,order)-W::one, order)) ];
linsolve(equations, unknowns);
\end{Mexin} %$

\subsection{Looking for diagonal Steenrod conjugated operator, as
  linear operators}

The approach used here is much more efficient, and yields equations
that are much simpler to analyse. Key idea: an operator which is a
series in the euler operators acts diagonally on monomials (that is
sends a monomial to a scalar multiple of this monomial); reciprocally,
any such diagonal operator can be expressed as a series in the euler
operators.

So, in the expressions f = \sum x_i f_i, one may as well implement
each f_i as a linear functions that sends a monomial $x^m$ to
$f_i(m)x^m$.

The equations are much simpler: each coefficient f_i(m) will appear
only in a few of them (the equations in the f_i's are graded wrt the
degree of m). Typically $\sum (-1)^\mu/\mu! t^\mu$ is replaced by the
function which zero on every monomial except $1$; the coefficient of
this function won't propagate in each degree.

To get simpler coefficients, we assume that $(1+qt_i)$ is invertible,
(that is $q$ is not of the form $-\frac1d$), and write the searched
operator as $f:=\sum x_i (1+qt_i) f_i$. We do similarwise for the
operators $g$, $h1$ and $h2$ on the right side.

\begin{Mexin}
reset(): export(experimental::Steenrod)://debugOff():
vars := [x,y,z]: order := 4: //q:=1:
P  := Dom::DistributedPolynomial(vars, Dom::ExpressionField()):
makeUnknown :=
proc(identifier: DOM_IDENT)
    option escape;
begin
    proc(m: combinat::integerVectors, i: Type::PosInt)
	option remember;
    begin
	(identifier.i)[op(m)]
    end;
end_proc:
makeF :=
proc(f: DOM_PROC)
    local i;
    option escape;
begin
    operators::makeLinear
      (m -> _plus(f(m,i)*P([[1,m]]) $ i=1..nops(vars)),
       Source=P, ImageSet=P):
end_proc:
P1 := makeF((m,i)->(1+q*op(m,i))*P(vars[i])):
P2 := makeF((m,i)->(1+q*op(m,i))*P(vars[i])^2):
FU := makeUnknown(hold(f)):
GU := makeUnknown(hold(g)):
H1U:= makeUnknown(hold(h1)):
H2U:= makeUnknown(hold(h2)):
F  := makeF((m,i)->FU (m,i)*(1+q*op(m,i))*P(vars[i])):
// Remark: if $f P1 = P1 g$, then one can easily prove that
// g1([m1,m2,m3]) = f1([m1+1,m2,m3])
// G  := makeF((m,i)->FU (m,i)*(1+q*op(m,i))*P(vars[i])):
G  := makeF((m,i)->FU ((m2:=m;m2[i]:=m2[i]+1;m2),i)*(1+q*op(m,i))*P(vars[i])):
H1 := makeF((m,i)->H1U(m,i)*P(vars[i])):
H2 := makeF((m,i)->H2U(m,i)*(1+q*op(m,i))*P(vars[i])^2):
equations :=
  [ (coeff(F(P1(m))-P1(G(m))), coeff(F(P2(m))-P1(H2(m))-P2(H1(m))))
    $ m in map(combinat::integerVectors::list(d, nops(vars)), m->P([[1,m]]))
    $ d = 0..order-1 ]:
unknowns :=
  [ op(valuesOfRememberTable(x)) $ x in [ H1U, H2U, FU ]]:
result := linsolve(equations, unknowns):
res2 := table(op(select(result, not has, [(g.i, h1.i, h2.i) $i=1..nops(vars)])));
\end{Mexin} %$

Conclusion:
\begin{gather}
H1=0\\
g_i ([m_1,\dots,m_i,\dots,m_n]) = f_i([m_1,\dots,m_i+1,\dots,m_n])\\
h2_i([m_1,\dots,m_i,\dots,m_n]) = f_i([m_1,\dots,m_i+2,\dots,m_n])\\
f_1(m)=\dots=f_n(m)\\
\end{gather}

Not that f has conjugates w.r.t. both P1 and P2!

\begin{conjecture}
  Inside $\K[x_i][[t_i]]$, the Steenrod-conjugated operators of degree
  $1$ are of the form $P_1.f + g$, where $f\in \K[[t_i]]$, and $g$ is
  a degree one operator which kills any monomial of degree $\geq2$.
\end{conjecture}
In short, considering such Steenrod conjugated operators does not
yield new interesting operators that stabilize the hits.




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Embedding the dual of the steenrod algebra inside
  quasi-symmetric functions}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{Mexin}
reset():
S:=experimental::SteenrodAlgebra():
QSym:=examples::QuasiSymmetricFunctions():
alias(Sp=S::powersum::domainWrapper):
alias(Dp=S::powersum::dual::domainWrapper):
export(operators):

embed:=
proc(d: Type::NonNegInt, algebra = QSym::m : DOM_DOMAIN)
    option remember;
    local p,i;//, vars, coprod, result;
begin
    // Compute the embeding of the algebra generators of lower degree,
    // and define the embeding of the generator of degree d as a generic
    // element of algebra
    p:=[embed(i, algebra) $i=1..d-1, algebra::genericElement(hold(a)@op, d)];
    vars := revert([hold(a)(op(p)) $ p in combinat::compositions(d)]);
    // Build by recurence the coproduct of Dp[d], as an element of algebra
    coprod := _plus(op(t,1) * tensor(algebra::mult(p[i] $ i in part) $ part in op(t,2))
                    $ t in poly2list(coproduct(Dp[d])));
    // Identify the coproduct of p[d] with coprod, and solve the system
    eqns := map(select([coeff(coproduct(p[d])-coprod)], not iszero), normal);
    result := linsolve(eqns, vars);
    mapcoeffs(p[d], subs, result);
end:
\end{Mexin} %$

(Note: for d=6, MuPAD finds a solution with solve, but not with
linsolve if one does not map normal on the coefficients)

The most natural ?
\begin{Mexin}
mapcoeffs(embed(6), factor@subs, (f(i)=i,a(i)=1) $i=1..6);
\end{Mexin}
This one has no denominators:
\begin{Mexin}
mapcoeffs(embed(6), factor@subs, (f(i)=a(i),a(i)=1) $i=1..6);
\end{Mexin}

\section{Looking for descending operators}

Attempt 1:

bracket :=(x,y) -> x*y-y*x:
Pk := k-> _plus(W(x^k)*W(x*D(x))    $ x in W::XVariables):
Dk := k-> _plus(W(x*D(x))*W(D(x))^k $ x in W::XVariables):
toEuler := // Tries to express a weyl element as a polynomial in the
           // euler operator
proc(p)
    local res;
begin
    res := solve([coeff(p - _plus(a[k]*euler^k $k=0..W::order(p)))]);
    if nops(res)=0 then FAIL else
       subs(_plus(a[k]*hold(e)^k $k=0..W::order(p)), op(res))
    end_if;
end;


\section{Computing the leading terms of the hits}

hitsLTerms :=
proc(n, d)
begin
vars := [hold(x).i $i=1..n]:
R := Dom::ExpressionField():
W := Dom::WeylAlgebra_Table(vars, R):
S := (experimental::SteenrodAlgebra(R))::p:
P := examples::MultivariatePolynomials(vars, R):
(hold(P).i :=S::toWeyl(W, S([i]))) $i=1..3:
act := (s, p) -> P::fromExpr((W::onPoly(s,expr(p)))):

P::basis(d): // P::tterm(p) is the lexicographically largest monomial of p!!!
M := _concat(Dom::Matrix(R)(0,0),
	     if d>=1 then P::sparseVector(act(P1, p)) $ p in P::basis(d-1) else null(); end_if,
             if d>=2 then P::sparseVector(act(P2, p)) $ p in P::basis(d-2) else null(); end_if):
map(sort([op(op(M::dom::gaussElim(M, ColumnElimination), 4))]), P@P::unrankBasis):
Mg := op(M::dom::gaussElim(M, ColumnElimination),1):

// A triangular basis of the image
basis  := select([P::fromSparseVector(Mg::dom::col(Mg,i)) $ i=1..Mg::dom::matdim(Mg)[2]],
                 not iszero):
return(map(basis, P::tterm));
end_proc:

hitsLTermsGen :=
proc(n,dmax)
    local lt, lv, d, t, v, v2, isNew;
begin
    lt := [];
    lv := [];
    for d from 1 to dmax do
        for t in hitsLTerms(n, d) do
            v := t::dom::lsupport(t);
	    isNew := TRUE;
	    for v2 in lv do
		if min(op(zip(v,v2,_subtract))) >= 0 then
		    isNew := FALSE;
		    break;
		end_if;
            end_for;
	    if isNew then
		lt := append(lt,t);
		lv := append(lv,v);
	    end_if;
	end_for;
    end_for;
    lt;
end_proc:

%$

\section{}

export(experimental::Steenrod):
setuserinfo(spair,0):

    spair2 :=
    proc(i,j,m)
	local q;
    begin
	q := p[j]*m;
	r1[i] := mapcoeffs(r1[i]/lcoeff(p[i])-r1[j]/lcoeff(q), factor);
	r2[i] := mapcoeffs(r2[i]/lcoeff(p[i])-r2[j]/lcoeff(q), factor);
	p [i] := mapcoeffs(p [i]/lcoeff(p[i])-q    /lcoeff(q), factor):
    end_proc;

vars := [X,Y,Z,T]: tvars := [tx,ty,tz,tt]: maxDegree := 6:
Q:=PolynomialsWithSkewThetaCoefficients(vars, tvars, Dom::ExpressionField()):
// Ensures that the lexicographic term order is used.
for d from 1 to maxDegree do
    map(revert(Q::basisIndices(d)), Q::term):
end_for:
[x,y,z,t] := [Q([0$i-1,1,0$nops(vars)-i]) $i=1..nops(vars)]:
p  := table(i=experimental::Steenrod::qPk(Q,i) $i=1..2):
r1 := table(1=Q::one,2=Q::zero):
r2 := table(1=Q::zero,2=Q::one):
for v in vars do
    spair2(2,1,Q(v)):
end_for:

set(2,2,(q*tx+1)*(q*tt+1)):
// p[2] est symmétrique par rapport aux variables 2,3,4,...
prog::test(expr(p[2]) - subs(expr(p[2]), [Z=T, T=Z, tz=tt, tt=tz]), 0):
prog::test(expr(p[2]) - subs(expr(p[2]), [Z=Y, Y=Z, tz=ty, ty=tz]), 0):

set(3,2,x):
for i from 2 to nops(vars) do
    for j from i to nops(vars) do
	spair2(3, 1, Q(vars[i])*Q(vars[j])):
    end_for;
end_for:
for i from 2 to nops(vars) do
    spair2(3, 2, Q(vars[i])):
end_for:
// p[3] est symmétrique par rapport aux variables 3,4,...
set(3,3,(q*tt+1)*(q*ty+1)*(q + q*tx + q*ty + 2)*(2*q + q*tx + q*ty + 2)/q):
prog::test(a - subs(a, [Z=T, T=Z, tz=tt, tt=tz]),0):

bilan() ;

set(4,3,y):
for i from 3 to nops(vars) do
    for j from i to nops(vars) do
	spair2(4, 2, Q(vars[i])*Q(vars[j])):
    end_for;
end_for:
// En fait, c'est à regrouper avec le p[3]*y ci-dessus;
// mais l'ordre sur les termes inverse l'ordre des opérations
for i from 3 to nops(vars) do
    spair2(4,3,Q(vars[i])):
end_for;

bilan() ;


Variation: Note: seems to require MuPAD 3.0.0 due to a bug in
Factored? (maybe _concat?) when computing q*Q(tx)

export(experimental::Steenrod):
vars := [X,Y,Z,T]: maxDegree := 6:
K := Dom::ExpressionField(normal):
K::Name := "K":
operators::overloaded::declareConversion(Dom::ExpressionField(),
                                         K, K):
Q:=PolynomialsWithSkewThetaCoefficients(vars, K):
// Ensures that the lexicographic term order is used.
for d from 1 to maxDegree do
    map(revert(Q::basisIndices(d)), Q::term):
end_for:
[x,y,z,t] := [Q([0$i-1,1,0$nops(vars)-i]) $i=1..nops(vars)]:
p  := table(i=experimental::Steenrod::qPk(Q,i) $i=1..2):
setuserinfo(Any,1):
p[2] := Q::remOneStep  (p[2],p[1]):
p[2] := Q::remSymmetric(p[2],p[1]):
p[2] := p[2] * (1+q+q*tX):
//Info: factor: T*(-(q*tT + 1)/(q + q*tX + 1)) + Y*(-(q*tY + 1)/(q + q*tX + 1)) + Z*(-(q*tZ + 1)/(q + q*tX + 1))
// p[2] est symmétrique par rapport aux variables 2,3,4,...
prog::test(expr(p[2]) - subs(expr(p[2]), [Z=T, T=Z, tZ=tT, tT=tZ]), 0):
prog::test(expr(p[2]) - subs(expr(p[2]), [Z=Y, Y=Z, tZ=tY, tY=tZ]), 0):

p[3] := p[2]*x:
p[3] := Q::remSymmetric(p[3], p[1]):
// Info: factor: Y^2 c + ...
// Info: factor: Y*Z c + ...
p[3] := Q::remOneStep(p[3], p[2]):
// Info: factor: Y c
p[3] := Q::remSymmetric(p[3], p[2]):
// Info: factor: Z c + ...
p[3] :=  p[3] * (1+q*tX)  *(2*q + 1+q*tX) *  (2*q + 1+q*tX + 1+q*tY) / q:
// p[3] est symmétrique par rapport aux variables 3,4,...
prog::test(a - subs(a, [Z=T, T=Z, tZ=tT, tT=tZ]),0):

p[4] := p[3]*y:
p[4] := Q::remSymmetric(p[4],p[2]):
p[4] := Q::remOneStep(p[4],p[3]):
p[4] := Q::remSymmetric(p[4],p[3]):






export(experimental::Steenrod):
vars := [X,Y,Z,T]: maxDegree := 6:
K := Dom::ExpressionField(normal):
K::Name := "K":
operators::overloaded::declareConversion(Dom::ExpressionField(),
                                         K, K):
Q:=PolynomialsWithSkewQThetaCoefficients(vars, K):
[x,y,z,t] := [Q([0$i-1,1,0$nops(vars)-i]) $i=1..nops(vars)]:
p  := table(i=experimental::Steenrod::qPk(Q,i) $i=1..2):
q:=1;
setuserinfo(Any,1):
p[2] := Q::remOneStep  (p[2],p[1]):
p[2] := Q::remSymmetric(p[2],p[1]):
p[2] := p[2] * (q+tX):
//Info: factor: T*(-(q*tT + 1)/(q + q*tX + 1)) + Y*(-(q*tY + 1)/(q + q*tX + 1)) + Z*(-(q*tZ + 1)/(q + q*tX + 1))
// p[2] est symmétrique par rapport aux variables 2,3,4,...
prog::test(expr(p[2]) - subs(expr(p[2]), [Z=T, T=Z, tZ=tT, tT=tZ]), 0):
prog::test(expr(p[2]) - subs(expr(p[2]), [Z=Y, Y=Z, tZ=tY, tY=tZ]), 0):

p[3] := p[2]*x:
p[3] := Q::remSymmetric(p[3], p[1]):
// Info: factor: Y^2 c + ...
// Info: factor: Y*Z c + ...
p[3] := Q::remOneStep(p[3], p[2]):
// Info: factor: Y c
p[3] := Q::remSymmetric(p[3], p[2]):
// Info: factor: Z c + ...
p[3] :=  p[3] * tX*  (2*q + tX + tY) / q:
// p[3] est symmétrique par rapport aux variables 3,4,...
prog::test(a - subs(a, [Z=T, T=Z, tZ=tT, tT=tZ]),0):

p[4] := p[3]*y:
p[4] := Q::remSymmetric(p[4],p[2]):
p[4] := Q::remOneStep(p[4],p[3]):
p[4] := Q::remSymmetric(p[4],p[3]):


p[5] := p[3]*x;
Q::rem(p[5],p[1]) ;
Info: factor: (Y*Z^2)*(1/(tX + tX*tY)*(tY*tZ^2 - 3*tY^2*tZ - tY^3*tZ - tY^2*tZ^2 + 4*tY*tZ + tX\
*tY*tZ^2 + tX^2*tY*tZ + 5*tX*tY*tZ))
Info: factor: (T*Y*Z)*(1/(tX + tX*tY)*(2*tT*tY*tZ - 2*tT*tY^2*tZ + 2*tT*tX*tY*tZ))
Info: factor: (T^2*Y)*(1/(tX + tX*tY)*(tT^2*tY - 3*tT*tY^2 - tT*tY^3 - tT^2*tY^2 + 4*tT*tY + tT\
*tX^2*tY + tT^2*tX*tY + 5*tT*tX*tY))
Info: factor: Z^3*(-1/tX*(16*tZ + 2*tX*tZ^2 + tX^2*tZ + tY*tZ^2 + 8*tX*tZ + 4*tY*tZ + 8*tZ^2 + \
tZ^3 + tX*tY*tZ))
Info: factor: (T*Z^2)*(-1/tX*(3*tT*tZ^2 + 8*tT*tZ + 2*tT*tX*tZ + tT*tY*tZ))
Info: factor: (T^2*Z)*(-1/tX*(3*tT^2*tZ + 8*tT*tZ + 2*tT*tX*tZ + tT*tY*tZ))
Info: factor: T^3*(-1/tX*(16*tT + tT*tX^2 + 2*tT^2*tX + tT^2*tY + 8*tT*tX + 4*tT*tY + 8*tT^2 + \
tT^3 + tT*tX*tY))
p[5]:=Q::rem(p[5],p[1]):
Info: factor: (Y*Z^2)*(1/(tX + tX*tY)*(tY*tZ^2 - 3*tY^2*tZ - tY^3*tZ - tY^2*tZ^2 + 4*tY*tZ + tX\
*tY*tZ^2 + tX^2*tY*tZ + 5*tX*tY*tZ))
Info: factor: (T*Y*Z)*(1/(tX + tX*tY)*(2*tT*tY*tZ - 2*tT*tY^2*tZ + 2*tT*tX*tY*tZ))
Info: factor: (T^2*Y)*(1/(tX + tX*tY)*(tT^2*tY - 3*tT*tY^2 - tT*tY^3 - tT^2*tY^2 + 4*tT*tY + tT\
*tX^2*tY + tT^2*tX*tY + 5*tT*tX*tY))
Info: factor: Z^3*(-1/tX*(16*tZ + 2*tX*tZ^2 + tX^2*tZ + tY*tZ^2 + 8*tX*tZ + 4*tY*tZ + 8*tZ^2 + \
tZ^3 + tX*tY*tZ))
Info: factor: (T*Z^2)*(-1/tX*(3*tT*tZ^2 + 8*tT*tZ + 2*tT*tX*tZ + tT*tY*tZ))
Info: factor: (T^2*Z)*(-1/tX*(3*tT^2*tZ + 8*tT*tZ + 2*tT*tX*tZ + tT*tY*tZ))
Info: factor: T^3*(-1/tX*(16*tT + tT*tX^2 + 2*tT^2*tX + tT^2*tY + 8*tT*tX + 4*tT*tY + 8*tT^2 + \
tT^3 + tT*tX*tY))
p[5]:=Q::rem(p[5],p[2]):
Info: factor: Z^2*(-1/(2*tX + tX*tY + tX^2)*(4*tZ + tX*tZ^2 + tX^2*tZ - tY*tZ^2 - tY^2*tZ + 5*t\
X*tZ - 3*tY*tZ + tZ^2))
Info: factor: (T*Z)*(-(2*tT*tZ + 2*tT*tX*tZ - 2*tT*tY*tZ)/(2*tX + tX*tY + tX^2))
Info: factor: T^2*(-1/(2*tX + tX*tY + tX^2)*(4*tT + tT*tX^2 + tT^2*tX - tT*tY^2 - tT^2*tY + 5*t\
T*tX - 3*tT*tY + tT^2))
p[5]:=Q::rem(p[5],p[3]):
Info: factor: Z*(1/(2*tX^2*tZ - tX*tY^2 + tX^3*tZ - 2*tX*tY + 2*tX^2 + tX^3 - tX*tY^2*tZ - 2*tX\
*tY*tZ)*(8*tZ - tX^2*tZ + 4*tY*tZ^2 + 2*tY^2*tZ - tX^2*tZ^2 + 4*tX*tZ + 14*tY*tZ + 2*tZ^2 + 3*t\
X*tY*tZ^2 + tX*tY^2*tZ + tX^2*tY*tZ + 11*tX*tY*tZ))
Info: factor: T*(1/(tX^3*tY - tX*tY^3 - 5*tX*tY^2 + 2*tX^2*tZ + tX^3*tZ - tX^2*tY^2 - 6*tX*tY +\
 6*tX^2 + 5*tX^3 + tX^4 - tX*tY^2*tZ - 2*tX*tY*tZ)*(2*tT*tY^2 - tT*tX^3 - 3*tT*tX^2 - 2*tT*tX +\
 4*tT*tY + 2*tT*tZ + 2*tT*tX*tY^2 + tT*tX^2*tY - tT*tX^2*tZ + 5*tT*tX*tY + 4*tT*tY*tZ + 3*tT*tX\
*tY*tZ))
p[5]:=Q::rem(p[5],p[4]):
Info: factor: 1*(-1/(6*tX + 4*tX*tY^2 + 18*tX^2*tY + 3*tX^3*tY + 9*tX^2*tY^2 + tX^2*tY^3 + tX^3\
*tY^2 + 10*tX*tY + 10*tX^2 + 2*tX^3)*(12*tY + 20*tX*tY^2 + 10*tX*tY^3 + tX*tY^4 + 2*tX^2*tY^2 +\
 tX^2*tY^3 + 8*tX*tY + 26*tY^2 + 14*tY^3 + 2*tY^4))
>> bilan() ;
1	X
2	Y^2
3	Y*Z^2
4	T^2*Y*Z
5	Z^4



export(experimental::Steenrod):
vars := [X,Y,Z,T]:
K := Dom::ExpressionField(normal):
K::Name := "K":
operators::overloaded::declareConversion(Dom::ExpressionField(),
                                         K, K):
Q:=PolynomialsWithSkewQThetaCoefficients(vars, K, 1):
//[x,y,z,t] := [Q([0$i-1,1,0$nops(vars)-i]) $i=1..nops(vars)]:
sys  := [ experimental::Steenrod::Pk(Q,i) $i=1..2 ]:
setuserinfo(Any,1):
debugOf():
Q::groebner(sys);




\begin{remark}
  $x$ is the dual of $\partial$ on polynomials, but not on Laurent
  polynomials (unless $<1/x|1/x>=\infty$):

  $1 = <1|1> = <x 1/x | 1> \ne <1/x | \partial 1> = 0 <1/x | 1/x>$,
\end{remark}

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: t
%%% ispell-local-dictionary: american
%%% End: 
