//Calculation of the Ekedahl-Oort type of a plane curve
//Magma code written by Wieb Bosma and Ben Moonen
//with changes by an anonymous referee
//Version of February 16, 2022

HWtriple := function(k, d, f)
  p := Characteristic(k);
  P<X0,X1,X2> := Parent(f);
  g := (d-1)*(d-2) div 2;

  Md := MonomialsOfDegree(P, d-3);  // Md[i] =m_i

  MA := MatrixAlgebra(k, g);
  APhi := MA ! 0;
  F := f^(p-2); // need (p-2)nd power also later
  for i := 1 to g do
    for j := 1 to g do
      APhi[i][j] := MonomialCoefficient(Md[i]*f*F,Md[j]^p * X0^(p-1)*X1^(p-1)*X2^(p-1));
    end for;
  end for;

  if IsZero(APhi) then return "superspecial", 0, 0, 0;
  elif IsInvertible(APhi) then return "ordinary", 0, 0, 0;
  end if;
  kappa := Basis(Kernel(Transpose(APhi))); // kappa[i] = kappa_i
  h := #kappa;

  m := -3*d+3;
  n := -2*d+2;

  // clumsy basis Tm:
  BTm := [];
  efm := -m; // exponent to make Laurentpolys poly
  for i in [m..-1] do
    for j in [m..-1] do
      l := m-i-j;
      if l le -1 then
        Append(~BTm, X0^(i+efm)*X1^(j+efm)*X2^(l+efm));
      end if;
    end for;
  end for;

  BTn := [];
  efn := efm;
  for i in [n..-1] do
    for j in [n..-1] do
      l := n-i-j;
      if l le -1 then
        Append(~BTn, X0^(i+efn)*X1^(j+efn)*X2^(l+efn));
      end if;
    end for;
  end for;

  df :=[]; // partial derivatives of f
  for j := 0 to 2 do
    Append(~df, Derivative(f, P.(j+1)));
  end for;

  KM := KMatrixSpace(k,#BTm,3* #BTn);
  VS := VectorSpace(k, #BTm);
  Multdf := KM ! 0;
  for i := 1 to #BTm do
    for j := 1 to #BTn do
      Multdf[i][j] := MonomialCoefficient(df[1]*BTm[i], BTn[j]);
      Multdf[i][#BTn + j] := MonomialCoefficient(df[2]*BTm[i], BTn[j]);
      Multdf[i][2* #BTn + j] := MonomialCoefficient(df[3]*BTm[i], BTn[j]);
    end for;
  end for;
  KerMultdf := Kernel(Multdf);  // should be generated by u:
  BKerMultdf := Basis(KerMultdf); // should have single element!
  if #BKerMultdf ge 2 then print "something's wrong! (computation of the space U)";
  end if;
  utilde := &+[ BKerMultdf[1][i]*BTm[i] : i in [1..#BTm] ]; // shifted version of the elt u


  Mon := MonomialsOfDegree(P, 2*d-3);
  c := (2*d-1)*(d-1); // number of monoms of degree 2d-3

  KMgxc := KMatrixSpace(k,g,c);
  Bmatr := KMgxc ! 0;
  Cmatr := KMgxc ! 0;
  for i := 1 to c do
    for j := 1 to g do
      Bmatr[j][i] := MonomialCoefficient(Mon[i]*Md[j]*utilde,X0^(3*d-4)*X1^(3*d-4)*X2^(3*d-4));
      Cmatr[j][i] := MonomialCoefficient(F*Mon[i],X0^(p-1)*X1^(p-1)*X2^(p-1)*Md[j]^p);
    end for;
  end for;

  KappaMatr := Matrix(k,h,g,kappa);
  check, APsiT, dummy := IsConsistent(Bmatr,KappaMatr * Cmatr);
  if check
    then APsi := Transpose(APsiT);
    else print "something's wrong (solving A(Psi))";
  end if;
  return "interesting", APhi, kappa, APsi;
end function;


function DieudMod(k, d, APhi, kappa, APsi)
  g := (d-1)*(d-2) div 2;
  Vg :=VectorSpace(k, g);
  Bg := Basis(Vg);
  Einit := ExtendBasis(kappa, Vg);
  h := #kappa;
  E := Einit[h+1..#Einit] cat Einit[1..h];
  VV := VectorSpaceWithBasis(E);
  EC := [];
  for j := 1 to g do
    Append(~EC, Coordinates(VV, Vg.j));
  end for;

  MS := KMatrixSpace(k,2*g,g);
  AF := MS ! 0;
  for j := 1 to g do
    ee := EC[j];
    for r := 1 to g do
      AF[r][j] := &+[ ee[mu]*APhi[r][Index(Bg,E[mu])] : mu in [1..g-h]];
    end for;
    for r := g+1 to 2*g do
      AF[r][j] :=&+[ ee[g-h+nu]*APsi[2*g+1-r][nu]: nu in [1..h]];
    end for;
  end for;

  AF := HorizontalJoin(AF, Parent(AF)!0);
  return AF;
end function;



EOSeq:=function(k, d, AF)
  g := (d-1)*(d-2) div 2;
  frobk := FrobeniusMap(k);
  TAF := Transpose(AF);
  MA := MatrixAlgebra(k, g);
  zero := MA ! 0;
  one := MA ! 1;
  SymForm := VerticalJoin(
    HorizontalJoin(zero, ReverseColumns(one)),
    HorizontalJoin(ReverseColumns((-1)*one), zero));
  VSym := VectorSpace(k, 2*g, SymForm);
  BSym := Basis(VSym);
  table_B := [ Parent({@ VSym| @}) | ];
  // initial
  table_B[1] := {@ VSym | @};
  table_B[2*g+1] := {@ BSym[j] : j in [1..2*g]  @};
  table_B[g+1] := {@ VSym ! Transpose(ColumnSubmatrix(AF, i, 1)) : i in [1..g] @};
  table_f := [ Integers() | ];
  // initial:
  table_f[1] := 0;
  table_f[2*g+1] := g;
  numberfs := 2;

  while numberfs lt 2*g+1 do
    found := false;
    for i := 2 to 2*g do
      if IsDefined(table_B, i) and not IsDefined(table_f, i) then
        found := true;
        ifound := i;
        break;
      end if;
    end for;
    if found then
      // Step (2)
      Bj :=[];
      for j := 1 to ifound-1 do
        bj := table_B[ifound][j];
        sigmabj := VSym ! [ frobk(bj[n]) : n in [1..2*g] ];
        Append(~Bj, sigmabj*TAF);
      end for;
      SBj := sub< VSym | Bj >;
      D := Dimension(SBj);
      table_f[ifound] := D;
      numberfs +:= 1;
      if IsDefined(table_B, D+1) then
        continue;
      end if;
      // Step (3)
      table_B[D+1] := {@ b : b in Basis(SBj) @};
      OS := OrthogonalComplement(VSym, SBj);
      table_B[2*g-D+1] := {@ b : b in Basis(OS) @};
    else
      // Step (4)
      while not numberfs eq 2*g+1 do
        for i := 1 to 2*g do
          if not IsDefined(table_f, i) then
            a := i;
            for j := i+1 to 2*g+1 do
              if IsDefined(table_f, j) then
                b := j;
                break i;
              end if;
            end for;
          end if;
        end for;
        //table_f;
        //table_B;
        if IsDefined(table_f, a-1) and table_f[a-1] eq table_f[b] then
          for m := a to b-1 do
            table_f[m] := table_f[a-1];
            numberfs +:= 1;
          end for;
        else
        // table_f[a-1] + (b-a+1) eq table_f[b]  !
          for m := a to b-1 do
            table_f[m] := table_f[a-1] + (m-a+1);
            numberfs +:= 1;
          end for;
        end if;
      end while;
    end if;
  end while;
return table_f;
end function;

WeylGrElt := function(k, d, AF)
  g := (d-1)*(d-2) div 2;
  table_f:=EOSeq(k, d, AF);
  // Step (5)
  Jg := []; Ig := []; w := [];
  for j := 2 to 2*g+1 do
    if table_f[j] eq table_f[j-1] then
      Append(~Jg, j-1);
    else
      Append(~Ig, j-1);
    end if;
  end for;
  for j := 1 to g do
    w[Jg[j]] := j;
    w[Ig[j]] := g+j;
  end for;
  return w, table_f;
end function;


intrinsic EOType(f::RngMPolElt) -> .
  {The Weyl group element and the final type of the p-kernel of the Jacobian of the plane curve defined by f}
  k:= CoefficientRing(Parent(f));
  p:= Characteristic(k);
  P<X0,X1,X2> := Parent(f);
  if not IsHomogeneous(f) then
    return "the defining polynomial must be homogeneous";
    else
      d := TotalDegree(f);
      g := (d-1)*(d-2) div 2;
      if IsDivisibleBy(d,p) then
        return "the degree of the polynomial must not be divisible by p";
      else
        S := ProjectiveSpace(P);
        C := Curve(S,f);
        if IsSingular(C) then
          return "the curve defined by your polynomial is singular";
        else
          Status, APhi, kappa, APsi := HWtriple(k, d, f);
          if Status eq "ordinary" then print "The curve defined by your polynomial is ordinary. \n The Ekedahl-Oort type (the Weyl group permutation, and the final type): ";
		return Sym(2*g)![i:i in [1..2*g]] , [i:i in [0..g]] cat [g : i in [1..g]];
          elif Status eq "superspecial" then print "The curve defined by your polynomial is superspecial. \n  The Ekedahl-Oort type (the Weyl group permutation, and the final type): ";
		return Sym(2*g)!([i: i in [g+1..2*g]] cat [i: i in [1..g]]) , [0:i in [0..g]] cat [i: i in [1..g]];
          else
            AF := DieudMod(k, d, APhi, kappa, APsi);
            w, eo  := WeylGrElt(k, d, AF);
            print "The Ekedahl-Oort type is (the Weyl group permutation, and the final type): ";
            return Sym(2*g)!w, eo;
          end if;
        end if;
      end if;
  end if;
end intrinsic;