# Calculate the projective symmetry group with full pruning in
# the generic case.

with(LinearAlgebra):

d:=Dimension(V)[1]; n:=Dimension(V)[2];

#P:=HermitianTranspose(V).V:
P:=V:
P:=MatrixInverse(P,method=pseudo).P;
P:=evalf(P); evalf(P^2-P);
Pflag:=Vector(n,j->j):
# reorder P
rP:=P[]:
for j from 1 to n do for k from 1 to n do
  rP[j,k]:=P[Pflag[j],Pflag[k]];
end do; end do:
Q:=P[]:
tol:=10^(-3);

#oneproducts:=false;
#twoproducts:=false;

k:=0:
baseflag:=Vector(n,j->j):
kflags:={baseflag[]}:
#newflag:=Vector(n+1);
for k from 1 to n do
  oldflags:=kflags:
  kflags:={}:
  for j in oldflags do
    for l from k to n do
      possible:=true;
      newflag:=j[]; newflag[k]:=j[l]; newflag[l]:=j[k];
      nf:=newflag:
    # Pruning
# 1-products
      if oneproducts then
      if abs(rP[k,k]-Q[nf[k],nf[k]]) > tol
        then possible:=false; 
      end if;
      end if;
# 2-products
      if twoproducts then
      for j1 from 1 to k-1 do
      if abs(rP[j1,k]*rP[k,j1]-Q[nf[j1],nf[k]]*Q[nf[k],nf[j1]]) > tol
        then possible:=false; break;
      end if;
      end do;
      end if;
# 3-products
      for j1 from 1 to k-1 do for j2 from j1+1 to k-1 do
      if abs(rP[j1,j2]*rP[j2,k]*rP[k,j1]-Q[nf[j1],nf[j2]]
*Q[nf[j2],nf[k]]*Q[nf[k],nf[j1]]) > tol
        then possible:=false; break;
      end if;
      end do; if not possible then break; end if; end do;
      if possible then
        kflags:=kflags union {newflag[]};
      end if;
    end do;
  end do;
  print(k,nops(kflags));
end do:

nops(kflags);