# File: QuatAlg.gap
# Description: GAP functions for computing fundamental
# regions for the action of the unit group of a 
# quaternion algebra over the Rationals.
# Author: Assaf Wool
# email: assafwool@yahoo.com
# Date: July 2007

# Define the global quaternion algebra record
QuatGlob:=rec();

# Function to analyze D and initialize the quaternion object
findOrderType:=function(D)
  local pr,i,len,p,q,r,ok;

  if D<1 then 
    Print("ERROR: discriminant ",D," must be positive\n");
    return 0;    
  fi;

  pr:=FactorsInt(D);
# check for multiple primes
  len:=Length(pr);
  i:=1;
  while i<len do
    if pr[i]=pr[i+1] then
      Print("ERROR: discriminant divisible by prime power ",pr[i],"\n");
      return 0;
    fi;
    i:=i+1;
  od;

# check for bad number of primes 
  if RemInt(len,2)<>0 or len=0 then
    Print("ERROR: discriminant has ",len," divisors\n");
    return 0;    
  fi; 

# find type, u, v
  QuatGlob:=rec();
  QuatGlob.disc:=D;

# find volume, 2 and 3 torsion, genus and number of generators
  QuatGlob.vol:=1/6;
  QuatGlob.n2:=1;
  QuatGlob.n3:=1;
  for p in pr do
    QuatGlob.vol:=QuatGlob.vol*(p-1);
    if p>2 then 
      QuatGlob.n2:=QuatGlob.n2*(1-Legendre(-1,p));
    fi;
    if p>3 then 
      QuatGlob.n3:=QuatGlob.n3*(1-Legendre(-3,p));
    fi;
    if p=2 then 
      QuatGlob.n3:=QuatGlob.n3*2;
    fi;
  od;
  QuatGlob.twog:=2-QuatGlob.n2/2-QuatGlob.n3*2/3+
    QuatGlob.vol;
  QuatGlob.numgen:=QuatGlob.twog+QuatGlob.n2+QuatGlob.n3;
  if QuatGlob.numgen>QuatGlob.twog then
    QuatGlob.numgen:=QuatGlob.numgen-1;
  fi;

  ok:=0;
  p:=pr[1];
  q:=pr[2];
  if p=2 and len=2 then
    if Legendre(-1,q)=-1 then
      QuatGlob.u:=1;
      QuatGlob.v:=q;
      QuatGlob.type:=1;
      QuatGlob.den:=2;
      ok:=1;
    elif Legendre(-2,q)=-1 then
      QuatGlob.u:=2;
      QuatGlob.v:=q;
      QuatGlob.type:=2;
      QuatGlob.den:=2;
      ok:=1;
    fi;
  fi;
  if p<>2 and len=2 then
    if Legendre(-1,p)=-1 and Legendre(-1,q)=-1 then
      QuatGlob.u:=1;
      QuatGlob.v:=D;
      QuatGlob.type:=2;
      QuatGlob.den:=2;
      ok:=1;
    elif Legendre(-p,q)=-1 and Legendre(q,p)=-1 then
      QuatGlob.u:=p;
      QuatGlob.v:=q;
      QuatGlob.type:=2;
      QuatGlob.den:=2;
      ok:=1;
    elif Legendre(-q,p)=-1 and Legendre(p,q)=-1 then
      QuatGlob.u:=q;
      QuatGlob.v:=p;
      QuatGlob.type:=2;
      QuatGlob.den:=2;
      ok:=1;
    fi;
  fi;
  if ok=1 then
    QuatGlob.alg:=QuaternionAlgebra(Rationals,
      -QuatGlob.u,QuatGlob.v);
    QuatGlob.basis:=Basis(QuatGlob.alg);
    return 1;
  fi;

# general types 
  r:=5;
  while true do 
    if Legendre(-1,r)=1 then
