freeze;


/******************************************************************
* Calculating the automorphism group of a hyperbolic lattice      *
* See J. Opgenorth, "Dual Cones and the Voronoi Algorithm", [Opg] *
*     M.H. Mertens, "Automorphism Groups of Hyperbolic Lattices", *
*                                                         [Mert]  *
*                                                                 *
* Michael H. Mertens, July 2013                                   *
******************************************************************/



PerfectPointsFormat:=recformat<
  ConeTestVectors,
  Points,
  Lat,
  directions,
  Neighbours,
  Stabilizers,
  Connectors,
  AutGrp,
  Watson,
  Complete
>;


ZZ:=Integers();
QQ:=Rationals();
// RR:=RealField(20);


// First we have to be able to determine whether a given matrix
// is hyperbolic or not


function SylvesterType(A)
  // The Sylvester type of A is given as a tuple (p,-q), 
  // where p denotes the number of positive, q the number of 
  // negative eigenvalues of A.
  
  if not (IsSymmetric(A) and (Type(BaseRing(A)) eq RngInt or Type(BaseRing(A)) eq FldRat)) then
    error "A must be a real symmetric matrix";
  end if;    
  D,_,_:=OrthogonalizeGram(A);
  Spec:=[D[i,i] : i in [1..Nrows(A)]];
  p:=#[ s : s in Spec | s gt 0];
  q:=#[ t : t in Spec | t lt 0];
  return <p,-q>;
end function;


function IsHyperbolic( A )
  // Returns true, if a symmetric matrix is of Sylvester type 
  // (n-1,-1),  false otherwise 
  t:=SylvesterType(A);
  n:=NumberOfRows(A);
  return t eq <n-1,-1>;
end function;


function IsIntegralMatrix( A )
  if Type(BaseRing(A)) eq RngInt then 
    return true;
  elif Type(BaseRing(A)) eq FldRat then
    ADen:=[Denominator(a) : a in Eltseq(A)];
    l:=LCM(ADen);
    return l eq 1;
  else
    error "Can not determine whether argument is integral or not";
  end if;
end function;



is_even := function(A)
  return GCD(Diagonal(A)) mod 2 eq 0;
end function;


function split_odd(A)  //Nach Kneser, Satz 26.3 Seite 107
  n0:=NumberOfRows(A);
  if n0 eq 1 then 
    return A[1,1], 0; 
  end if;
  even:=true; 
  for j in [1..n0] do 
    if GF(2) ! A[j][j] ne 0 then
      i:=j; even:=false;
      break;
    end if; 
  end for;

  if even then print "Fehler"; end if;
  C := RMatrixSpace(Integers(8),n0-1,n0) ! 0;
  //u := Integers() ! (Integers(32) ! A[i,i]^-1);
  u:=A[i][i]; //=u^-1 mod 8

  for j in [1..n0-1] do
    if j lt i then
      C[j,j] := 1;
      C[j,i] := -A[j,i]*u ;
    else
      C[j,j+1] := 1;
      C[j,i] := -A[j+1,i]*u ;
    end if;
  end for;
    
  B := C*A*Transpose(C);
  even:=true;
  for j in [1..n0-1] do 
    if GF(2) ! B[j][j] ne 0 then
      i1:=j; even:=false;
      break;
    end if; 
  end for;

  if even then
    for j in [1..n0-1] do 
      for k in [j+1..n0-1] do 
        if GF(2) ! B[j][k]  ne 0 then
          j0:=j; k0:=k;
          break j;
        end if;    
      end for;
    end for;
    I:=MatrixRing(Integers(8),n0) ! 0; 

    for j in [1..n0] do
      I[1][j] := C[j0][j];
    end for;

    I[1][i] +:= 1; 

    for k in [2..n0] do
      for j in [1..n0] do
        I[k][j]:=C[k-1][j];
      end for;
    end for;
     
    A := I*A*Transpose(I);
    C := RMatrixSpace(Integers(8),n0-1,n0) ! 0;
    u:=A[1][1]; //=u^-1 mod 8

    for j in [1..n0-1] do
      C[j,j+1] := 1;
      C[j,1] := -A[j+1,1]*u ;
    end for;

    B := C*A*Transpose(C);
  end if;
  return u, B;
end function;


function twoDiagonalization(A) 
  n:=NumberOfRows(A);
  A8:=MatrixRing(Integers(8),n) ! A;
  AD:=Parent(A) ! 0;
  B:=A8;

  for i in [1..n] do  
    u,B := split_odd(B); 
    AD[i][i]:=Integers() ! u; 
  end for;

  return AD;
end function;


function two_adic_genus_symbol(A, bound)
  error if not IsSymmetric(A), "A is not a symmetric matrix.";
  p:=2;
  n := NumberOfRows(A);
  F_n := MatrixAlgebra(GF(p), n); // F_p^{n \times n}
  V_n := VectorSpace(GF(p), n); // F_p^n
  Z_n := MatrixAlgebra(Integers(), n); // Z^{n \times n}
  Z8_n := MatrixAlgebra(Integers(8), n); // (Z/8Z)^{n \times n}
  A := ChangeRing(A, Integers()); // allows A to be formally 
                                  // defined over Rationals() 
				  // upon calling this routine
  // Divide A by the highest power of p dividing all its elements:
  
  g := Content(A);  // gcd of all elements of A
  v := Valuation(g, p);
  A := A div (p^v);
  A_p := F_n!A; // A over F_2
  C := Kernel(A_p);

  // Otherwise, we split off the Jordan constituent of the highest power of 2:  
  if Dimension(C) eq 0 then
    e0 := KroneckerSymbol(Determinant(A), 2);
    if is_even(A) then
      s := 0; // "even"
      o := 0; // oddity not applicable for even lattices, stored as 0
    else
      s := 1; // "odd"
      // An odd matrix can be diagonalized over the 2-adic integers. We can then determine
      // the oddity by computing the trace of Z/8Z:
      AD := twoDiagonalization(A);
      o := Integers(8)!Trace(AD); 
    end if;
    return [[v, n, e0, s, Integers()!o]];						    
  else					
    T := ChangeRing(BasisMatrix(C), Integers());
    B_p := BasisMatrix(Complement(V_n, C));
    B := ChangeRing(B_p, Integers());
    AB := B*A*Transpose(B);
    ABI:=Adjoint(AB);
    det:=Determinant(AB);
    deti:=Integers() ! ((Integers(p^bound) ! det) ^-1);
    X:=B*A*Transpose(T);
    KOR:= deti*Transpose(X)*ABI*X;
    J2 := T * A * Transpose(T)-KOR;
    n0 := Rank(B_p);

    if is_even(AB) then
      s:=0;
      o:=0;
    else
      s:=1;
      //to compute the oddity diagonalize the matrix over the p-adic integers and 
      //then compute the trace over Z mod 8
      AD:=twoDiagonalization(AB);
      o:=Integers(8) ! Trace(AD);
    end if;

    return [[v,n0,KroneckerSymbol(Determinant(AB),p),s,Integers()!o]] cat $$(p^v*J2,bound);
  end if;
