/* MAGMA code newDL
   R.R.Bruner 2011
   */



/* Take an F_2 module H^*(X) in Ext_1.x format and compute the F_2 module
   which is the cohomology of the N skeleton the extended power
   D_r(X) using the fact that it is the weight r part of QX.  */

/* set N and r, and let file be the name of the file containing the definition
   of H^*(X), as in:

   N := 7;
   r := 8;
   file := "S0";
   load newDL;

/* Write the results in the same format needed by Ext.1.x in the file named
   D_<r><file>to<N>mod, e.g., D_8S0to7mod with the values of N,r and file above.
*/



BC := func<n,i | GF(2)!Binomial(n,i) >;


/* ------ */

function read_f2_mod(file)
  input := StringToIntegerSequence(Read(file));
  degs := input[2..1+input[1]];
  ind := 2+input[1];
  ops := [[] : x in degs];
  while ind le #input do
    ops[1+input[ind]][input[ind+1]] :=
        [1+input[i] : i in [ind+3..ind+2+input[ind+2]]];
    ind +:= 3 + input[ind+2];
  end while;
  if ind ne 1+#input then
    print "ind not exact: ind ended at ",ind, "#input = ",#input;
  end if;

  /* degree check */
  okay := true;
  for g in [1..#degs] do
    for i in [1..#ops[g]] do
      if IsDefined(ops[g],i) then
/*      print g,i,ops[g]; */
      for v in ops[g][i] do
        if not degs[g]+i eq degs[v] then
	   okay := false;
	   printf "\nDegree error: Sq^%o of gen %o contains %o; degree %o+%o ne %o\n",
	       i,g,v,i,degs[g],degs[v]; 
	end if;
      end for;
      end if;
    end for;
  end for;

  /* Now take hom dual of the Steenrod ops */
  dualops := [[] : x in degs];
  for i in [1..#dualops] do
     for j in [1..#ops[i]] do
        if IsDefined(ops[i],j) then
	    for k in ops[i][j] do
	       if IsDefined(dualops[k],j) then
                  dualops[k][j] cat:= [i];
	       else
	          dualops[k][j] := [i];
	       end if;
	    end for;
	end if;
     end for;
  end for;

return degs,dualops,okay;
end function;



/* ------ */

function GenWt(L)
 return 2^(#L-1);
end function;


/* ------ */

function GenDeg(L,degs)
  if #L gt 1 then
     return &+[L[i] : i in [1..#L-1]] + degs[L[#L]];
  else
     return degs[L[1]];
  end if;
end function;



/* ------ */

function MonWt(L)
  return &+[GenWt(x) : x in L];
end function;



/* ------ */

function MonDeg(L,degs)
  return &+[GenDeg(x,degs) : x in L];
end function;



/* ------ */

function prods(w,n,L,degs)   /* returns list of ordered lists of elements of L
                           of weight wt and degree le n.  Assumes L is sorted so
			   that wt is non-decreasing.    The list degs is the list
			   of degrees of generators needed by GenDeg.
			   */
  if w eq 0 then
         return [[]];
  elif #L eq 0 or
     w lt GenWt(L[1]) or
     n lt Min([GenDeg(x,degs) : x in L]) then
        return [];
  else
        return  [[L[1]] cat x 
	         : x in prods(w-GenWt(L[1]),
		              n-GenDeg(L[1],degs),
		              [y : y in L | GenWt(y) le w - GenWt(L[1]) 
			                 and GenDeg(y,degs) le n-GenDeg(L[1],degs)],
			      degs)]
	        cat prods(w,n,L[2..#L],degs);
  end if;
end function;



/* Gen - an operation Q^i1 Q^i2 ... Q^ik(gen x)
         represented as [i1 i2 ... ik x]
   Mon - a product of such, 
         represented as a list of the Gens in the product,
	 normal form has the gens in increasing order.
	 These form the basis.
   Elt - a sum of Mons,
         represented as a list of the Mons in the sum
*/




/* arithmetic ---------------------------------------------------- */


/* ------ */

function MonNF(x,G)  /* reorder terms in the monomial x so that they
                           are in increasing order, w/ G specifying the order */
  return [G[i] : i in Sort([Index(G,g) : g in x])];
end function;



/* ------ */

function EltSumAux(L1,L2)         /* add the Elts L1 and L2,
                                      assuming L1 is not redundant, so
				      no terms in it will cancel. */
  val := {x : x in L1};
  for x in L2 do
    val := val sdiff {x};
  end for;
  return [x : x in val];
end function;



/* ------ */

function EltSum(L)         /* add the Elts in the list L */
  LL := &cat L;
  val := {};
  for x in LL do
    val := val sdiff {x};
  end for;
  return [x : x in val];
end function;



/* ------ */

function EltProd(L,M,gens)     
                 /* multiply Elts L and M, 
		    using list gens of Gens to order monomials  */
  return EltSum([[MonNF(x cat y,gens) : x in L] : y in M]);
end function;

function EltToList(x,L)
  return [Index(L,i) : i in x];
end function;


/* DL operations    ---------------------------------------------------- */

forward DL;

function DLGen(a,G,gens,degs)
                          /* Q^a of Gen G = [i1 i2 ... g].
			    Needs degrees of Base gens in degs and
			    list of Gens in gens.
			    Enforces excess requirements too, using degs.
			    Returns an Elt             */
  if #G eq 0 then
     print "The trivial generator appeared in DLGen\n";
     return false;
  elif #G eq 1 then
     if a lt degs[G[1]] then
       return [];
     elif a eq degs[G[1]] then
       return [[ [G[1]],[G[1]] ]];
     else
       return [[ [a] cat G]];
     end if;
  else
     /* Really ought to check that G is Admissable somewhere, no? */
     gd := GenDeg(G,degs);
     if a lt gd then
        return [];
     elif a eq gd then
        return [[G,G]];
     elif a le 2*G[1] then
        return [[ [a] cat G]];
     else
         return EltSum([DL(a+G[1]-t, DLGen(t,G[2..#G],gens,degs),gens,degs)
                        : t in [Ceiling(a/2)..a-G[1]-1]
		        | BC(t-G[1]-1,2*t-a) ne 0]);
     end if;
  end if;
end function;



/* ------ */

function DLMon(a,M,gens,degs)  
                        /* Q^a of Mon M, returns an Elt  
			 Needs Gens gens for MonNF to work and 
			 Base_degs degs for GenDeg   */
   if #M eq 0 then
     if a eq 0 then
       return [[]];
     else
       return [];
     end if;
   end if;
   dd := MonDeg(M,degs);
   if a lt dd then
      return [];
   elif a eq dd then
      return [MonNF(M cat M,gens)];
   end if;
   val := [];
   for i in [0..a] do 
       val := EltSumAux(val,
                        EltProd(DLGen(i,M[1],gens,degs),
		                DLMon(a-i,M[2..#M],gens,degs),
			        gens));
   end for;
   return val;
end function;



/* ------ */

function DL(a,L,gens,degs)    
                        /* Q^a of Elt L, 
			 Needs Gens gens for MonNF to work and 
			 Base_degs degs for GenDeg   */
  return EltSum([DLMon(a,x,gens,degs) : x in L]);
end function;

/* Squaring operations ------------------------------------------------- */

forward Sq;
forward SqMon;

function SqGen(a,G,gens,degs,ops)
  if #G eq 0 then
     printf "Empty generator in SqGen";
     return false;
  elif a eq 0 then
     return [[G]];
  elif #G eq 1 then
     if IsDefined(ops[G[1]],a) then
       return [ [[g]] : g in ops[G[1]][a]];
     else
       return [];
     end if;
  else
    return EltSum([DL(G[1]-a+t,SqGen(t,G[2..#G],gens,degs,ops),gens,degs)
                   : t in [0..Truncate(a/2)] 
		   | BC(G[1]-a,a-2*t) ne 0]);
  end if;
end function;

function SqMon(a,M,gens,degs,ops)
  if #M eq 0 then
     printf "Empty generator in SqMon";
     return false;
  elif a eq 0 then
     return [M];
  elif #M eq 1 then
    return SqGen(a,M[1],gens,degs,ops);
  else
    return EltSum([ EltProd(SqGen(i,M[1],gens,degs,ops),
                            SqMon(a-i,M[2..#M],gens,degs,ops),gens)
	            : i in [0..a]]);
  end if;
end function;

function Sq(a,x,gens,degs,ops)
  return EltSum([SqMon(a,m,gens,degs,ops) : m in x]);
end function;



/* Now get specific ---------------------------------------------------- */

/* N := 7;
   r := 8;
   file := "S0";
*/
Base_degs, Base_ops, okay := read_f2_mod(file);
print Base_ops;


Gens := [i : i in [1..#Base_degs] | Base_degs[i] le N];

newGens := &cat[ [[i,j] : i in [Base_degs[j]+1..N-Base_degs[j]]]
                 : j in Gens];
Gens := [[i] : i in Gens];


lastwt := 2;
wt := 4;
while wt le r do
  Gens := Gens cat newGens;
  newGens := &cat[ [[i] cat x 
                  : i in [1+GenDeg(x,Base_degs)..Min(N-GenDeg(x,Base_degs),2*x[1])]]
                  : x in newGens];
  lastwt := wt;
  wt +:= wt;
end while;

Gens cat:= newGens;

printf "%o generators\n",#Gens;

/* Gens is now sorted in order of increasing weight.  If it weren't we'd have
   to sort it for prods to work properly.  */

Mons_list := prods(r,N,Gens,Base_degs);

printf "%o monomials\n",#Mons_list;

/* Now we sort Elts by degree */

Mons := Sort(Mons_list,func<x,y | MonDeg(x,Base_degs)-MonDeg(y,Base_degs)>);

printf "Mons sorted\n";

/*
exclude := [i : i in [1..#Mons] | Mons[i][1] eq [1] and Mons[i][2] eq [1]];
Mons := [Mons[i] : i in [1..#Mons] | not i in exclude];
*/

Min_deg := MonDeg(Mons[1],Base_degs);
Max_deg := MonDeg(Mons[#Mons],Base_degs);

Mon_degs := [MonDeg(x,Base_degs) : x in Mons];
Mon_ops := [[] : x in Mons];

for ind in [1..#Mons] do
printf "ops on monomial %o\n",ind;
  x := Mons[ind];
  for i in [1..Mon_degs[ind]-Min_deg] do
     for y in EltToList(SqMon(i,x,Gens,Base_degs,Base_ops),Mons) do
	if y ne 0 then
	if IsDefined(Mon_ops[y],i) then
            Mon_ops[y][i] cat:= [ind];
	else
            Mon_ops[y][i] := [ind];
	end if;
	end if;
     end for;
  end for;
end for;


/* Now write it out -------- */

SetOutputFile("D_" cat IntegerToString(r) cat file cat "to" cat IntegerToString(N) cat "mod");    
printf "%o\n\n",#Mons;

count := 0;
for x in Mon_degs do
printf "%o",x;
if count eq 9 then
   coint := 0;
   printf "\n";
else
  printf " ";
end if;
end for;

printf "\n\n";

for ind in [1..#Mons] do
  for i in [1..#Mon_ops[ind]] do
    if IsDefined(Mon_ops[ind],i) then
      printf "%o %o %o",ind-1,i,#Mon_ops[ind][i];
      for j in Mon_ops[ind][i] do
        printf " %o",j-1;
      end for;
      printf "\n";
    end if;
  end for;
end for;

  UnsetOutputFile();