# this question takes care of both odd and even cases
      if ForAll(pr,x->Legendre(x,r)=-1) then
        QuatGlob.u:=D;
        QuatGlob.v:=r;
        QuatGlob.type:=3;
        QuatGlob.den:=2*r;
        QuatGlob.aux:=RootMod(-D,r);
        QuatGlob.alg:=QuaternionAlgebra(Rationals,
          -QuatGlob.u,QuatGlob.v);
        QuatGlob.basis:=Basis(QuatGlob.alg);
        return 1;
      fi;
    fi;
    r:=NextPrimeInt(r);
  od;
  Print("ERROR: bad type, shouldn't be here\n");
  return 0; 
end;

# Function that finds all pairs (x,y) with norm of type x^2+Uy^2
# between min and max, appends them to the global lists and
# sorts by norm. 
findComplexNorms:=function(U, max, min)
  local x,y,n,order,tmpl;

# initialize u 
# QuatGlob.u:=U;
# Initialize if starting from 0
  if min=0 then
    QuatGlob.x:=[];
    QuatGlob.y:=[];
    QuatGlob.norm:=[];
  fi;

# Find the norms
  x:=0;
  while x*x<=max do
    y:=0;
    n:=x*x+U*y*y;
    while n<=max do
      if n>=min and (QuatGlob.type<>3 or RemInt(n,QuatGlob.v)=0) then
        Append(QuatGlob.x,[x]);
        Append(QuatGlob.y,[y]);
        Append(QuatGlob.norm,[n]);
      fi;
      y:=y+1;
      n:=x*x+U*y*y;
    od;
    x:=x+1;
  od; 

# Sort the lists according to norm
  order:=Sortex(QuatGlob.norm);
  tmpl:=ShallowCopy(QuatGlob.x);
  QuatGlob.x:=Permuted(tmpl,order);
  tmpl:=ShallowCopy(QuatGlob.y);
  QuatGlob.y:=Permuted(tmpl,order);

end;

# Function that finds all quaternion elements in the 
# algebra B(-U,V), in a maximal order defined by type.
# It uses the (x,y) pairs computed and sorted by norm.
findQuatElem:=function(V,type)
  local pos,i,maxNm,tmpNm,a,b,c,d,tmpL,size,j,t;

# initialize v
# QuatGlob.v:=V;
# Initialize the element array
  QuatGlob.elem:=[];

# make sure type is 1-3
  if type<>1 and type<>2 and type<>3 then
    Print("Error, type must be 1 or 2 or 3\n");
    return;
  fi;

  maxNm:=QuatGlob.norm[Length(QuatGlob.norm)];
  i:=1;
  while QuatGlob.norm[i]*V<maxNm do
    tmpNm:=QuatGlob.norm[i]*V+QuatGlob.den*QuatGlob.den;
    pos:=PositionSorted(QuatGlob.norm,tmpNm);
    while pos<=Length(QuatGlob.norm) and QuatGlob.norm[pos]=tmpNm do
      a:=QuatGlob.x[pos];
      b:=QuatGlob.y[pos];
      c:=QuatGlob.x[i];
      d:=QuatGlob.y[i];

      tmpL:=[[a,b,c,d]];
      size:=1;

# expand all plus-minus signs necessary
      if (a<>0 and b<>0) then
        j:=1;
        while j<=size do
          Add(tmpL, [tmpL[j][1],-tmpL[j][2],tmpL[j][3],tmpL[j][4]]);
          j:=j+1;
        od;
        size:=size*2;
      fi;      
      if (c<>0) then
        j:=1;
        while j<=size do
          Add(tmpL, [tmpL[j][1],tmpL[j][2],-tmpL[j][3],tmpL[j][4]]);
          j:=j+1;
        od;
        size:=size*2;
      fi;      
      if (d<>0) then
        j:=1;
        while j<=size do
          Add(tmpL, [tmpL[j][1],tmpL[j][2],tmpL[j][3],-tmpL[j][4]]);
          j:=j+1;
        od;
        size:=size*2;
      fi;