end function;


function p_adic_genus_symbol(A, p,bound)
  error if not IsSymmetric(A), "A is not a symmetric matrix";
  error if not IsPrime(p), "p is not prime";

  // parameter "bound": should be >= 3 + p-valuation of the maximal elementary divisor of 
  // (passed as a parameter since p_adic_genus is usually called multiple times for the
  // same matrix and calculating elementary divisors is expensive)

  if p eq 2 then 
    return two_adic_genus_symbol(A, bound); 
  end if;
  
  n := NumberOfRows(A);
  F_n := MatrixAlgebra(GF(p), n); // F_p^{n \times n}
  V_n := VectorSpace(GF(p), n); // F_p^n
  Z_n := MatrixAlgebra(Integers(), n); // Z^{n \times n}
  A := ChangeRing(A, Integers()); // allows A to be formally defined over Rationals() 
                                  // upon calling this routine
  v := Valuation(Content(A), p);
  A := A div (p^v);
  A_p := F_n!A; // A over F_p 
  C := Kernel(A_p);
  // If the kernel of A mod p is trivial, A consists only of a zero-valuated Jordan constituent.
  // Otherwise, we split off the Jordan constituent of the highest p-power:

  if Dimension(C) eq 0 then
    e0 := KroneckerSymbol(Determinant(A), p);
    // printf "Only one Jordan constituent, stopping iteration. p^v=%o.\n", p^v;
    return [[v, n, e0]];
  else
    T := ChangeRing(BasisMatrix(C), Integers());
    B_p := BasisMatrix(Complement(V_n, C));
    A_p := B_p * A_p * Transpose(B_p);
    B := ChangeRing(B_p, Integers());
    AB := B*A*Transpose(B);
    ABI:=Adjoint(AB); 
    det:=Determinant(AB);
    deti:=Integers() ! ((Integers(p^bound) ! det) ^-1);
    X:=B*A*Transpose(T); 
    KOR:= deti*Transpose(X)*ABI*X;
    J2 := T * A * Transpose(T)-KOR;
    n0 := Rank(A_p);
    if IsSquare(Determinant(A_p)) then
      e0 := 1;
    else
      e0 := -1;
    end if;
    // printf "Iterating with p^v=%o, matrix A =\n%o\n,
    //  matrix A_p =\n%o\nmatrix J2 = \n%o\nbasis of kernel = \n%o\nbasis
    //  of complement = \n%o\n\n",p^v,A,A_p,J2,BasisMatrix(C),
   
    return [[v, n0, e0]] cat $$(p^v * J2, p,bound);
  end if;
end function;


