# STEENROD
# A Maple Package for Computing with the Steenrod Algebra  
> 
# 
# by:  Ken Monks
#      Department of Mathematics
#      University of Scranton
#      Scranton, PA  18510
#      (717) 941-6101
#      monks@scranton.edu
# 
# MBasis routine by: 
#      Vince Giambalvo 
#      University of Connecticut 
#      Storrs, CT 06269
#      vince@math.uconn.edu
#  
# Version 10.3
# [Last Revised - March 2007 - made the pretty printers global]
# Source Code:
> restart;
> 
> Steenrod:=module()
> 
> local MakeProc;
> 
> global Sq,xi,T,Q,`print/T`,`print/&x`;
> 
> export init,IsInAn,MBasis,`&rlex`,SqDickson,tri,V,v,
>        NULLSpace,ReduceToBasis,ExtendToBasis,FindPivotColumns,
>        Equal,Inverse,ZERO,id,delta,w,s1,s2,omega,nu,alpha,
>        Alpha,multinom,Parts,IntSolve,Top,Bottom,Nil,ChangeBasisMatrix,
>        AdmisAnBasis,AnBasis,Filtration,May,DEGREE,Excess,
>        K,ActOn,chi,Coproduct,`&^`,`&x`,`&xx`,`&*`,MilnorMult,
>       
> DualMult,TensorMilnorMult,TensorDualMult,Gamma,Flip,BinaryChart,
>        Adem;
> 
> init:=proc()
>    local MakeProc,Zero,One,ListAtoms;
>    global
> Sq,xi,T,Q,`type/MilnorBasis`,`type/Milnor`,`type/Dual`,`type/P_s`,
>        `type/Inert`,`type/MonomialSqn`,`type/List`,`type/TensorDual`,
>        `type/TensorMilnor`,`type/AdmissibleBasis`,`type/Admissible`,
>       
> `convert/List`,`convert/Dual`,`convert/P_s`,`convert/Admissible`,
>        `convert/Milnor`,`convert/Inert`,`print/T`,`print/&x`;
>    
>    print("Steenrod package: by Ken Monks, Version 10.3"):
>    print("MBasis routine by Vince Giambalvo"):
>    printf("\n");
> 
>     Zero:=identical(0);
>     One:=identical(1);
>    `type/MilnorBasis`:={One,specfunc(nonnegint,Sq)}:
>    `type/Milnor`:={Zero,MilnorBasis,''`+`''(MilnorBasis)}:
>   
> `type/Dual`:={Zero,One,xi[posint],''`*`''(Dual),''`+`''(Dual),Dual^pos
> int}:
>   
> `type/P_s`:={Zero,One,X[posint],''`*`''(P_s),''`+`''(P_s),P_s^posint}:
>    `type/Inert`:={Milnor,specfunc(Inert,`&x`),''`+`''(Inert)}:
>   
> `type/MonomialSqn`:={Zero,One,''Sq''(posint),specfunc(MonomialSqn,`&x`
> )}:
>     ListAtoms:={Zero,One,list(nonnegint)};
>   
> `type/List`:={ListAtoms,specfunc(ListAtoms,T),list(specfunc(ListAtoms,
> T)),list(ListAtoms)}: 
>   
> `type/TensorDual`:={Dual,specfunc(TensorDual,T),''`+`''(TensorDual)}:
>   
> `type/TensorMilnor`:={Milnor,specfunc(TensorMilnor,T),''`+`''(TensorMi
> lnor)}:
>    `type/AdmissibleBasis`:=proc(a)
>       local i,A;
>       if not type(a,MonomialSqn) then RETURN(false) fi;
>       A:=[op(convert(a,List))];
>       for i to nops(A)-1 do 
>         if A[i]<2*A[i+1] then RETURN(false) fi;
>       od;
>       true;
>     end:
>    `type/Admissible`:={Zero,AdmissibleBasis,''`+`''(AdmissibleBasis)}:
> 
>    `convert/List`:=proc(a)
>       local m,ans,i,j;
> 
>       # don't convert Lists
>       if type(a,List) then RETURN(a)
>       
>       # convert sums to lists
>       elif type(a,`+`) then RETURN(map(`convert/List`,[op(a)]))
> 
>       # Monomials in Sq(n)
>       elif type(a,MonomialSqn) then
> RETURN([seq(op(op(i,a)),i=1..nops(a))]);
> 
>       # Milnor
>       elif type(a,Milnor) then RETURN([op(a)])
>       
>       # Dual or P_s
>       elif type(a,{Dual,P_s}) then
>         if nargs>1 then m:=max(args[2],max(op(map(op,indets(a)))))
>         else m:=max(op(map(op,indets(a))))
>         fi;
>         ans:=array([seq(0,i=1..m)]);
>         if type(a,`*`) then
>             for i in a do
>                 if type(i,`^`) then ans[op(op(1,i))]:=op(2,i) 
>                 elif type(i,indexed) then ans[op(i)]:=1
>                 fi;
>             od;
>         elif type(a,`^`) then ans[op(op(1,a))]:=op(2,a)
>         elif type(a,indexed) then ans[op(a)]:=1
>         fi;
>         RETURN(convert(ans,list));
>       elif type(a,{TensorMilnor,TensorDual}) then
>         RETURN(map(`convert/List`,a));
>       fi;
>       args;
>    end:
> 
>    `convert/Dual`:=proc(a::{List,Dual})
>       # don't convert Duals (this takes care of 0,1)
>       if type(a,Dual) then RETURN(a)
>       # handle nonsums
>       elif type(a,list(nonnegint)) then
> RETURN(mul(xi[i]^a[i],i=1..nops(a))) 
>       # if the List is a tensor of Lists
>       elif type(a,specfunc(List,T)) then RETURN(map(`convert/Dual`,a))
>       # otherwise it's a list of Lists
>       else RETURN(convert(map(`convert/Dual`,a),`+`) mod 2)
>       fi;
>    end:
> 
>    `convert/P_s`:=proc(a::{List,P_s})
>       # don't convert P_s
>       if type(a,P_s) then RETURN(a)
>       # handle nonsums
>       elif type(a,list(nonnegint)) then
> RETURN(mul(X[i]^a[i],i=1..nops(a)))  
>       # if the List is a tensor
>       elif type(a,specfunc(List,T)) then RETURN(map(`convert/P_s`,a))
>       # otherwise it's a list(List)
>       else RETURN(convert(map(`convert/P_s`,a),`+`) mod 2)
>       fi;
>    end:
> 
>    `convert/Milnor`:=proc(a::{Milnor,List,Inert})
>       # don't convert Milnors (handles 0,1 too)
>       if type(a,Milnor) then RETURN(a)
>       # do sums mod 2 (because of Inerts)
>       elif type(a,`+`) then RETURN(map(`convert/Milnor`,a) mod 2)
>       # do ordinary lists of nonnegint
>       elif type(a,list(nonnegint)) then RETURN(Sq(op(a)))
>       # Tensored Lists to TensorMilnor
>       elif type(a,specfunc(List,T)) then
> RETURN(map(`convert/Milnor`,a))
>       # list of Lists to Milnor
>       elif type(a,list(List)) then
> RETURN(convert(map(`convert/Milnor`,a),`+`) mod 2)
>       # Inert to Milnor
>       elif type(a,Inert) then
>         RETURN(eval(subs( `&x`=`&*` , a)) mod 2)
>       fi;
>    end:
> 
>    `convert/Inert`:=proc(a::{List,Inert})
>      # don't convert Inert
>      if type(a,Inert) then RETURN(a)
>      # do nonsums mod 2
>      elif type(a,list(nonnegint)) then RETURN(&x(op(map(Sq,a))))
>      # tensored lists
>      elif type(a,specfunc(List,T)) then RETURN(map(`convert/Inert`,a))
>      # list(Lists)
>      elif type(a,list(List)) then
> RETURN(convert(map(`convert/Inert`,a),`+`) mod 2)
>      fi;
>    end:
> 
>   
> `convert/Admissible`:=proc(a::{List,Milnor,MonomialSqn,`'+'`(MonomialS
> qn)})
>      local t;
>      # don't convert Admissible's
>      if type(a,Admissible) then RETURN(a)
>      # do sums mod 2, of course
>      elif type(a,`+`) then RETURN(map(`convert/Admissible`,a) mod 2)
>      # use Adem relations for MonomialSqn
>      elif type(a,MonomialSqn) then RETURN(Adem(a))
>      # Lists and Inerts to Milnor first
>      elif not type(a,Milnor) and type(a,{List,Inert}) then 
>         RETURN(convert(convert(convert(a,Inert),Milnor),Admissible));
>      # Milnors
>      elif type(a,Milnor) then
>         t:=Gamma(a);
>         RETURN(t+convert(a+convert(t,Milnor) mod 2,Admissible) mod 2)
>      fi;
>    end:
> 
>    `print/&x`:=proc() 
>       local i,ans,s;
>       ans:=``;
>       for i to nargs do
>         s:=convert(args[i],string);
>         ans:=``||ans||s;
>       od;
>       ans;
>    end:
> 
>    `print/T`:=proc() 
>       local i,ans,s;
>       if type([args],list(Dual)) then RETURN('T'(args))fi;
>       ans:=convert(args[1],string);
>       for i from 2 to nargs do
>         s:=convert(args[i],string);
>         ans:=``||ans||`@`||s;
>       od;
>       ans;
>    end:
> 
>   Sq:=proc()
>      local i;
>      if not type([args],list(integer)) then 
>         ERROR(`invalid arguments`)
>      fi;
>      # convenient to say it is 0 if any index is negative
>      if not type([args],list(nonnegint)) then RETURN(0) fi;
>      for i from nargs by -1 to 1 while args[i]=0 do od;
>      if i=0 then RETURN(1)
>      else RETURN('procname'(args[1..i]))
>      fi;
>   end:
> 
>   xi[0]:=1:
> 
>   T := proc()
>     local i,j,ans,newargs,f;
> 
>     newargs:=[args];
>     # T() = 0
>     if nargs<1 then RETURN(0) fi;
>     # check argument types
>     if not type(newargs,list({TensorMilnor,TensorDual,List})) then 
>         ERROR(`invalid arguments`)
>     fi;
>     # T(a) = a
>     if nargs=1 then RETURN(args) fi;
>     # zero=0
>     if member(0,{args}) then RETURN(0) fi;
>     # associative
>     newargs:=map(proc(x)
>                      if type(x,function) and op(0,x)=`T` then
>                          op(x)
>                      else
>                          x
>                      fi;
>                  end,
>                  newargs);
>     # multilinearity
>     for i to nops(newargs) while not type(newargs[i],`+`) do od;
>     if i<=nops(newargs) then
>        
> f:=subs(a1=op(newargs[1..i-1]),a2=op(newargs[i+1..nops(newargs)]),
>                 proc(x)
>                     `T`(a1,x,a2)
>                 end); 
>         RETURN(map(f,newargs[i]) mod 2);
>     fi;
>     'procname'(op(newargs));    
>   end:
> 
>   Q := proc(k,s)
>     if k=s then RETURN(1) fi:
>     if s<0 then RETURN(0) fi;
>     Q(k-1,s-1)^2+V(k)*Q(k-1,s);
>   end:
> 
> end: # of init
> 
> IsInAn:=proc(a::Milnor,n::nonnegint)
>     local i;
>     if type(a,`+`) then
>         for i in [op(a)] do
>             if not IsInAn(i,n) then RETURN(false) fi;
>         od;
>     else
>         for i from 1 to nops(a) do
>             if not op(i,a)<2^(n+2-i) then RETURN(false) fi;
>         od;
>     fi;
>     true;
> end:
> 
> `&rlex`:=proc(R::{MilnorBasis,List,MonomialSqn},
>               S::{MilnorBasis,List,MonomialSqn})
>     local m,n,i,r,s;
>     r:=R; s:=S;
>     if type(r,{MilnorBasis,MonomialSqn}) then r:=convert(r,List) fi;
>     if type(s,{MilnorBasis,MonomialSqn}) then s:=convert(s,List) fi;
>     if type(r,0) and not type(s,0) then RETURN(true) fi;
>     if type(r,1) and not type(s,{0,1}) then RETURN(true) fi;
>     m:=nops(r); n:=nops(s);
>     r:=[op(r)]; s:=[op(s)];
>     if m<n then RETURN(true) fi;
>     if m=n then
>         for i from m by -1 to 1 while r[i]=s[i] do od;
>         if i>0 and r[i]<s[i] then RETURN(true) fi;
>     fi;
>     false;
> end:
> 
> Gamma:=proc(a::Milnor)
>     local ans,i,dig; 
>     if a=0 or a=1 then RETURN(a)
>     elif type(a,`+`) then ERROR(`can't handle sums yet`)
>     else
>         ans:=NULL; dig:=0;
>         for i from nops(a) by -1 to 1 do
>             dig:=op(i,a)+2*dig;
>             ans:=dig,ans;
>         od;
>         RETURN(convert([ans],Inert));
>     fi;
> end:
> 
> Adem:=proc(a::{MonomialSqn,`+`(MonomialSqn)})
>     local i,A,B,C,S;
>     options remember;
>     # don't change if already Admissible
>     if type(a,Admissible) then RETURN(a) 
>     # do sums term by term
>     elif type(a,`+`) then RETURN(map(procname,a) mod 2)
>     # Monomials in Sq(n)
>     else
>         # find the first non admissible pair
>         for i to nops(a)-1 while op(1,op(i,a))>=2*op(1,op(i+1,a)) do
> od;
>         A:=op(1,op(i,a));  B:=op(1,op(i+1,a));
>         S:=add(multinom(B+C-A-1,A-2*C)*(Sq(A+B-C) &x Sq(C)),
>                C=max(A-B+1,0)..floor(A/2)) mod 2;
>         RETURN(Adem(eval(subsop(i=(S),i+1=NULL,a)) mod 2));
>     fi;
> end:
> 
> Flip:=proc(a::{TensorMilnor,TensorDual})
>     if type(a,{identical(0),identical(1)}) then RETURN(a)
>     elif type(a,TensorMilnor) then 
>         RETURN(convert(convert(a,List),Dual))
>     elif type(a,TensorDual) then 
>         RETURN(convert(convert(a,List),Milnor)) 
>     fi;
> end:
> 
> BinaryChart:=proc(R::MilnorBasis)
>     local i,j,m,r:
>     r:=[op(convert(R,List))];
>     m:=nops(r);
>     for j from max(op(map(omega,r)))-1 by -1 to 0 do
>         for i to m do
>             if j<omega(r[i]) or j=0 then 
>                 printf(`%1d `,Alpha(j,r[i]))
>             else
>                 printf(`  `)
>             fi
>         od;
>         printf(`\n`);
>     od;
> end:
> 
> `&*`:=proc()
>     if type([args],list(TensorMilnor)) then
> RETURN(TensorMilnorMult(args))
>     elif type([args],list(TensorDual)) then
> RETURN(TensorDualMult(args))
>     else ERROR(`All arguments must have type TensorMilnor or
> TensorDual`)
>     fi;
> end:
> 
> TensorMilnorMult:= proc(R::TensorMilnor,S::TensorMilnor)
>     local ans,n,i;
>     if type([args],list(Milnor)) then RETURN(MilnorMult(args)) fi; 
>     if nargs>2 then 
>         ans:=args[1] &* args[2]; 
>         RETURN(procname(ans, args[3..nargs])) 
>     fi;
>     # from here down there are only two args
>     if type(S,`+`) then 
>         RETURN(map(unapply('TensorMilnorMult'(R,x),x),S) mod 2) 
>     fi;
>     if type(R,`+`) then RETURN(map(TensorMilnorMult,R,S) mod 2) fi;
>     #  check tensor lengths
>     if not type(R,specfunc(Milnor,T)) or 
>        not type(S,specfunc(Milnor,T)) then
>         ERROR(`Invalid arguments`);
>     fi;
>     n:=nops(R);
>     for i from 1 to nargs do 
>         if nops(args[i])<>n then 
>             ERROR(`All Tensors must have the same dimension`) 
>         fi;
>     od;
>     # Multiply them
>     T(op(zip((x,y)->x &* y,[op(R)],[op(S)])));
> end:
> 
> TensorDualMult:= proc(R::TensorDual,S::TensorDual)
>     local ans,n,i;
> 
>     if type([args],list(Dual)) then RETURN(DualMult(args)) fi; 
>     if nargs>2 then 
>         ans:=args[1] &* args[2];
>         RETURN(procname(ans, args[3..nargs]))
>     fi;
>     # from here down there are only two args
>     if type(S,`+`) then 
>         RETURN(map(unapply('TensorDualMult'(R,x),x),S) mod 2) 
>     fi;
>     if type(R,`+`) then RETURN(map(TensorDualMult,R,S) mod 2) fi;
>     #  check tensor lengths
>     if not type(R,specfunc(Dual,T)) or not type(S,specfunc(Dual,T))
> then
>         ERROR(`Invalid arguments`);
>     fi;
>     n:=nops(R);
>     for i from 1 to nargs do 
>         if nops(args[i])<>n then 
>             ERROR(`All Tensors must have the same dimension`) 
>         fi;
>     od;
>     # Multiply them
>     T(op(zip((x,y)->x * y,[op(R)],[op(S)])));
> end:
> 
> DualMult:= proc()
>     expand(convert([args],`*`) mod 2)
> end:
> 
> MilnorMult:= proc(R::Milnor,S::Milnor)
>     local lr,ls,lm,M,x,y,i,j,k,n,T,FOUND,ans,SUM;
>     if nargs>2 then 
>         ans:=args[1] &* args[2]; 
>         RETURN(procname(ans, args[3..nargs])) 
>     fi;
>     if type(S,`+`) then
>         RETURN(map(unapply('MilnorMult'(R,x),x),S) mod 2)
>     fi;
>     if type(R,`+`) then RETURN(map(MilnorMult,R,S) mod 2) fi;
>     if R=0 or S=0 then RETURN(0) fi;
>     if R=1 then RETURN(S) elif S=1 then RETURN(R) fi;
>     lr:=nops(R); ls:=nops(S); lm:=lr+ls; ans:=0;
>     #       Initialize the matrix            
>     M:=array(0..lr,0..ls,[seq([seq(0,i=0..ls)],j=0..lr)]);
>     for i to lr do M[i,0]:=op(i,R) od;
>     for i to ls do M[0,i]:=op(i,S) od;          
>     #       Main Computation Loop    
>     FOUND:=true;
>     while FOUND do
>         #    Check the matrix diagonals
>         for n to lm while 
>             multinom(seq(M[n-i,i],i=max(n-lr,0)..min(ls,n)))=1 do
>         od;
>         #   If matrix was good...compute T   
>         if n>lm then
>           T:=array(1..lm,[seq(0,i=1..lm)]);
>             for n to lm do
>                 T[n]:=add(M[n-j,j],j=max(n-lr,0)..min(ls,n))
>             od;
>             # Truncate T and save it
>             for i from lm by -1 to 0 while T[i]=0 do od;
>             ans:=ans+Sq(seq(T[j],j=1..i));
>         fi;
>         #   Find the position to increment   
>         FOUND:=false;
>         for x to lr while not FOUND do
>             SUM:=M[x,0];
>             for y from 1 to ls while not FOUND do
>                 # is there enough to the left of it to increment it 
>                 if SUM>=2^y then
>                     # if all zeros are above it then it's no good 
>                     for k from 0 to x-1 while M[k,y]=0 do od;
>                     if k<x then
>                         # found a candidate 
>                         FOUND:=true;
>                         #    ...and make the new matrix
>                         for i to x-1 do
>                             M[i,0]:=op(i,R);
>                             for j to ls do 
>                                 M[0,j]:=M[0,j]+M[i,j];
>                                 M[i,j]:=0; 
>                             od;
>                         od;
>                         for j to y-1 do
>                             M[0,j]:=M[0,j]+M[x,j];
>                             M[x,j]:=0;
>                         od;
>                         M[0,y]:=M[0,y]-1;
>                         M[x,y]:=M[x,y]+1;
>                         M[x,0]:=SUM-2^y;
>                     else
>                         SUM:=SUM+M[x,y]*2^y;
>                     fi;
>                 else
>                     SUM:=SUM+M[x,y]*2^y;
>                 fi;
>             od;
>         od;
>     od;
>     ans mod 2;
> end:
> 
> `&^` := proc(a::{TensorMilnor,TensorDual},p::nonnegint)
>     options remember;                  
>     if p=0 then RETURN(1) 
>     elif p=1 then RETURN(a)
>     else RETURN(`&*`(a,`&^`(a,p-1)));
>     fi;
> end:
> 
> `&x`:=proc()
>     local i,j,ans,newargs,f;
>     # this routine is almost exactly like tensor product, T...
>     newargs:=[args];
>     if nargs<1 then 
>         ERROR(`&x product must have at least one argument.`); 
>     fi;
>     if nargs=1 then RETURN(args) fi;
>     # zero=0
>     if member(0,newargs) then RETURN(0) fi;
>     # identity=1
>     for i to nops(newargs) do
>         if type(newargs[i],1) then
>            
> RETURN(&x(op(newargs[1..i-1]),op(newargs[i+1..nops(newargs)])))
>         fi;
>     od;
>     # associative
>     newargs:=map(proc(x)
>                      if type(x,function) and op(0,x)=`&x` then
>                          op(x)
>                      else
>                          x
>                      fi;
>                  end,
>                  newargs);
>     # multilinearity
>     for i to nops(newargs) while not type(newargs[i],`+`) do od;
>     if i<=nops(newargs) then
>        
> f:=subs(a1=op(newargs[1..i-1]),a2=op(newargs[i+1..nops(newargs)]),
>                 proc(x)
>                     `&x`(a1,x,a2)
>                 end); 
>         RETURN(map(f,newargs[i]) mod 2);
>     fi;
>     'procname'(op(newargs));
> end:
> 
> `&xx` := proc(a::Inert,p::nonnegint)                  
>     if p=0 then RETURN(1)
>     elif p=1 then RETURN(a)
>     else RETURN(`&x`(a,`&xx`(a,p-1)) mod 2)
>     fi;
> end:
> 
> Coproduct:=proc(a::{Milnor,Dual})
>     local i;
>     # 0 and 1
>     if type(a,0) then RETURN(0) fi;
>     if type(a,1) then RETURN(T(1,1)) fi;
>     # Milnor first
>     if type(a,Milnor) then
>         if type(a,`+`) then RETURN(map(Coproduct,a) mod 2)
>         else 
>             RETURN(convert([seq(T(Sq(op(i[1])),Sq(op(i[2])) ),
>                             i=Parts([op(a)]))],`+`) );
>         fi;
>     # now Dual
>     else
>         if type(a,`+`) then RETURN( expand(map(Coproduct,a)) mod 2)
>         elif type(a,`*`) then 
>             RETURN( expand(`&*`( op( map(Coproduct,[op(a)]) ) ) ) mod
> 2 )
>         elif type(a,`^`) then 
>             RETURN( expand(Coproduct(op(1,a)) &^ op(2,a)) mod 2 )
>         elif type(a,xi[posint]) then 
>             RETURN(add(T(xi[op(a)-i]^(2^i),xi[i]),i=0..op(a)))
>         fi;
>     fi;
>     'procname(args)';
> end:
> 
> chi := proc(a::{Milnor,Dual})
>     local f,cp;
>     options remember;
>     # handle 0,1 first
>     if type(a,{0,1}) then RETURN(a) 
>     # distribute sums and Dual products
>     elif type(a,{`+`,`*`}) then RETURN(expand(map(chi,a)) mod 2)
>     # Dual powers
>     elif type(a,`^`) then RETURN(expand(chi(op(1,a))^op(2,a)) mod 2)
>     # xi[n]
>     elif type(a,xi[posint]) then
>         RETURN(expand( add( xi[op(a)-i]^(2^i)*chi(xi[i]),
>                             i=0..op(a)-1) ) mod 2)
>     # Milnors
>     elif type(a,Milnor) then
>         # make a function to convert T(a,b) to a &* chi(b)
>         f:=x-> op(1,x) &* chi(op(2,x));
>         # the appropriate terms from the coproduct
>         cp:=Coproduct(a)-T(1,a);
>         if type(cp,`+`) then RETURN(map(f,cp) mod 2)
>         else RETURN(f(cp) mod 2)
>         fi;
>     fi;
>     'procname(args)';
> end:
> 
> ActOn:=proc(a::Milnor,p::P_s)
>     local m,P,i,cp,s,k;
> 
>     P:=expand(p) mod 2;
>     if type(P,`+`) then
>         RETURN(expand( map( unapply('ActOn'(a,x),x),P) ) mod 2)
>     fi;
>     if type(a,1) then RETURN(P) fi;
>     if type(a,`+`) then RETURN(expand(map(ActOn,a,P)) mod 2) fi;
>     # P is not a sum, so it is a monomial
>     if type(P,X[integer]) then
>         s:=op(1,P);
>         RETURN(multinom(1-Excess(a),op(a))*X[s]^(1+DEGREE(a)))
>     fi;
>     if type(P,`^`) then
>         k:=op(2,P); s:=op(1,op(1,P));
>         RETURN(multinom(k-Excess(a),op(a))*X[s]^(k+DEGREE(a)))
>     fi;
>     if type(P,`*`) then
>         m:=nops(P);
>         cp:=Coproduct(a);
>         s:=seq(ActOn(op(1,i),op(1,P))*ActOn(op(2,i),P/op(1,P)),i=cp);
>         RETURN(expand(convert([s],`+`)) mod 2);
>     fi;
> end:
> 
> K := proc(R::Milnor,a::Inert)
>     local cp,s,i;
>     if type(R,`+`) then 
>         ERROR(`Kristensen operators are only defined for Basis
> elements.`)
>     fi;
>     if type(a,`+`) then RETURN(map(unapply('K'(R,x),x),a) mod 2) fi;
>     if type(R,0) then RETURN(0) fi;
>     if type(R,1) then RETURN(a) fi;
>     if type(a,specfunc(nonnegint,Sq)) then
>         if nops(R)>nops(a) then RETURN(0)
>         else RETURN(Sq(op(zip((x,y)->x-y,[op(a)],[op(R)],0)) ) );
>         fi;
>     fi;
>     if type(a,specfunc(Milnor,`&x`)) then
>         cp:=map(Flip,[op(Coproduct(Flip(R)))]);
>         s:=seq(K(op(1,i),op(1,a)) &x
> K(op(2,i),&x(op(2..nops(a),a))),i=cp);
>         RETURN(convert([s],`+`) mod 2);
>     fi;
>     'procname(args)'
> end:
> 
> Excess:=proc(a::{Milnor,Inert})
>     local f,s,i,t,ans;
>   
>     if type(a,0) then ERROR(`The excess of 0 is undefined`) fi;
>     # not the most efficient way to do inerts, but hey...
>     if not type(a,Milnor) then RETURN(Excess(convert(a,Milnor))) fi;
>     # Now Milnors
>     if type(a,`+`) then
>         if member(1,[op(a)]) then RETURN(0) fi;
>         min(op(map(unapply('convert([op(i)],`+`)',i),[op(a)])))
>     else
>         convert([op(a)],`+`)
>     fi
> end:
> 
> DEGREE := proc(a::{TensorMilnor,TensorDual})
>     local i,ans;
>     # sums
>     if type(a,`+`) then
>         ans:=map(DEGREE,{op(a)});
>         if nops(ans)=1 then RETURN(op(ans)) else RETURN(ans) fi;
>     fi;
>     # 0,1 have degree 0
>     if type(a,{0,1}) then RETURN(0) fi;
>     # a pure Milnor
>     if type(a,Milnor) then
>         RETURN(add((2^i-1)*op(i,a),i=1..nops(a)))
>     fi;
>     # a pure Dual - trick: convert to a Milnor via list first
>     if type(a,Dual) then 
>         RETURN(DEGREE(convert(convert(a,List),Milnor)))
>     fi;
>     # a Tensored Milnor or Dual that's not Milnor or Dual
>     if type(a,{TensorMilnor,TensorDual}) then
>         RETURN(add(DEGREE(op(i,a)),i=1..nops(a)))
>     fi;
>     'procname(args)'
> end:
> 
> May:=proc(a::Milnor)
>     local x,i,ans;
> 
>     if type(a,{0,1}) then RETURN(0) fi;
>     if type(a,`+`) then
>         RETURN(min(op(map(May,{op(a)}))));
>     else
>         RETURN(add(alpha(op(i,a))*i,i=1..nops(a)))
>     fi;        
> end:
> 
> Filtration:=proc(R::Milnor)
>     local g,x,y;
>     g:=convert(Gamma(R),Milnor);
>     if g=R then 
>           RETURN(0) 
>     else
>           x:=g+R mod 2;
>           if type(x,`+`) then y:=[op(x)]
>           else y:=[x];
>           fi;
>           RETURN(max(op(map(Filtration,y) ) )+1);
>     fi;
> end:
> 
> AnBasis:=proc(g::nonnegint)
>     local c,m,i,n;
>     if nargs=1 then n:=infinity
>     else n:=args[2]
>     fi;
>     if not type(n,{nonnegint,infinity}) then 
>         ERROR(`Invalid value of n`); 
>     fi;
>     if type(n,infinity) then 
>         if g<2 then m:=1
>         else m:=ceil(evalf(ln(g)/ln(2)))
>         fi;
>         c:=[seq(2^i-1,i=1..m)];
>         RETURN(op(sort(map(x->Sq(op(x)),IntSolve(g,c)),`&rlex`)))
>     fi;
>     m:=n+1;
>     c:=[seq(2^i-1,i=1..m)];
>    
> op(sort(select(IsInAn,map(x->Sq(op(x)),IntSolve(g,c)),n),`&rlex`));
> end:
> 
> AdmisAnBasis:=proc(g::nonnegint)
>     local n;
>     if nargs=1 then n:=infinity
>     else n:=args[2]
>     fi;
>     op(map(Gamma,[AnBasis(g,n)]));
> end:
> 
> macro(CBM=`ChangeBasisMatrix`);
> 
> ChangeBasisMatrix:=proc(g::nonnegint)
>     local BA,BM,M,n,i,j,x,y,k;
> 
>     BA:=[AdmisAnBasis(g)];
>     BM:=[AnBasis(g)];
>     n:=nops(BA);
>     M:=matrix(n,n,0);
>     for i from 1 to n do
>         x:=BA[i];
>         y:=convert(x,Milnor);
>         for j from 1 to n do
>             M[j,i]:=coeff(y,BM[j]);
>         od;
>     od;    
>     eval(M);
> end:
> 
> Nil:=proc(a::Milnor,t::posint)
>     local n,s,st;
>     s:=a; st:=time();
>     for n from 1 while s<>0 and time()-st<t*60 do s:=s &* a; od;
>     if s<>0 then 
>         printf(`%3d -- Gave up after %d minutes %d seconds\n`,
>                op(a),
>                iquo(trunc(time()-st),60),
>                irem(trunc(time()-st),60));
>     else printf(`%3d --%3d\n`,op(a),n)
>     fi;
> end:
> 
> Bottom:=proc(A::array) op(1,op(2,eval(A))) end:
> 
> Top:=proc(A::array) op(2,op(2,eval(A))) end:
> 
> IntSolve:=proc(g::nonnegint,c::list(posint))
>     local t,i,j,k,ans;
> 
>     t:=nops(c);
>     if t=1 then
>         if g mod c[1] = 0 then RETURN([g/c[1]])
>         else RETURN([]);
>         fi;
>     else
>         ans:=NULL;
>         for i from 0 to iquo(g,c[1]) do
>             for j in IntSolve(g-i*c[1],c[2..t]) do
>                 ans:=[i,op(j)],ans;
>             od;
>         od;
>     fi;
>     [ans];
> end:
> 
> Parts:=proc(R)
>     local m,i,j,ans;
>     m:=nops(R);
>     if m=1 then RETURN([seq([[i],[R[1]-i]],i=0..R[1])])
>     else
>         ans:=NULL;
>         for j in Parts(R[2..m]) do
>             ans:=ans,seq([[i,op(j[1])],[R[1]-i,op(j[2])]],i=0..R[1]);
>         od;
>     fi:
>     RETURN([ans]):
> end:
> 
> multinom:=proc()
>     local f,L,S;
>     L:=[args];
>     if not type(L,list(integer)) then ERROR(`Invalid arguments`); fi;
>     # if any arguments are negative, then return zero
>     if not type(L, list(nonnegint)) then RETURN(0) fi;
>     f:=n->(n-(n mod 2))/2;
>     S:=convert(L,set);
>     if nops(S)=1 and op(1,S)=0 then RETURN(1)
>     elif convert(L mod 2,`+`) > 1 then RETURN(0)
>     else multinom(op(map(f,L)))
>     fi;
> end:
> 
> Alpha:=proc(i::integer,n::integer);
>     RETURN((n-(n mod 2^i))/2^i mod 2)
> end:
> 
> alpha:=proc(n::integer)
>     local n1;
>     options remember;
>     if n=0 then RETURN(0):
>     else
>         n1:=n mod 2:
>         n1+alpha((n-n1)/2):
>     fi:
> end:
> 
> nu:=proc(n::integer)
>     options remember;
> 
>     if n=0 then RETURN(0):
>     elif n mod 2 = 1  then RETURN(0):
>     else
>         nu(n/2)+1:
>     fi:
> end:
> 
> omega:=proc(n::integer)
>     if n=0 then RETURN(0)
>     elif type(n,odd) then RETURN(omega((n-1)/2)+1);
>     else RETURN(omega(n/2)+1);
>     fi;
> end:
> 
> s1:=proc(u)
>     options remember;
>     if type(u,even) then RETURN(0)
>     else RETURN(s1((u-1)/2)+1)
>     fi;      
> end:
> 
> s2:=u->s1(u+2^s1(u)):
> 
> w:=proc(R::list(integer))
>     local i;
>     add(2^(i-1)*R[i],i=1..nops(R));
> end:
> 
> delta:=proc(x,y) if x=y then 1 else 0 fi end:
> 
> id:=proc(n) if n=0 then RETURN(NULL) else matrix(n,n,delta) fi end:
> 
> ZERO:=proc(n,m) if n=0 or m=0 then RETURN(NULL) else matrix(n,m,0) fi
> end:
> 
> Equal:=proc(A,B)
>     local a1,a2,b1,b2,i;
> 
>     a1:=Bottom(A); a2:=Top(A);
>     b1:=Bottom(B); b2:=Top(B);
>     if a1<>b1 or a2<>b2 then RETURN(false) fi;
>     for i from a1 to a2 do
>         if A[i]<>B[i] then RETURN(false) fi;
>     od;
>     true;
> end:
> 
> Inverse:=proc(M)
>     local n,r,m;
> 
>     n:=linalg[rowdim](M);
>     Gaussjord(M,r) mod 2; m:=linalg[coldim](M);
>     if r<>n or m<>n then ERROR(`Matrix not invertible:`,eval(M)) fi;
>     linalg[submatrix](Gaussjord(linalg[augment](M,id(n))) mod
> 2,1..n,n+1..2*n);
> end:
> 
> FindPivotColumns:=proc(B)
>     local i,j,m,n,ans;
> 
>     m:=linalg[rowdim](B);
>     n:=linalg[coldim](B);
>     ans:=[];
>     for i from 1 to m do 
>         for j from 1 to n while B[i,j]=0 do od;
>         if j<=n then ans:=[op(ans),j] fi;
>     od;
>     ans;
> end:
> 
> ExtendToBasis:=proc(B)
>     local m,A,L;
> 
>     m:=linalg[rowdim](B);
>     A:=linalg[augment](B,id(m));
>     L:=FindPivotColumns(Gaussjord(A) mod 2);
>     linalg[submatrix](A,1..m,L);
> end:
> 
> ReduceToBasis:=proc(B)
>     local m,L;
>     m:=linalg[rowdim](B);
>     L:=FindPivotColumns(Gaussjord(B) mod 2);
>     if nops(L)=0 then RETURN(NULL) fi;
>     linalg[submatrix](B,1..m,L);
> end:
> 
> NULLSpace:=proc(B)
>     local A,L,CL,m,n,r,C1,C2,i,i1,i2,V;
> 
>     m:=linalg[rowdim](B); n:=linalg[coldim](B);
>     A:=Gaussjord(B,r) mod 2; # A is m x n
>     if n=r then RETURN(NULL) fi;
>     L:=FindPivotColumns(A);
>     CL:=sort(convert({seq(i,i=1..n)} minus convert(L,set),list));
>     C1:=linalg[submatrix](A,1..m,CL); # C1 is m x n-r
>     C2:=id(n-r); # C2 is n-r x n-r of course
>     i1:=1; i2:=1; V:=array(1..n);
>     for i to n do
>         if member(i,L) then
>             V[i]:=linalg[subvector](C1,i1,1..n-r); i1:=i1+1;
>         else
>             V[i]:=linalg[subvector](C2,i2,1..n-r); i2:=i2+1;
>         fi;
>     od;
>     linalg[stackmatrix](seq(eval(V[i]),i=1..n));
> end:
> 
> v:=proc(L)
>     local i,ans;
>     add(L[i]*X[i],i=1..nops(L))+X[nops(L)+1];
> end:
> 
> V:=proc(k)
>     local ans,i,Z2,x;
> 
>     if k=1 then RETURN(X[1]) fi:
>     ans:=1;
>     Z2:=[0,1];
>     x:=combinat[cartprod]([seq(Z2,i=1..k-1)]);
>     while not x[finished] do 
>       i:=x[nextvalue]();
>       ans:=ans*v(i);
>     od;      
>     ans;
> end:
> 
> tri:=proc(u,v) if v<2^s1(u) then true else false fi end:
> 
> SqDickson:=proc(i,x,k)
>     local s,t,a,kset,r;
>     options remember;
>     if  i=0 then x
>     elif type(x,`+`) then map(unapply(Sq(i,s),s),x) mod 2
>     elif type(x,`*`) then 
>         a:=0;
>         for t from 0 to i do
>             a:=a+Sq(i-t,op(1,x))*Sq(t,convert([op(2..nops(x),x)],`*`)
> );
>         od;
>         expand(a) mod 2;
>     elif type(x,`^`) then
>         a:=0;
>         for t from 0 to i do
>             a:=a+Sq(i-t,op(1,x))*Sq(t,op(1,x)^(op(2,x)-1));
>         od;
>         expand(a) mod 2;
>     elif type(x,indexed) then
>         s:=op(1,x);
>         if i=2^k-2^s then RETURN(op(0,x)[s]^2) fi;
>         for r from 0 to s do
>             for t from s+1 to k do
>                 if i=2^k-2^t+2^s-2^r then
>                     RETURN(op(0,x)[t]*op(0,x)[r])
>                 fi;
>             od;
>         od;
>         0
>     else 'Sq'(i,x)
>     fi;
> end:
> 
> MBasis:= proc(r::integer, alims, exchk, ipow::integer,
> initize::integer) 
>     local t, i, j, k, toppow, lims, ans, pow, rr, qr, chk; 
> 
>     description " MBasis(deg,n,ex) constructs the Milnor basis for the
> Steenrod Algebra "
>     "in degree deg. If n is supplied, it constructs a basis for An. " 
>     "If ex is supplied, it gives a basis for the elements of excess at
> most ex."
>     " To specify ex without n, use MBasis(deg,infinity,ex)" ;
> 
>     if nargs < 5  then  
>       t := ilog[2](r + 1); 
>       pow := 2^t - 1;
>       if nargs < 3 then chk := infinity else chk := exchk; fi;
>       if nargs < 2 then 
>         lims := infinity 
>       else	   
>         lims := 2^(alims + 2 - t);
>       fi; 
>     else	
>       pow := ipow; 
>       t := initize; 
>       lims := alims;
>       chk:= exchk;     
>     fi;	
>     toppow := floor(lims) - 1;
>     rr := irem(r, pow,'qr');   
>     if (pow = 1) then    
>       `if`(((rr = 0) and ((qr <= toppow) and (qr <= chk)) ), 
>           RETURN([[qr]]), RETURN([])); 
>     fi;
>     ans := [];
>     for i from min(qr, toppow) by -1 to 0 do 
>       if (chk < i) then next; fi;
>       ans := [op(ans), op(map((a,b)->[op(a),b], 
>               MBasis(r - i*pow, lims*2, chk - i, (pow - 1)/2, t -
> 1),i))];
>     od;
>     RETURN(ans);
> end:
> 
> end module:
# Save to Personal Lib:
> savelib('Steenrod');
> 
>