# filter elements according to order type
      if (type=1 and RemInt(a,2)=RemInt(b,2) and 
          RemInt(a,2)=RemInt(c,2) and 
          RemInt(a,2)=RemInt(d,2)) or 
         (type=2 and RemInt(a,2)=RemInt(c,2) and 
          RemInt(b,2)=RemInt(d,2)) then 
        Append(QuatGlob.elem,tmpL/QuatGlob.den);
      fi;
      if type=3 then
        for t in tmpL do
          if RemInt(t[2],QuatGlob.v)=0 and
             RemInt(t[2]-t[4],2)=0  and
             RemInt(t[3]-QuatGlob.aux*(t[4]-t[2]),
                    QuatGlob.v)=0 and
             RemInt(t[1]-t[3]-QuatGlob.aux*(t[2]-t[4]),
                    2*QuatGlob.v)=0 then
            Add(QuatGlob.elem,t/QuatGlob.den);
          fi;
        od;
      fi;

      pos:=pos+1;
    od;
    i:=i+1;
  od;
end;

# function that finds a list of non-redundant 
# generators, using the element list. It finds
# all elements whose centers are not in any
# isometric circle of a previous element 
findQuatGenerators:=function()
  local i,len,n,Nm,Nm2,Nmg,first,Last,c,d,g,h,flag;

# initialize generators list
  QuatGlob.gen:=[];

# compute all norms
  len:=Length(QuatGlob.elem);
  Nm:=[1..len];
  Apply(Nm, i->QuatGlob.elem[i][3]^2+QuatGlob.u*QuatGlob.elem[i][4]^2);
# find first non zero norm, and find last index for each norm
  first:=PositionNonZero(Nm);
  Nm2:=0;
  n:=0;
  Last:=[];
  for i in [first..len] do
    if Nm[i]>Nm2 then
      n:=i-1;
      Nm2:=Nm[i];
    fi;
    Last[i]:=n;
  od;

# go over elements and add to generators if
# no other element contains the center in its
# isometric circle.
  n:=first;
  while n<=len do
    h:=QuatGlob.elem[n];

# first go over the current generators, they have the best
# chance of reducing the norm, this should improve
# speed
    flag:=1;
    for g in QuatGlob.gen do
      Nmg:=g[3]*g[3]+QuatGlob.u*g[4]*g[4];
      if Nmg<Nm[n] then
        c:=g[1]*h[3]-g[3]*h[1]+QuatGlob.u*(g[2]*h[4]-g[4]*h[2]);
        d:=g[1]*h[4]-g[4]*h[1]-g[2]*h[3]+g[3]*h[2];
        Nm2:=c*c+QuatGlob.u*d*d;
        if Nm2<Nm[n] then
          flag:=0;
          break;
        fi;
      fi;
    od;

    if flag=1 then
      i:=first;
      while i<=Last[n] do
        g:=QuatGlob.elem[i];
        c:=g[1]*h[3]-g[3]*h[1]+QuatGlob.u*(g[2]*h[4]-g[4]*h[2]);
        d:=g[1]*h[4]-g[4]*h[1]-g[2]*h[3]+g[3]*h[2];
        Nm2:=c*c+QuatGlob.u*d*d;
        if Nm2<Nm[n] then
          flag:=0;
          break;
        fi;
        i:=i+1;
      od;
      if flag=1 then
        Add(QuatGlob.gen, h);
      fi;
    fi;
    n:=n+1;
  od;
# if no torsion then make sure generator list is symmetric
# by adding symmetric generators and unique.
  if QuatGlob.n2=0 and QuatGlob.n3=0 then
    len:=Length(QuatGlob.gen);
    for i in [1..len] do
      g:=QuatGlob.gen[i];
      Add(QuatGlob.gen,[g[1],g[2],-g[3],-g[4]]);
    od;
    QuatGlob.gen:=Unique(QuatGlob.gen);
  fi;
end;

# Function that compares directions for sorting 
# along the circle
cmpDir:=function(d1,d2)
  local v1,v2;

# case of different y signs
  if d1[3]>=0 and d2[3]<0 then
    return true;
  fi;
  if d1[3]<0 and d2[3]>=0 then
    return false;
  fi;
# case of y=0
  if d1[3]=0 then 
    return (d1[2]>0 and (d2[3]<>0 or d2[2]<0));
  fi;
  if d2[3]=0 then
    return (d2[2]<0 and (d1[3]<>0 or d1[2]>0));
  fi;