function genus_symbol(A)
  if not(A in RMatrixSpace(IntegerRing(), NumberOfRows(A), NumberOfColumns(A)))
     or (NumberOfRows(A) ne NumberOfColumns(A)) or not(IsSymmetric(A)) then    
    error "Argument gram_matrix must be a symmetric positive definite integral matrix.";
  end if;

  divisors := [d[1]: d in Factorization(Determinant(A))| d[1] gt 2];
  g := ElementaryDivisors(A);
  h := Factorization(g[#g]);
  sym := [<2, p_adic_genus_symbol(A, 2, Valuation(g[#g], 2)+3)>];
 
  for t in h do
    if t[1] ne 2 then
      Append(~sym, <t[1], p_adic_genus_symbol(A, t[1], t[2])>);
    end if;
  end for;

  return sym;
end function;



// Watson process

function WatsonGram(A, p)
  // Valuation(Determinant(A), p);
  D:= p * ChangeRing(A, QQ)^-1;
  L:= Lattice(D);
  P:= (L meet StandardLattice(Ncols(A)));
  T:= Matrix(QQ, BasisMatrix(P));
  AA:= p^-1*T * A * Transpose(T);
  // Valuation(Determinant(AA), p);
  return ChangeRing(AA,ZZ), T;
end function;


function WatsonProcess(A)
  n:=Nrows(A);
  B:=A;
  C:=IdentityMatrix(QQ,n);
  sym:=genus_symbol(A);

  for p in sym do
    q:=p[1];
    a:=Max([p[2][i][1] : i in [1..#p[2]]]);
    while a ge 2 do
      B,T:=WatsonGram(B,q);
      C:=T*C;
      a-:=2;
    end while;
    if Valuation(Determinant(B),q) gt Valuation(Determinant(A),q) then
      B,T:=WatsonGram(B,q);
      C:=T*C;
    end if;
  end for;

  return B, C;
end function;


function Get_y( A : record:=rec< PerfectPointsFormat | > )
  // finds a vector y s.t. yA^-1y^tr<0 for a rational matrix A
  if not assigned record`ConeTestVectors then
    A:=ChangeRing(A,QQ);
    ok, Ainv:=IsInvertible(A);
    if not ok then
      error "Error: A must be invertible.";
    end if;
    D,T,_:=OrthogonalizeGram(A);
    SpecNeg:=[ i : i in [1..Nrows(A)] | D[i,i] lt 0];
    if SpecNeg eq [] then
      error "Error: There is no such y.";
    end if;
    y1:=Rows(Transpose(T)^-1)[SpecNeg[1]];
    y2:=Rows(T)[SpecNeg[1]];
    return [y1,y2];
  else return record`ConeTestVectors;
  end if;
end function;


function xAxtr(x,A)
  x:=ChangeRing(x,QQ);
  A:=ChangeRing(A,QQ);
  return (Vector(x) * A * Matrix(1, Eltseq(x)))[1];
end function;



function xperp(x)
  if Type(BaseRing(x)) eq FldRat then
    xDen:=[Denominator(xi) : xi in Eltseq(x)];
    r:=LCM(xDen);
    x:=ChangeRing(r*x, ZZ);
  elif Type(BaseRing(x)) eq RngInt then
    x:=Vector(RemoveRowContents(Matrix(x)));
  else
    error "Argument must be a rational vector.";
  end if;
  H:= Basis(NullSpace(Matrix(ZZ, 1,Eltseq(x))));
  assert Max(ElementaryDivisors(Matrix(H))) eq 1;
  return H;
end function;


// This is an exact copy of the corresponding function in the 
// normalizer package:
// G acts through h: SxG->S. We find generators for the stabilizer 
// in G of some X \in S. G is usually infinite! But the orbit of X is not.

function MyStabGens(G, h, X : IsFinite:=false)
  n:= Ngens(G);
  Orbit:= {@ X @};
  Path:= [ Id(G) ];
  Stab:= [];
  if IsFinite then 
    S:=sub< G | Stab >;
    g:=#G;
    s:=1; 
  end if;
  i:= 0;
  complete:= false;
  while i lt #Orbit and not complete do
    i+:= 1;
    for j in [1..n] do
      x:= Path[i] * G.j;
      Xx:= h(X, x);
      idx:= Index(Orbit, Xx);
      if idx eq 0 then
        Include(~Orbit, Xx);
        Append(~Path, x);
        if IsFinite then
          complete:=2*s*#Orbit gt g;
        end if;
      else
        x:= x * Path[idx]^-1;
        assert h(X, x) eq X;
        if not IsOne(x) and x notin Stab and x^-1 notin Stab then
          Append(~Stab, x);
          if IsFinite then
            S:=sub< G | Stab>;
            s:=#S;
            complete:= 2*s*#Orbit gt g;
          end if; 
        end if;
      end if;
    end for;
  end while;

  return Stab;
end function;



// now we start:

/************************************************************************
* For given integer x in the cone V_1^{>0} one can calculate the set    *
* of D-minimal vectors for x by the following method:                   *
* A^-1 defines a positive definite scalar product on the  hyperplane H  *
* perpendicular (w.r.t the standard scalar product) to x. Define the    *
* lattice  L= H\cap 1/-N(x)ZZ^n endowed with the scalar product induced *
* by A^-1. Here N(x)=xAx^tr. For lambda=2/-N(x),3/-N(x),.. calculate    *
* the set M_lambda of lattice points y with yA^-1y^tr leq -lambda^2N(x) *
* As soon as this set is not empty (which is the case at least for      *
* lambda=1), the set {-lambda xA + y | y in S_L}\cap ZZ^n is the desired*
* answer.                                                               *
************************************************************************/

function DMinimalVectors(x, A : record:=rec< PerfectPointsFormat | >, test:=false )
	record`ConeTestVectors:=Get_y(A : record:=record);
  y1:=record`ConeTestVectors[1];
  y2:=record`ConeTestVectors[2];
  if not (Type(BaseRing(x)) eq RngInt or 
          Type(BaseRing(x)) eq FldRat) then
    error "First argument must be a rational vector";
  elif not (IsIntegralMatrix(A) and IsHyperbolic(A)) then
    error "Second argument must be an integral hyperbolic matrix.";
  elif not (xAxtr(x,A) lt 0 and 
           (ChangeRing(x,QQ),ChangeRing(y1,QQ)) gt 0) then
    error "First argument does not belong to the dual cone.";
  end if;
  // First rescale x s.t. it becomes an integer vector
  // Note that the D-minimum of a vector x is rescaled by lambda
  // and the set of D-minimal vectors remains the same if one replaces
  // x by lambda*x, lambda > 0
  
  if Type(BaseRing(x)) eq FldRat then
    xSeq:=Eltseq(x);
    xDen:=[Denominator(xi) : xi in xSeq];
    r:=Lcm(xDen);
  else
    r:=1;
  end if;
  x:=ChangeRing(r*x, ZZ);
  t:=Gcd(Eltseq(x));
  x:=Vector(x div t);
  y2:=ChangeRing(y2,QQ);
  H:=xperp(x);
  x:=ChangeRing(x, QQ);
  A:=ChangeRing(A, QQ);
  Ainv:=A^-1;
  N:=-xAxtr(x,A);
  B:=1/N*ChangeRing(Matrix(H),QQ);
  L:=LatticeWithGram(B*Ainv*Transpose(B));
  ZZn:=StandardLattice(Nrows(A));

  Long:= N ge 3000 or Determinant(GramMatrix(L)) ge 10000;
  if Long then
    N, Determinant(GramMatrix(L));
    print "Warning: This might take a while.";
  end if;

  // Of course Ainv is not positive definite on RR^n, 
  // but it is indeed on H.
  l:=0;
  S:=[Vector(ChangeRing(L!0,QQ)*B)];
  M:=[];
  while M eq [] do
    l+:=1;
    lambda1:=(l-1)/N;
    lambda2:=l/N;
    Snew:=lambda1 eq 0 select ShortVectors(L, lambda2^2*N) else
                              ShortVectors(L,lambda1^2*N,lambda2^2*N);
    S1:=[Vector(ChangeRing(s[1],QQ)*B) : s in Snew];
    S2:=[-s : s in S1];
    S cat:=S1 cat S2/* cat [Vector(ChangeRing(L!0,QQ)*B)]*/;
    if test then
      Stest:=ShortVectors(L, lambda2^2*N);
      Stest1:=[Vector(ChangeRing(s[1],QQ)*B) : s in Stest];
      Stest2:=[-s : s in Stest1];
      Stest:=Stest1 cat Stest2 cat [Vector(ChangeRing(L!0,QQ)*B)];
      assert Seqset(S) eq Seqset(Stest);
    end if;

    if Long then
      #S;
    end if;

    for i in [1..#S] do
      y:=S[i];
      d:=y-lambda2*x*A;
      if d in ZZn and (y2,ChangeRing(d,QQ)) ge 0 then
        Append(~M, d);
      end if;
    end for;
  end while;
 
  mu:=(ChangeRing(x,QQ),ChangeRing(M[1],QQ));
  return M, t*mu/r;
end function;



function DMinimalVectors2(x, A : record:=rec< PerfectPointsFormat | > ) 
    	 record`ConeTestVectors:=Get_y(A : record:=record);
  y1:=record`ConeTestVectors[1];
  y2:=record`ConeTestVectors[2];
  if not (Type(BaseRing(x)) eq RngInt or 
          Type(BaseRing(x)) eq FldRat) then
    error "First argument must be a rational vector";
  elif not (IsIntegralMatrix(A) and IsHyperbolic(A)) then
    error "Second argument must be an integral hyperbolic matrix.";
  elif not (xAxtr(x,A) lt 0 and 
           (ChangeRing(x,QQ),ChangeRing(y1,QQ)) gt 0) then
    error "First argument does not belong to the dual cone.";
  end if;
  // First rescale x s.t. it becomes an integer vector
  // Note that the D-minimum of a vector x is rescaled by lambda
  // and the set of D-minimal vectors remains the same if one replaces
  // x by lambda*x, lambda > 0
  
  if Type(BaseRing(x)) eq FldRat then
    xSeq:=Eltseq(x);
    xDen:=[Denominator(xi) : xi in xSeq];
    r:=Lcm(xDen);
  else
    r:=1;
  end if;
  x:=ChangeRing(r*x, ZZ);
  t:=Gcd(Eltseq(x));
  x:=Vector(x div t);
  y2:=ChangeRing(y2,QQ);
  H:=xperp(x);
  x:=ChangeRing(x, QQ);
  A:=ChangeRing(A, QQ);
  Ainv:=A^-1;
  N:=-xAxtr(x,A);
  B:=1/N*ChangeRing(Matrix(H),QQ);
  L:=LatticeWithGram(B*Ainv*Transpose(B));
  ZZn:=StandardLattice(Nrows(A));

  Long:= N ge 3000 or Determinant(GramMatrix(L)) ge 10000;
  if Long then
    N, Determinant(GramMatrix(L));
    print "Warning: This might take a while.";
  end if;

  // Of course Ainv is not positive definite on RR^n, 
  // but it is indeed on H.
  l:=0;
  M:=[];
  while M eq [] do
    l+:=1;
    lambda2:=l/N;
    Stest:=ShortVectors(L, lambda2^2*N);
    Stest1:=[Vector(ChangeRing(s[1],QQ)*B) : s in Stest];
    Stest2:=[-s : s in Stest1];
    S:=Stest1 cat Stest2 cat [Vector(ChangeRing(L!0,QQ)*B)];
    
    if Long then
      #S;
    end if;

    for i in [1..#S] do
      y:=S[i];
      d:=y-lambda2*x*A;
      if d in ZZn and (y2,ChangeRing(d,QQ)) ge 0 then
        Append(~M, d);
      end if;
    end for;
  end while;
 
  mu:=(ChangeRing(x,QQ),ChangeRing(M[1],QQ));
  return M, mu/r*t;
end function;



// We do almost the same to calculate the set of D-vectors with
// (x, d) le C, C > 0.
function DShortVectors(x, A, C : 
			 record:=rec<PerfectPointsFormat | >, test:=false)
  if not C gt 0 then
    error "Bound must be positive.";
  end if;   
  
  record`ConeTestVectors:=Get_y(A : record:=record);
  y1:=record`ConeTestVectors[1];
  y2:=record`ConeTestVectors[2];

  if not (Type(BaseRing(x)) eq RngInt or 
          Type(BaseRing(x)) eq FldRat) then
    error "First argument must be a rational vector";
  elif not (Type(BaseRing(A)) eq RngInt and IsHyperbolic(A)) then
    error "Second argument must be an integral hyperbolic matrix.";
  elif not (xAxtr(x,A) lt 0 and 
           (ChangeRing(x,QQ),ChangeRing(y1,QQ)) gt 0) then
    error "First argument does not belong to the dual cone.";
  end if;
  // First rescale x s.t. it becomes an integer vector
  
  if Type(BaseRing(x)) eq FldRat then
    xSeq:=Eltseq(x);
    xDen:=[Denominator(xi) : xi in xSeq];
    r:=Lcm(xDen);
  else
    r:=1;
  end if;
  x:=ChangeRing(r*x, ZZ);
  t:=Gcd(Eltseq(x));
  x:=Vector(x div t);
  H:=xperp(x);
  x:=ChangeRing(x, QQ);
  A:=ChangeRing(A, QQ);
  Ainv:=A^-1;
  N:=-xAxtr(x,A);
  B:=1/N*ChangeRing(Matrix(H),QQ);
  L:=LatticeWithGram(B*Ainv*Transpose(B));
  ZZn:=StandardLattice(Nrows(A));
  // Of course Ainv is not positive definite on RR^n, 
  // but it is indeed on H.
  l:=0;
  M:=[];
  mu:=0;
  S:=[Vector(ChangeRing(L!0,QQ)*B)];
  while mu lt r*C/t do
    l+:=1;
    lambda1:=(l-1)/N;
    lambda2:=l/N;
    Snew:=lambda1 eq 0 select ShortVectors(L, lambda2^2*N) else
                              ShortVectors(L,lambda1^2*N,lambda2^2*N);
    S1:=[Vector(ChangeRing(s[1],QQ)*B) : s in Snew];
    S2:=[-s : s in S1];
    S cat:=S1 cat S2 /*cat [Vector(ChangeRing(L!0,QQ)*B)]*/;
    if test then
      Stest:=ShortVectors(L, lambda2^2*N);
      Stest1:=[Vector(ChangeRing(s[1],QQ)*B) : s in Stest];
      Stest2:=[-s : s in Stest1];
      Stest:=Stest1 cat Stest2 cat [Vector(ChangeRing(L!0,QQ)*B)];
      assert Seqset(S) eq Seqset(Stest);
    end if;
    
    for i in [1..#S] do
      y:=S[i];
      d:=y-lambda2*x*A;
      if d in ZZn and (ChangeRing(y2,QQ),d) ge 0 then
        Append(~M, d);
        mu2:=(ChangeRing(x,QQ),ChangeRing(d,QQ));
        if mu2 gt mu then
          mu:=mu2;
        end if;
      end if;
    end for;
  end while;
  return M;
end function;


function DShortVectors2(x, A, C : 
         record:=rec<PerfectPointsFormat | >)   
  if not C gt 0 then
    error "Bound must be positive.";
  end if;   
  
  record`ConeTestVectors:=Get_y(A : record:=record);
  y1:=record`ConeTestVectors[1];
  y2:=record`ConeTestVectors[2];

  if not (Type(BaseRing(x)) eq RngInt or 
          Type(BaseRing(x)) eq FldRat) then
    error "First argument must be a rational vector";
  elif not (Type(BaseRing(A)) eq RngInt and IsHyperbolic(A)) then
    error "Second argument must be an integral hyperbolic matrix.";
  elif not (xAxtr(x,A) lt 0 and 
           (ChangeRing(x,QQ),ChangeRing(y1,QQ)) gt 0) then
    error "First argument does not belong to the dual cone.";
  end if;
  // First rescale x s.t. it becomes an integer vector
  
  if Type(BaseRing(x)) eq FldRat then
    xSeq:=Eltseq(x);
    xDen:=[Denominator(xi) : xi in xSeq];
    r:=Lcm(xDen);
  else
    r:=1;
  end if;
  x:=ChangeRing(r*x, ZZ);
  t:=Gcd(Eltseq(x));
  x:=Vector(x div t);
  H:=xperp(x);
  x:=ChangeRing(x, QQ);
  A:=ChangeRing(A, QQ);
  Ainv:=A^-1;
  N:=-xAxtr(x,A);
  B:=1/N*ChangeRing(Matrix(H),QQ);
  L:=LatticeWithGram(B*Ainv*Transpose(B));
  ZZn:=StandardLattice(Nrows(A));
  // Of course Ainv is not positive definite on RR^n, 
  // but it is indeed on H.
  l:=0;
  M:=[];
  mu:=0;
  while mu lt r*C/t do
    l+:=1;
    lambda2:=l/N;
    Stest:=ShortVectors(L, lambda2^2*N);
    Stest1:=[Vector(ChangeRing(s[1],QQ)*B) : s in Stest];
    Stest2:=[-s : s in Stest1];
    S:=Stest1 cat Stest2 cat [Vector(ChangeRing(L!0,QQ)*B)];
    
    for i in [1..#S] do
      y:=S[i];
      d:=y-lambda2*x*A;
      if d in ZZn and (ChangeRing(y2,QQ),d) ge 0 then
        Append(~M, d);
        mu2:=(ChangeRing(x,QQ),ChangeRing(d,QQ));
        if mu2 gt mu then
          mu:=mu2;
        end if;
      end if;
    end for;
  end while;
  return M;
end function;



function GetConeVector( A : record:=rec< PerfectPointsFormat | > )
  if not (IsHyperbolic(A) and IsIntegralMatrix(A)) then
    error "Argument must be an integral hyperbolic matrix";
  end if;
  record`ConeTestVectors:=Get_y(A : record:=record);
  y1:=record`ConeTestVectors[1];  

  n:=Nrows(A);
  ZZn:=StandardLattice(n);
  Bound:=0;
  x:=ZZn!0;
  null:=x;
  
  while x eq null do
    Bound+:=1;
    S:=Bound eq 1 select [Vector(s[1]) : s in ShortVectors(ZZn,Bound)] else [Vector(s[1]) : s in ShortVectors(ZZn, Bound-1,Bound)];
    S:=S cat [-s : s in S];
    for y in S do
      Ny:=xAxtr(y,A);
      if Ny lt 0 then
        sign:=1;
        if (ChangeRing(y,QQ),ChangeRing(y1,QQ)) lt 0 then
          sign:=-1;
        end if;
	//We don't have to take the first best point in the cone, maybe we can find a (close to) perfect one just by trial and error:
	if Ny ge -500 then
          M,mu:=DMinimalVectors(sign*y,A);
	else
	  M:=[];
	end if;
		if #M ge n or Bound gt 10 then
		  x:=sign*y;
        end if;
      end if;
    end for;      
  end while;

  return ChangeRing(x, ZZ);
end function;



function IsDPerfect( x, A : M:=[], 
         record:=rec<PerfectPointsFormat | > )
  if M eq [] then
    M,_:=DMinimalVectors(x, A : record:=record);
  end if;
  V:=VectorSpace(QQ, Ncols(x));
  U:=sub< V | M >;
  return Dimension(U) eq Dimension(V); 
end function;



function Directions( x, A : M:=[], 
         record:=rec<PerfectPointsFormat | >)
  if not IsDPerfect(x, A : record:=record) then
    error "First argument must be a D-perfect vector.";
  end if;
  n:=Nrows(A);
  if M eq [] then
    M, mu:=DMinimalVectors(x,A);
  else
    mu:=(ChangeRing(x,QQ),ChangeRing(M[1],QQ));
  end if;
  T:=ToricLattice(n);
  M:=[T!v : v in M];
  C:=Cone(M);
  dirs:={Vector(Eltseq(r)) : r in Inequalities(C)};
  directions:={ };
  // r is now a direction of x, but maybe r is blind:
  for r in dirs do
    y1:=Get_y(A : record:=record)[1];
    if not (xAxtr(r,A) le 0 and 
	       (ChangeRing(r,QQ),ChangeRing(y1,QQ)) ge 0) then
        // we want the directions to be integral vectors
          if Type(BaseRing(r)) eq FldRat then
            rDen:=[Denominator(ri) : ri in Eltseq(r)];
            l:=Lcm(rDen);
          else 
            l:= 1;
          end if;
          r1:=ChangeRing(l*r, ZZ);
          d:=Gcd(Eltseq(r1));
          r1:=Vector(r1 div d);
          Include(~directions, r1);
        else
          print "Found blind direction:"; 
          r;
      end if;
  end for;
  return SetToSequence(directions);
end function;


// This function adapts the function vor_neighbor from the 
// normalizer package
function VoronoiNeighbour1(A, x, r : Limit:=10^6, 
         record:=rec<PerfectPointsFormat | >)
  record`ConeTestVectors:=Get_y(A : record:=record);
  Mx,mux:=DMinimalVectors2(x,A : record:=record);
  // rescale x s.t. it becomes integral
  if Type(BaseRing(x)) eq FldRat then
    xSeq:=Eltseq(x);
    xDen:=[Denominator(xi) : xi in xSeq];
    l:=Lcm(xDen);
  else
    l:=1;
  end if;
  if Type(BaseRing(r)) eq FldRat then
    rSeq:=Eltseq(r);
    rDen:=[Denominator(ri) : ri in rSeq];
    lr:=Lcm(rDen);
  else 
    lr:=1;
  end if;
  
  t1:=Get_y(A : record:=record)[1];
  // Just to make sure nothing awkward happens...
  x:=ChangeRing(x*l, ZZ);
  mux:=ZZ!(mux*l);
  r:=ChangeRing(lr*r, ZZ);
  xr:= ChangeRing(x,QQ);
  rr:= ChangeRing(r,QQ);

  // Calculate a vector y=lm*x+rm*r which lies in the cone and 
  // whose D-minimum is smaller than mux
 
  lu:= 1; ru:= 0;
  lo:= 0; ro:= 1;
  found:= false;
  counter:=0;
  repeat
    lm:=lu + lo;
    rm:=ru + ro;
    y:=lm*xr + rm*rr;
    if xAxtr(y,A) lt 0 and (y,ChangeRing(t1,QQ)) gt 0 then
      My, muy:=DMinimalVectors2(y,A);
      if muy lt lm*mux then
        found:= true;
      else
        lu:= lm; ru:= rm;
      end if;
    else 
      lo:= lm; ro:= rm;
    end if;
    counter +:=1;
  until found or counter ge Limit; 

  if counter ge Limit then // This should never happen
    error "After",Limit," trials there was no suiting point.";
  end if;

  M:=DShortVectors2(y,A,lm*mux : record:=record);
  M:= [ x: x in M | x notin Mx ];
  assert  forall{d: d in M | (Vector(r)*Matrix(1,Eltseq(d)))[1] lt 0 };
  lambda:=Min([ (mux - ZZ!(Vector(x)*Matrix(1,Eltseq(d)))[1])
          /(Vector(r)*Matrix(1,Eltseq(d)))[1] : d in M ]);
  p:=Numerator(lambda);
  q:=Denominator(lambda);
  y:=q*x+p*r;
  My, muy:=DMinimalVectors2(y,A : record:=record);
  
  if IsDPerfect(x,A : M:=Mx, record:=record) then
    assert IsDPerfect(y,A : M:=My, record:=record);
  end if;
  return y, My,muy;
end function;


// This function is supposed to do the same as VoronoiNeighbour1, 
// but uses a slightly different method

function VoronoiNeighbour2(A,x,r : Limit:=10^3, 
         record:=rec<PerfectPointsFormat | >,
   ShortVectors:=[], test:=false)
  record`ConeTestVectors:=Get_y(A : record:=record);
  if xAxtr(x,A) ge 0 or (ChangeRing(x,QQ),record`ConeTestVectors[1]) le 0 then
    error "2nd argument does not belong to the cone.";
  end if;
  // rescale x s.t. it becomes integral
  if Type(BaseRing(x)) eq FldRat then
    xSeq:=Eltseq(x);
    xDen:=[Denominator(xi) : xi in xSeq];
    l:=Lcm(xDen);
  else
    l:=1;
  end if;
 
  if Type(BaseRing(r)) eq FldRat then
    rSeq:=Eltseq(r);
    rDen:=[Denominator(ri) : ri in rSeq];
    lr:=Lcm(rDen);
  else
    lr:=1;
  end if;

  t1:=record`ConeTestVectors[1];
  yA:= ChangeRing(record`ConeTestVectors[1], QQ);
  
  // Just to make sure nothing awkward happens...
  x:=ChangeRing(l*x, ZZ);
  tx:=Gcd(Eltseq(x));
  x:=Vector(x div tx);
  Mx,mux:=DMinimalVectors(x,A : test:=test);
  r:=ChangeRing(lr*r, ZZ);
  tr:=Gcd(Eltseq(r));
  r:=Vector(r div tr);
  xr:= ChangeRing(x,QQ);
  rr:= ChangeRing(r,QQ);

  // Calculate a vector y=lm*x+rm*r which lies in the cone and 
  // whose D-minimum is smaller than lm*mux
  lu:= 1; ru:= 0;
  lo:= 0; ro:= 1;
  lm:=0;
  found:= false;
  counter:=0;
  repeat
    counter+:=1;
    lm:=lu + lo;
    rm:=ru + ro;
    y:=lm*xr + rm*rr;
    if xAxtr(y,A) lt 0 and (y,yA) gt 0 then
    My, muy:=DMinimalVectors(y,A : record:=record, test:=test);
      if muy lt lm*mux then
        found:= true;
      else
        lu:= lm; ru:= rm;
      end if;
    else 
      lo:= lm; ro:= rm;
    end if;
  until found or counter ge Limit;
  //print "Found neighbour after ",counter,"trials.";
  if counter ge Limit then // This should never happen
    error "After",Limit," trials there was no suiting point.";
  end if;
  //t0:=Realtime();
  M:=DShortVectors(y,A,lm*mux : record:=record , test:=test);
  //t1:=Realtime(t0);
  //print "Needed",t1,"seconds to compute DShortvectors of y with norm",xAxtr(y,A)," of length le",lm*mux;
  M:= [ x: x in M | x notin Mx ];
  assert  forall{d: d in M | (Vector(r)*Matrix(1,Eltseq(d)))[1] lt 0 };
  lambda:=Min([ (mux - ZZ!(Vector(x)*Matrix(1,Eltseq(d)))[1])
          /(Vector(r)*Matrix(1,Eltseq(d)))[1] : d in M ]);
  p:=Numerator(lambda);
  q:=Denominator(lambda);
  y:=q*x+p*r;
  d:=Gcd(Eltseq(y));
  y:=Vector(y div d);
  My, muy:=DMinimalVectors(y,A : record:=record, test:=test);
  
  if IsDPerfect(x,A : M:=Mx, record:=record) then
    assert IsDPerfect(y,A : M:=My, record:=record);
  end if;

  return y, My,muy;
end function;


// L : Basis of Q^n
function action(L, g)
  B:= Matrix([ ChangeRing(l,QQ) * ChangeRing(g, QQ) : l in L ]);
  d:= Denominator(B);
  M:= HermiteForm( Matrix(ZZ, d*B) );
  return [ ChangeRing(M[i], QQ) / d : i in [1..Nrows(M)] ];
end function;



function PerfectStab(x,A)
  if Type(BaseRing(x)) eq FldRat then
    xDen:=[Denominator(xi) : xi in Eltseq(x)];
    l:=Lcm(xDen);
  else
    l:=1;
  end if;
  
  x:=ChangeRing(l*x,ZZ);
  n:=Nrows(A);
  ZZn:=Basis(StandardLattice(n));
  ZZn:=[ChangeRing(Vector(l),QQ) : l in ZZn];
  H:=Matrix(xperp(x*A));
  C:=H*A*Transpose(H);
  L:=LatticeWithGram(C);
  Aut:=GeneratorsSequence(AutomorphismGroup(L));
  Aut:=[DiagonalJoin(g, Matrix(1,[1])) : g in Aut ];
  B:=ChangeRing(VerticalJoin(H, Matrix(x)), QQ);
  GxGens:=[B^-1*g*B : g in Aut ];
  Gx:=sub< GL(n,QQ) | GxGens>;

  // Gx is not a unimodular group in general, but the stabilizer
  // Sx of a D-perfect point can be found as a subgroup via orbit 
  // algorithm. Sometimes, there is a shortcut if the generators of 
  // Gx are all integral.

  if forall{ g : g in GxGens | IsIntegralMatrix(g) } then
    Sx:=sub< GL(n,ZZ) | GxGens >;
    Sxrat:=Gx;
  else    
   SxGens:=MyStabGens(Gx, action, ZZn : IsFinite:=true);
   assert forall{ s : s in SxGens | s*ChangeRing(A,QQ)*Transpose(s)
                  eq ChangeRing(A,QQ)};
   Sx:=sub< GL(n,ZZ) | SxGens >;
   Sxrat:=sub< Gx | SxGens >;
  end if;

  return Sx, Transversal(Gx, Sxrat), LLLGramMatrix(L);
end function;


function PerfectConnector(x, y, A : Trans:={@ @})
  if Type(BaseRing(x)) eq FldRat or 
     Type(BaseRing(y)) eq FldRat then
    xDen:=[Denominator(xi) : xi in Eltseq(x)];
    lx:=Lcm(xDen);
    yDen:=[Denominator(yi) : yi in Eltseq(y)];
    ly:=Lcm(yDen);
  else
    lx:=1;
    ly:=1;
  end if;

  x:=ChangeRing(lx*x,ZZ);
  y:=ChangeRing(ly*y,ZZ);

  // rescale x and y s.t. xAx^tr=yAy^tr
  xAx:=-ZZ!(xAxtr(x,A));
  yAy:=-ZZ!(xAxtr(y,A));
  l:=Gcd(xAx,yAy);
  okx,lx:=IsSquare(yAy div l);
  oky,ly:=IsSquare(xAx div l);
  if not (okx and oky) then
    return false,_;
  end if;

  x:=(ZZ!lx)*x;
  y:=(ZZ!ly)*y;

  Hx:=Matrix(xperp(x*A));
  Cx:=Hx*A*Transpose(Hx);
  Lx:=LatticeWithGram(Cx);

  Hy:=Matrix(xperp(y*A));
  Cy:=Hy*A*Transpose(Hy);
  Ly:=LatticeWithGram(Cy);

  ok,T := IsIsometric(Lx,Ly);

  if not ok then 
    return false, _; 
  end if;
  
  T:=DiagonalJoin(ChangeRing(T,QQ),Matrix(QQ,1,[1]));
  Hx:=ChangeRing(VerticalJoin(Hx,Matrix(x)),QQ);
  Hy:=ChangeRing(VerticalJoin(Hy,Matrix(y)),QQ);
  TT:=Hy^-1*T*Hx;
  assert TT*ChangeRing(A,QQ)*Transpose(TT) eq ChangeRing(A,QQ);
  Sx:=Trans;
  ok:=exists(w){TT*ChangeRing(g,QQ) : g in Sx |
                IsIntegralMatrix(TT*ChangeRing(g,QQ))};  
  if ok then 
    w:=ChangeRing(w,ZZ);
    assert y*w eq x;
    assert w*A*Transpose(w) eq A;
    return ok, w;
  else 
    print "No integral isometry between ", x,y;
    return ok, _;
  end if;
end function;


// returns an integral D-perfect vector and 
// a list of its D-minimal vectors.
function FirstPerfect( A : method:=2, Limit:=10^6, record:=rec<PerfectPointsFormat | >, test:=false)
   record`ConeTestVectors:=Get_y(A: record:=record);
  x:=GetConeVector(A : record:=record);
  xr:=ChangeRing(x,QQ);
  M, mu:=DMinimalVectors(x,A : record:=record);
  MM:=Matrix(M);
  Mperp:=Basis(Nullspace(Transpose(MM)));
  y1:=record`ConeTestVectors[1];
  y1r:=ChangeRing(y1,QQ);
  Ar:=ChangeRing(A,QQ);
  while not (Mperp eq []) do
    print "Starting to compute neighbours";
    size:=#Mperp;
    z:=Vector(Mperp[1]);
    i:=1;
    
    zr:=ChangeRing(z,QQ);
    
    if xAxtr(z,A) le 0 and (zr,y1r) ge 0 then
      z:=-z; // z is not in the closure of the first cone.
    end if;
    
    if method eq 1 then 
      x,M,mu:= VoronoiNeighbour1(A,x,z : Limit:=Limit, record:=record);
    else 
      x,M,mu:=VoronoiNeighbour2(A,x,z : Limit:=Limit, record:=record, test:=test);
    end if;
    xr:=ChangeRing(x,QQ);
    MM:=Matrix(M);
    Mperp:=Basis(Nullspace(Transpose(MM)));
    assert #Mperp lt size;
  end while;
  return x, M;
end function;


// This is a copy of the corresponding function in the 
// normalizer-package

function OrbitReps(Stab, Dirs)
  R:=IndexedSet(Dirs);
  Sx:=GeneratorsSequence(Stab);
  Ps:= [ []^^#Sx ];
  for i in [1..#R] do
    for j in [1..#Sx] do
      X:= R[i] * Sx[j];
      Append(~Ps[j], Index(R, X));
    end for;
  end for;

  Orbs:= Orbits( sub< Sym(#R) | Ps > );
  return R[ [ Rep(x) : x in Orbs ] ];
end function;


intrinsic GetPerfectPoints(A:: Mtrx : Limit:=10^6, max:=Infinity(), method:=2, 
	       PerfectPoints:=rec< PerfectPointsFormat | >, test:=false)
          ->[]
{compute all D-perfect points with D-minimum 1 up to action 
 of Aut(A)^tr/<- id>}
  
  if not (IsHyperbolic(A) and IsIntegralMatrix(A)) then
    error "Argument must be an integral hyperbolic matrix.";
  end if;

  A:=ChangeRing(A,ZZ);

  pp:= assigned PerfectPoints`Points select PerfectPoints else 
                rec< PerfectPointsFormat | Complete:=false>;

  if assigned pp`Points and 
    (assigned pp`Complete and pp`Complete) then 
    return pp`Points;
  end if;
  pp`Watson:=WatsonProcess(A) eq A; 
  pp`ConeTestVectors:=Get_y(A);
  t0:=Realtime();
  x, Mx:=FirstPerfect(A : record:=pp, method:=method, Limit:=Limit, test:=test);
  t1:=Realtime(t0);
  print "Got first perfect point after",t1,"seconds. It has ", #Mx ," D-minimal vectors";
  pp`Points:=[x];
  pp`directions:=[Directions(x,A : M:=Mx, record:=pp)];
  print "Found ",#pp`directions[1]," directions of first point";
  Sx,_,Lx:=PerfectStab(x,A);
  pp`Lat:=[Lx];
  pp`Stabilizers:=[Sx];
  pp`Connectors:=[];
  pp`Neighbours:=[];
  
  // If there are no directions, then there is just one perfect point
  // and we can leave
  
  if #pp`directions[1] eq 0 then
    S,_,L:=PerfectStab(x,A);
    Append(~pp`Stabilizers, S);
    Append(~pp`Lat, L);
    pp`AutGrp:=pp`Stabilizers;
    pp`Complete:=true;
    return pp;
  end if;

  inspected:=0;
  
  while inspected lt #pp`Points and #pp`Points le max do
    inspected +:= 1;
    print "Testing number ",inspected;
    x:=pp`Points[inspected];
    Sx:=pp`Stabilizers[inspected];
    // Append(~pp`Stabilizers, Sx);
    // Append(~pp`Lat, Lx);
    for r in OrbitReps(Sx, pp`directions[inspected]) do
      if method eq 1 then
        y,My:=VoronoiNeighbour1(A,x,r : Limit:=Limit, record:=pp);
      else 
      y,My:=VoronoiNeighbour2(A,x,r : Limit:=Limit, record:=pp, test:=test);
      end if;
      Sy,Trans,Ly:=PerfectStab(y,A);

      for i in [1..#pp`Points] do
        ok,T:=PerfectConnector(y, pp`Points[i],A : Trans:=Trans);
        if ok then
          Append(~pp`Connectors, <inspected, i, y, T>);
          continue r;
        end if;
      end for;
      
      Append(~pp`Points,y);
      Append(~pp`Stabilizers,Sy);
      Append(~pp`Lat, Ly);
      Append(~pp`Neighbours, <inspected, #pp`Points>);
      Append(~pp`directions, Directions(y,A: record:=pp, M:=My));
      print #pp`Points, ": Found point with ",#pp`directions[#pp`directions],"directions";
    end for;
  end while;
  
  pp`Complete:= inspected eq #pp`Points;
  
  if pp`Complete then
    pp`AutGrp:= &cat [GeneratorsSequence(Stab) : 
             Stab in pp`Stabilizers] cat [T[4] : T in pp`Connectors | 
             T[4] ne IdentityMatrix(ZZ, Nrows(A))];
  end if;

  PerfectPoints:=pp;
  
  return PerfectPoints;
end intrinsic;


// Now put everything together and obtain the automorphism 
// group of A (modulo -I_n)

intrinsic AutGroup(A:: Mtrx : 
          PerfectPoints:=rec< PerfectPointsFormat | >,
          Watson:=true,Limit:=10^6, max:=Infinity(),
	  method:=2)
          -> GrpMat
{Calculate the automorphism group of an integral 
 hyperbolic lattice represented by its Gram matrix.}

  if not (IsHyperbolic(A) and IsIntegralMatrix(A)) then
    error "Argument must be an integral hyperbolic matrix.";
  end if;
  
  n:=Nrows(A);
  GLn:=GL(n,ZZ);
  GLnQ:=GL(n,QQ);
  ZZn:=Basis(StandardLattice(n));
  ZZn:=[ChangeRing(Vector(l),QQ) : l in ZZn];

  if (assigned PerfectPoints`Complete and PerfectPoints`Complete) and 
      assigned PerfectPoints`AutGrp then
    Gens:=[ChangeRing(g, ZZ) : g in PerfectPoints`AutGrp];
    Aut:=sub< GLn | Gens cat [-GLn!1]>;
  else
    if Watson then
      B,T:=WatsonProcess(A);
      pp:=GetPerfectPoints(B : method:=method);
      print "Starting stabilizer computation.";
      Gens:=[T^-1*ChangeRing(g,QQ)*T : g in pp`AutGrp];
      AutB:=sub< GLnQ | Gens >;
      Stab:=MyStabGens(AutB,action,ZZn);
      Aut:=sub< GLn | [ChangeRing(g,ZZ) : g in Stab] cat [-GLn!1]>;
      assert forall{g : g in GeneratorsSequence(Aut) | g*A*Transpose(g) eq A};
    else
      pp:=GetPerfectPoints(A);
      Gens:=[ChangeRing(g,ZZ) : g in pp`AutGrp];
      Aut:=sub< GLn | Gens cat [-GLn!1]>;
    end if;
  end if;

  return Aut,#pp`Points;
end intrinsic;


intrinsic DVoronoiGraph(A::Mtrx :
	  Limit:=10^6, max:=Infinity(),
	  method:=2, PerfectPoints:=rec< PerfectPointsFormat | > ) 
          -> [],[],[]
{Compute the residue class graph Gamma_D/Aut(A) (see 
 e.g. [Opg, Def 1.10 and Lemma 2.1]).
 The first output are D-perfect points forming the vertex set of
 Gamma_D/Aut(A), the second gives a list of pairs symbolizing the
 (direct) neighbouring relations between the vertices and the third 
 one gives the connecting family together with the targets.}

  if not (IsHyperbolic(A) and IsIntegralMatrix(A)) then
    error "Argument must be an integral hyperbolic matrix.";
  end if;

  if assigned PerfectPoints`Complete and PerfectPoints`Complete and 
      assigned PerfectPoints`Points  and assigned PerfectPoints`Neighbours 
      and assigned PerfectPoints`Connectors then
  
    Vertices:=PerfectPoints`Points;
    Edges:=[<Vertices[e[1]],Vertices[e[2]]> : 
             e in PerfectPoints`Neighbours];
    Connects:=[<<Vertices[c[1]],Vertices[c[2]]>,c[4]> : 
               c in PerfectPoints`Connectors];
  else
    pp:=GetPerfectPoints(A : PerfectPoints:=PerfectPoints, 
		             Limit:=Limit, max:=max, method:=method);
    Vertices:=pp`Points;
    Edges:=[<Vertices[e[1]],Vertices[e[2]]> : 
             e in pp`Neighbours];
    Connects:=[<<Vertices[c[1]],Vertices[c[2]]>,c[4]> :
                c in pp`Connectors];
  end if;
  
  return Vertices, Edges, Connects;
end intrinsic;
