\\ Given a subgroup H in HNF and an element h, reduce h to a canonical 
\\ coset representative by subtracting elements of H
coset_reduce(H,h) = forstep(i = 4, 1, -1, h -= (h[i] \ H[i,i]) * H[,i]~); h;

\\ Put a subgroup into a canonical form
subgroup_reduce(H) =
{
  for(i = 1, 3, H[1][i,] %= 4; if(H[2] != [], H[2][i] %= 4));
  H[1][4,] %= 2; if(H[2] != [], H[2][4] %= 2);
  H[1] = mathnf(concat(H[1],matdiagonal([4,4,4,2])));
  if(H[2] != [], H[2] = coset_reduce(H[1],H[2]));
  H;
}

\\ Test whether G is a subgroup of H
issubgroup(G,H) =
{
  local(M);
  if(H[2] == [],
    if(G[2] == [],
      M = G[1];
      mathnf(concat(M, H[1])) == mathnf(H[1])
    , 0)
  , 
    M = G[1];
    if(G[2] != [], M = concat(M, G[2]~ - H[2]~));
    mathnf(concat(M, H[1])) == mathnf(H[1])
  );
}

\\ Test whether some conjugate of G is a subgroup of H
issubgroup_conj(G,H) =
{
  local M;
  if(H[2] == [],
    if(G[2] == [],
      mathnf(concat(G[1], H[1])) == mathnf(H[1])
    , 0)
  , 
    if(mathnf(concat(G[1], H[1])) == mathnf(H[1]),
      if(G[2] == [],
        1,
        M = concat(H[1], 2*matid(4));
        mathnf(concat(M, G[2]~-H[2]~)) == mathnf(M)
      )
    )
  );
}

\\ Swap two rows of the subgroup
swap_rows(H,i,j) =
{
  local(t);
  t = H[1][i,]; H[1][i,] = H[1][j,]; H[1][j,] = t;
  if(H[2] != [], t=H[2][i]; H[2][i] = H[2][j]; H[2][j] = t);
  H;
}

\\ Test whether some group equivalent to G is a subgroup of H
issubcase(G,H) =
{
  local(t, A, B, C);
  if(G[2] != [] && H[2] == [], return(0));
  if(matdet(G[1],1) % matdet(H[1],1) != 0, return(0));
  t = vector(24);
  t[1] = G;
  A = [-1,0,0,0;-1,1,0,0;-1,0,1,0;0,0,0,1];
  B = [1,-1,0,0;0,-1,0,0;0,-1,1,0;0,0,0,1];
  C = [1,0,-1,0;0,1,-1,0;0,0,-1,0;0,0,0,1];
  if(G[2] == [],
    t[7] = [A*G[1], []];
    t[13] = [B*G[1], []];
    t[19] = [C*G[1], []];
  ,
    t[7] = [A*G[1], (A*G[2]~)~];
    t[13] = [B*G[1], (B*G[2]~)~];
    t[19] = [C*G[1], (C*G[2]~)~];
  );
  for(i=0,3,
    t[6*i+2] = swap_rows(t[6*i+1],1,2);
    t[6*i+3] = swap_rows(t[6*i+2],2,3);
    t[6*i+4] = swap_rows(t[6*i+3],1,2);
    t[6*i+5] = swap_rows(t[6*i+4],2,3);
    t[6*i+6] = swap_rows(t[6*i+5],1,2);
  );
  for(i=1,24,if(issubgroup_conj(t[i],H), return(1)));
  0
}

\\ Return the element of t which comes first in the lex ordering 
lexmin(t) =
{
  local(n, j);
  n = matsize(t)[2];
  j = 1;
  for(i = 2, n, if(lex(t[i], t[j]) == -1, j = i));
  t[j];
}

\\ Return a canonical representative of the conjugacy class of H
conj_reduce(H) =
{
  for(i = 1, 3, H[1][i,] %= 4; if(H[2] != [], H[2][i] %= 4));
  H[1][4,] %= 2; if(H[2] != [], H[2][4] %= 2);
  H[1] = mathnf(concat(H[1],matdiagonal([4,4,4,2])));
  if(H[2] != [], H[2] = coset_reduce(mathnf(concat(H[1],2*matid(4))),H[2]));
  H;
}

\\ Return a canonical subgroup equivalent to H, by permuting a, b, c
make_unique1(H) =
{
  local(t);
  t = vector(6);
  t[1] = conj_reduce(H);
  t[2] = conj_reduce(swap_rows(t[1],1,2));
  t[3] = conj_reduce(swap_rows(t[2],2,3));
  t[4] = conj_reduce(swap_rows(t[3],1,2));
  t[5] = conj_reduce(swap_rows(t[4],2,3));
  t[6] = conj_reduce(swap_rows(t[5],1,2));
  lexmin(t);
}

\\ Return a canonical subgroup equivalent to H, up to projective changes
\\ of variables
make_unique(H) =
{
  local(A, B, C, t);
  A = [-1,0,0,0;-1,1,0,0;-1,0,1,0;0,0,0,1];
  B = [1,-1,0,0;0,-1,0,0;0,-1,1,0;0,0,0,1];
  C = [1,0,-1,0;0,1,-1,0;0,0,-1,0;0,0,0,1];
  t = vector(4);
  t[1] = make_unique1(H);
  if(H[2] == [],
    t[2] = make_unique1([A*H[1], []]);
    t[3] = make_unique1([B*H[1], []]);
    t[4] = make_unique1([C*H[1], []]);
  ,
    t[2] = make_unique1([A*H[1], (A*H[2]~)~]);
    t[3] = make_unique1([B*H[1], (B*H[2]~)~]);
    t[4] = make_unique1([C*H[1], (C*H[2]~)~]);
  );
  lexmin(t);
}