# case of same y signs, no y=0
  v1:=d1[2]/d1[3];
  v2:=d2[2]/d2[3];
  if v1>v2 then
    return true;
  fi;
  return false;
end;

# Function that finds inverse pairs of generators
# and arranges them around the circle. Special
# care must be given for the case of U=1, where i
# is a special generator.
arrangeGenerators:=function()
  local i,j,len,len2,g,dir,dir2,x,y,inv,map,gen,last,ng;

  dir:=[];
  i:=1;
  len:=Length(QuatGlob.gen);
  while i<=len do
    g:=QuatGlob.gen[i];
    x:=g[1]*g[3]-QuatGlob.u*g[2]*g[4];
    y:=g[1]*g[4]+g[2]*g[3];
    Add(dir,[i,x,y]);  
    i:=i+1;  
  od;
# sort along the circle
  Sort(dir, cmpDir);

# Must do special treatment for u=1, remove
# generators with the same direction, replace
# generators by equivalent ones with inverse
# in the right half circle
  if QuatGlob.u=1 then
    dir2:=[];
    if len>0 then
      Add(dir2,dir[1]);
    fi;
    for i in [2..len] do
      if cmpDir(dir[i-1],dir[i]) then
        Add(dir2,dir[i]);
      else
        QuatGlob.gen[dir[i][1]]:=[1,1,0,0];
      fi;
    od;
    dir:=dir2;

    last:=-[1,dir[1][2],dir[1][3]];
    i:=1;
    g:=QuatGlob.gen[dir[i][1]];
    while (not cmpDir(last,[1,dir[i][2],dir[i][3]])) do
      x:=-g[1]*g[3]-g[2]*g[4];
      y:=-g[1]*g[4]+g[2]*g[3];
# if inverse is in opposite half, multiply by i on the right
# and remember that u=1.
      if ((not cmpDir([1,x,y],last)) or cmpDir([1,x,y],-last)) then
        ng:=[-g[2],g[1],g[4],-g[3]];
        if ng[1]<0 or (ng[1]=0 and ng[2]<0) then
          QuatGlob.gen[dir[i][1]]:=-ng;
        else
          QuatGlob.gen[dir[i][1]]:=ng;        
        fi;
      fi;
# special case of first/last element of order 2
      g:=QuatGlob.gen[dir[i][1]];
      if ((i=1 and g[1]=0) or 
          (not cmpDir([1,dir[i][2],dir[1][3]],last) and g[1]=0)) then
        ng:=[-g[2],g[1],g[4],-g[3]];
        if ng[1]<0 or (ng[1]=0 and ng[2]<0) then
          QuatGlob.gen[dir[i][1]]:=-ng;
        else
          QuatGlob.gen[dir[i][1]]:=ng;        
        fi;
      fi;
      i:=i+1;  
      g:=QuatGlob.gen[dir[i][1]];
    od;
    while i<=Length(dir) do
      QuatGlob.gen[dir[i][1]]:=[1,1,0,0];
      i:=i+1;
    od;
  fi;

# find inverses
  inv:=[];
  for i in [1..len] do
    inv[i]:=0;
  od;
  for i in [1..len] do
    if inv[i]=0 then
      if QuatGlob.gen[i][1]=0 then
        inv[i]:=i;
      else
        for j in [i+1..len] do
          if (QuatGlob.gen[i][1] = QuatGlob.gen[j][1] and
              QuatGlob.gen[i][2] = -QuatGlob.gen[j][2] and
              QuatGlob.gen[i][3] = -QuatGlob.gen[j][3] and
              QuatGlob.gen[i][4] = -QuatGlob.gen[j][4]) then
            inv[i]:=j;
            inv[j]:=i;
            break;
          fi;
        od;
      fi;
    fi;
  od;

# create map for removing elements with no inverse
# and putting them in their order on the circle
  len2:=Length(dir);
  map:=[];
  j:=1;
# make room for generator i if u=1
  if QuatGlob.u=1 then
    j:=2;
  fi;
  for i in [1..len2] do 
    if inv[dir[i][1]]<>0 then
      map[dir[i][1]]:=j;
      j:=j+1;
    fi;
  od;

# Put sorted generators and inverse information in the 
# global object
  QuatGlob.imap:=[];
  gen:=[];
# put i in front if u=1
  if QuatGlob.u=1 then
    Add(QuatGlob.imap,1);
    Add(gen,LinearCombination(QuatGlob.basis,[0,1,0,0]));
  fi;
  for i in [1..len2] do
    if inv[dir[i][1]]<>0 then
      Add(QuatGlob.imap,map[inv[dir[i][1]]]);
      Add(gen,LinearCombination(QuatGlob.basis,QuatGlob.gen[dir[i][1]]));
    fi;
  od;
  QuatGlob.gen:=gen;
end;

# Function for finding the relators
findRelations:=function()
  local i,len,done,rel,g,curr,coef,ok,tor2,tor3,ng;

  QuatGlob.rel:=[];
  QuatGlob.grel:=[];

  done:=[];
  len:=Length(QuatGlob.gen);
# nothing to do if there are no generators
  if len<1 then
    return 0;
  fi;

  for i in [1..len] do
    done[i]:=0;
  od;
  
  ng:=0;
  for i in [1..len] do
    if QuatGlob.imap[i]<>i then
      ng:=ng+1/2;
    else
      ng:=ng+1;
    fi;
  od;
  tor2:=0;
  tor3:=0;
  ok:=1;
# find relations by using the inverse and circle order
  for i in [1..len] do
    if done[i]=0 then
      rel:=[i];
      g:=QuatGlob.gen[i];
# special care of order 2 generators, add relation
      coef:=Coefficients(QuatGlob.basis,g);
      if coef[1]=0 then
        Add(QuatGlob.rel,[i]);
        Add(QuatGlob.grel,g);
	tor2:=tor2+1;
      fi;
      done[i]:=1;
      curr:=QuatGlob.imap[i]-1;
      if curr<1 then
        curr:=len;
      fi;
      while done[curr]=0 do
        Add(rel,curr);
        g:=g*QuatGlob.gen[curr];
        done[curr]:=1;
        curr:=QuatGlob.imap[curr]-1;
        if curr<1 then
          curr:=len;
        fi;
      od;
# add relation to global object, and check if it is a real
# relation. 
      Add(QuatGlob.rel,rel);
      Add(QuatGlob.grel,g);
      coef:=Coefficients(QuatGlob.basis,g);
      if coef[1]=0 then
        tor2:=tor2+1;
      elif coef[1]=1/2 or coef[1]=-1/2 then
        tor3:=tor3+1;
      elif coef=[1,0,0,0] or coef=[-1,0,0,0] then
        ng:=ng-1;
      else
        ok:=0;
      fi;
    fi;
  od;
  QuatGlob.tor2:=tor2;
  QuatGlob.tor3:=tor3;
  if tor2=0 and tor3=0 then
    ng:=ng+1;
  fi;
  QuatGlob.ng:=ng;
  if ok=1 and 
    [tor2,tor3,ng]<>[QuatGlob.n2,QuatGlob.n3,QuatGlob.numgen] then
    ok:=0;
    Print("WARN: Thought to be impossible (bug maybe?)\n",
      [tor2,tor3,ng]," ",[QuatGlob.n2,QuatGlob.n3,QuatGlob.numgen],
      "\n");
  fi;
  return ok;
end;

# function for finding a torsion element using w=i
# in case the group has no torsion. 
findTorsionW:=function()
  local w,i,done,len,g,flag,curr,coef,numl;

  w:=QuatGlob.basis[2];

  done:=[];
  len:=Length(QuatGlob.gen);

  for i in [1..len] do
    done[i]:=0;
  od;

  i:=1;
  flag:=true;
  while flag and i<=len/2 do
    if done[i]=0 then

      g:=QuatGlob.gen[i];
      done[i]:=1;
      curr:=QuatGlob.imap[i];
      if curr>len/2 then
        g:=g*w;
        curr:=curr-len/2;
      fi;
      curr:=curr-1;
      if curr=0 then
        g:=g*w;
        curr:=curr+len/2;
      fi;
      while done[curr]=0 do

        g:=g*QuatGlob.gen[curr];
        done[curr]:=1;
        curr:=QuatGlob.imap[curr];
        if curr>len/2 then
          g:=g*w;
          curr:=curr-len/2;
        fi;
        curr:=curr-1;
        if curr=0 then
          g:=g*w;
          curr:=curr+len/2;
        fi;
      od;

      coef:=Coefficients(QuatGlob.basis,g);
      if coef[1]=0 then
        numl:=[NumeratorRat(coef[2]),NumeratorRat(coef[3]),NumeratorRat(coef[4])];
        QuatGlob.wtor:=g/Gcd(numl);
        return;
      fi;
    fi;
    i:=i+1;
  od;
  Print("WARN: can't find w-torsion point\n");
end;

# function for reducing redundant generators
reduceGens:=function()
  local i,len,inRel,j,coef,max,ind,norm,gc,
        flag,inv,rlen,tmpl,rel,newRel,map,
        tgen,tinv,trel,tgrel;

# go over relations, find where each generator is
# ignore order 2 generators relation
#  Print(QuatGlob.rel,"\n",QuatGlob.imap,"\n");
  len:=Length(QuatGlob.rel);
  inRel:=[];
  for i in [1..len] do
    coef:=Coefficients(QuatGlob.basis,QuatGlob.grel[i]);
    if not(Length(QuatGlob.rel[i])=1 and coef[1]=0) then
      for j in QuatGlob.rel[i] do
        inRel[j]:=i;
      od;
    fi;
  od;

  flag:=false;
  while (not flag) do
# find relevant generator with  
# maximal c^2+u*d^2 norm. Ignore generators
# whose inverse is in the same relation
    max:=0;
    ind:=0;
    for i in [1..len] do
      coef:=Coefficients(QuatGlob.basis,QuatGlob.grel[i]);
      if (coef=[1,0,0,0] or coef=[-1,0,0,0]) then
        for j in QuatGlob.rel[i] do
          if (QuatGlob.imap[j]=j or 
              inRel[QuatGlob.imap[j]]<>i) then
            gc:=Coefficients(QuatGlob.basis,QuatGlob.gen[j]);
            norm:=gc[3]*gc[3]+QuatGlob.u*gc[4]*gc[4];
            if norm>max then
              max:=norm;
              ind:=j;
            fi;
          fi;
        od;
      fi;
    od;
    if max>0 then
# update relations by substitution of gen[ind]
#      Print("removing: ",ind, " ", QuatGlob.gen[ind]);
      rel:=QuatGlob.rel[inRel[ind]];
      i:=Position(rel,ind);
      rlen:=Length(rel);
      tmpl:=Concatenation(rel{[i+1..rlen]},rel{[1..i-1]});
# find the location of the inverse gen, and substitute
      inv:=QuatGlob.imap[ind];
      rel:=QuatGlob.rel[inRel[inv]];
      i:=Position(rel,inv);
      rlen:=Length(rel);
      newRel:=Concatenation(rel{[1..i-1]},tmpl,rel{[i+1..rlen]});
# update inRel, and relations in QuatGlob
      QuatGlob.rel[inRel[ind]]:=[];
      QuatGlob.rel[inRel[inv]]:=newRel;
#      Print(" ", newRel);
      inRel[ind]:=0;
      for i in tmpl do
        inRel[i]:=inRel[inv];
      od;
      inRel[inv]:=0;
    else
      flag:=true;
    fi;
  od;
# find removed generators
  j:=1;
  map:=[];
  len:=Length(QuatGlob.gen);
  for i in [1..len] do
    if inRel[i]<>0 then
      Add(map,j);
      j:=j+1;
    else
      Add(map,0);
    fi;
  od;
# update generators and inverse
  tgen:=[];
  tinv:=[];
  for i in [1..len] do
    if inRel[i]<>0 then
      Add(tgen,QuatGlob.gen[i]);
      Add(tinv,map[QuatGlob.imap[i]]);      
    fi;
  od;
  QuatGlob.gen:=tgen;
  QuatGlob.imap:=tinv;
# update relations 
  trel:=[];
  tgrel:=[];
  i:=1;
  len:=Length(QuatGlob.rel);
  while i<=len do
    if Length(QuatGlob.rel[i])>0 then 
       Add(tgrel,QuatGlob.grel[i]);
       Apply(QuatGlob.rel[i],j->map[j]);
       Add(trel,QuatGlob.rel[i]);
    fi;
    i:=i+1;
  od;
  QuatGlob.grel:=tgrel;
  QuatGlob.rel:=trel;
end;

# find and print all vertices
findVertices:=function()
  local coef,len,i,j,g,cg,w,x,y,first,co,
        dir,tc,ind,pos,rel,rlen,norm,allDir,
        ok,num,den;

  len:=Length(QuatGlob.rel);
  norm:=1;
  allDir:=[];
  for i in [1..len] do
    g:=QuatGlob.grel[i];
    rel:=QuatGlob.rel[i];
    coef:=Coefficients(QuatGlob.basis,g);
    if coef=[1,0,0,0] or coef=[-1,0,0,0] then
      if len > 1 then
        Print("ERROR: more than one surface relation\n");
        return;
      fi;
      if not IsBound(QuatGlob.wtor) then
        Print("ERROR: no w-torsion point found\n");
        return; 
      fi;
# init g and find the starting point for the relation orbit
      g:=QuatGlob.wtor;
      coef:=Coefficients(QuatGlob.basis,g);
      norm:=coef[1]^2+QuatGlob.u*coef[2]^2-
            QuatGlob.v*(coef[3]^2+QuatGlob.u*coef[4]^2);
# be careful about the direction of the w-torsion
# point, sign of the real coefficient has to be considered
      if (coef[2]>0 and coef[2]^2*QuatGlob.u>norm) or
         (coef[2]<0 and coef[2]^2*QuatGlob.u<norm) then
        dir:=[1,-QuatGlob.u*coef[4],coef[3]];
      else
        dir:=[1,QuatGlob.u*coef[4],-coef[3]];
      fi;
      ind:=1;
      for j in [1..Length(QuatGlob.gen)] do
        tc:=Coefficients(QuatGlob.basis,QuatGlob.gen[j]);
        x:=tc[1]*tc[3]-QuatGlob.u*tc[2]*tc[4];
        y:=tc[1]*tc[4]+tc[2]*tc[3];
#        if cmpDir(dir,[1,x,y]) then
        if (not cmpDir([1,x,y],dir)) then
          ind:=j;
          break;
        fi;
      od;
      pos:=Position(QuatGlob.rel[1],QuatGlob.imap[ind]);
      rlen:=Length(QuatGlob.rel[1]);
      rel:=Concatenation(QuatGlob.rel[1]{[pos+1..rlen]},
           QuatGlob.rel[1]{[1..pos]});
    fi;
    if coef[1]=0 then
      for j in rel do
        cg:=Coefficients(QuatGlob.basis,g);
        if cg=[0,1,0,0] or cg=[0,-1,0,0] then
          Print(j,": 0 (2-tor)\n");
        else
          if (cg[2]>0 and cg[2]^2*QuatGlob.u>norm) or
             (cg[2]<0 and cg[2]^2*QuatGlob.u<norm) then
            allDir[j]:=[1,-QuatGlob.u*cg[4],cg[3]];
          else
            allDir[j]:=[1,QuatGlob.u*cg[4],-cg[3]];
          fi;
# extract constants before printing
          num:=GcdInt(NumeratorRat(cg[4]*QuatGlob.u),NumeratorRat(cg[3]));
          den:=LcmInt(DenominatorRat(cg[4]*QuatGlob.u),DenominatorRat(cg[3]));
          Print(j,": ",num/den/(cg[3]*cg[3]+QuatGlob.u*cg[4]*cg[4])," x ");
          Print("(",cg[2],"S(",QuatGlob.u,")");
          if cg[2]>0 then
            Print("-S(",norm,"))");
          else
            Print("+S(",norm,"))");
          fi;
          Print("/", "S(",QuatGlob.u*QuatGlob.v,") x ");
          Print("(",-QuatGlob.u*cg[4]/num*den);
          if cg[3]>=0 then
            Print("+");
          fi;
          Print(cg[3]/num*den,"S(-",QuatGlob.u,"))\n");
# print the relation coefficients too
	  Print(j,": rel " , cg[1], " ", cg[2], " ", cg[3], " ", cg[4],"\n");
        fi;
# shift relation cyclically
        g:=QuatGlob.gen[QuatGlob.imap[j]]*g*QuatGlob.gen[j];
      od;
    elif coef[1]=1/2 or coef[1]=-1/2 then
      for j in rel do
        cg:=Coefficients(QuatGlob.basis,g);
        if (cg[2]>0 and 4*cg[2]^2*QuatGlob.u>3) or
           (cg[2]<0 and 4*cg[2]^2*QuatGlob.u<3) then
          allDir[j]:=[1,-QuatGlob.u*cg[4],cg[3]];
        else
          allDir[j]:=[1,QuatGlob.u*cg[4],-cg[3]];
        fi;
        Print(j,": (",2*cg[2],"S(",QuatGlob.u,")");
        if cg[2]>0 then
          Print("-S(3))");
        else
          Print("+S(3))");
        fi;
        Print("/(",2*(cg[3]*cg[3]+QuatGlob.u*cg[4]*cg[4]),
              "S(",QuatGlob.u*QuatGlob.v,") x ");
        Print("(",-QuatGlob.u*cg[4]);
        if cg[3]>=0 then
          Print("+");
        fi;
        Print(cg[3],"S(-",QuatGlob.u,"))\n");
# print the relation coefficients too
	  Print(j,": rel " , cg[1], " ", cg[2], " ", cg[3], " ", cg[4],"\n");
# shift relation cyclically
        g:=QuatGlob.gen[QuatGlob.imap[j]]*g*QuatGlob.gen[j];
      od;
    else
      Print("ERROR: non-torsion relation ",rel,"\n");
      return;
    fi;
  od;
# verify that vertices are ordered as they should
  ok:=true;
  len:=Length(QuatGlob.gen);
  for i in [2..len-2] do
    if (not cmpDir(allDir[i],allDir[i+1])) then
      ok:=false;
      break;
    fi;
  od;
  if not ( (cmpDir(allDir[1],allDir[2]) and 
            cmpDir(allDir[len],allDir[1]) ) or
           (cmpDir(allDir[1],allDir[2]) and 
            cmpDir(allDir[len-1],allDir[len]) ) or 
           (cmpDir(allDir[len],allDir[1]) and
            cmpDir(allDir[len-1],allDir[len]) ) ) then
    ok:=false;
  fi;
  if (not ok) then
    Print("ERROR: vertices are not sorted along the circle\n");
  fi;
end;

# Runner function, take D and get a full object with 
# generators and fundamental region.
createQuatObj:=function(D)
  local ok,max,min;

  ok:=findOrderType(D);
  if ok=0 then 
    return 0;
  fi;
 
  max:=QuatGlob.disc*QuatGlob.disc;
  min:=0;
  ok:=0;
  while ok=0 do
    findComplexNorms(QuatGlob.u, max, min);
    findQuatElem(QuatGlob.v,QuatGlob.type);
    findQuatGenerators();
    arrangeGenerators();
    ok:=findRelations();
    if ok=0 then
      min:=max+1;
      max:=2*max;
      Print("not enough generators, increase max norm to ",max,"\n");
    fi;
  od;
# use w=i to find torsion element in group with w
# if group has no torsion.
  if QuatGlob.n2=0 and QuatGlob.n3=0 then
    findTorsionW();
  fi;
# find redundant generators until all relations are eliminated
# except torsion or surface relations
  reduceGens();
# find and print vertices of region
  findVertices();
  return 1;
end;