// Extracted from QuatRefCode.tex by MagmaCode 3.0 on Tue Oct 21 22:39:07 2025
// Documentation in QuatRefCode.pdf
//
// Don Taylor, The University of Sydney
// First version: 11 January 2021

qC := function(m)
  if m gt 2 then
    mm := IsEven(m) select m else 2*m;
    F<z> := CyclotomicField(2*mm);
    K<zz> := sub< F | z + z^-1 >; // the real subfield
    _<t> := PolynomialRing(K);
    C<ii> := ext< K | t^2 + 1 >;
    ok, f := IsIsomorphic(F,C);
    assert ok;
    A<i,j,k> := QuaternionAlgebra< K | -1, -1 >;
    es := IsEven(m) select Eltseq(f(z)^2) else Eltseq(f(z)^4);
    zeta := es[1]+es[2]*i;
  else
    A<i,j,k> := QuaternionAlgebra< Rationals() | -1, -1 >;
    zeta := (m eq 1) select Id(A) else -Id(A);
  end if;
  return MatrixGroup< 1, A | [zeta] >;
end function;

a2C := function(m)
  mm := case< m mod 4 | 0 : m, 2 : 2*m, default : 4*m>;
  F<z> := CyclotomicField(2*mm);
  q := mm div 4;
  a2 := z^-q + z^q;
  K<zz> := sub< F | z + z^-1 >;
  _<t> := PolynomialRing(K);
  C<ii> := ext< K | t^2 + 1 >;
  ok, f := IsIsomorphic(F,C);
  assert ok;
  A<i,j,k> := QuaternionAlgebra< K | -1, -1 >;
  es := case< m mod 4 | 0 : f(z)^2, 2 : f(z)^4, default : f(z)^8>;
  zeta := es[1]+es[2]*i;
  return MatrixGroup< 1, A | [zeta] >, a2;
end function;

a5C := function(m)
  if m mod 5 ne 0 then
    mm := IsEven(m) select 5*m else 10*m;
  else
    mm := IsEven(m) select m else 2*m;
  end if;
  F<z> := CyclotomicField(2*mm);
  q := 2*mm div 5;
  xi := z^q;
  a5 := -2*xi^2 - 2*xi^3 -1;
  K<zz> := sub< F | z + z^-1 >;
  _<t> := PolynomialRing(K);
  C<ii> := ext< K | t^2 + 1 >;
  ok, f := IsIsomorphic(F,C);
  assert ok;
  A<i,j,k> := QuaternionAlgebra< K | -1, -1 >;
  es := f(z)^(2*mm div m);
  zeta := es[1]+es[2]*i;
  return MatrixGroup< 1, A | [zeta] >, A!a5;
end function;

qD := function(m)
  if m gt 2 then
    mm := IsEven(m) select m else 2*m;
    F<z> := CyclotomicField(2*mm);
    K<zz> := sub< F | z + z^-1 >;
    _<t> := PolynomialRing(K);
    C<ii> := ext< K | t^2 + 1 >;
    ok, f := IsIsomorphic(F,C);
    assert ok;
    A<i,j,k> := QuaternionAlgebra< K | -1, -1 >;
    es := IsEven(m) select Eltseq(f(z)) else Eltseq(f(z)^2);
    zeta := es[1]+es[2]*i;
  else
    A<i,j,k> := QuaternionAlgebra< Rationals() | -1, -1 >;
    zeta := (m eq 1) select i^2 else i;
  end if;
  return MatrixGroup< 1, A | [zeta], [j] >;
end function;

a2D := function(m)
  mm := case< m mod 4 | 0 : m, 6 : 2*m, default : 4*m>;
  F<z> := CyclotomicField(2*mm);
  q := mm div 4;
  a2 := z^-q + z^q;
  K<zz> := sub< F | z + z^-1 >;
  _<t> := PolynomialRing(K);
  C<ii> := ext< K | t^2 + 1 >;
  ok, f := IsIsomorphic(F,C);
  assert ok;
  A<i,j,k> := QuaternionAlgebra< K | -1, -1 >;
  es := case< m mod 4 | 0 : f(z), 6 : f(z)^2, default : f(z)^4>;
  zeta := es[1]+es[2]*i;
  return MatrixGroup< 1, A | [zeta], [j] >, a2;
end function;

qValidParametersDC := function(m,ell)
  pp := [];
  if IsOdd(ell) then
    ok, m0 := IsDivisibleBy(m,ell);
    if ok then pp := [1] cat 
      [r : r in [3..m0 by 2] | 
        GCD(m0,r) eq 1 and m0 eq GCD(m0,(r+1) div 2)*GCD(m0,(r-1) div 2)];
    end if;
  else
    ok, n := IsDivisibleBy(2*m,ell);
    if ok then pp := [1] cat [r : r in [2..n div 2] | 
      GCD(n,r) eq 1 and GCD(n div GCD(n,r+1), n div GCD(n,r-1)) eq 1];
    end if;
  end if;
  return pp;
end function;

qArgCheck := function(m,ell,r)
  ok, n := IsDivisibleBy(2*m,ell);
  if ok then
    ok := (r eq 1) or ((1 lt r) and (r le (n div 2)) and (GCD(n,r) eq 1)
      and (GCD(n div GCD(n,r+1), n div GCD(n,r-1)) eq 1));
  end if;
  return ok;
end function;

qDC := function(m,ell,r)
  error if not qArgCheck(m,ell,r), "invalid argument(s)";
  n := 2*m div ell;
  K := qD(m);
  A := BaseRing(K);
  x := (K.1)[1,1]; y := (K.2)[1,1];
  gens := [Matrix(A,2,[0,1, 1,0]),Matrix(A,2,[x^n,0, 0,1]),Matrix(A,2,[x,0, 0,x^r]),
    Matrix(A,2,[y,0, 0,y^-1])];
  G := MatrixGroup< 2,A | gens >;
  assert #G eq 8*m*ell;
  return G;
end function;

a2DC := function(m,ell,r)
  error if not qArgCheck(m,ell,r), "invalid argument(s)";
  n := 2*m div ell;
  K, a2 := a2D(m);
  A := BaseRing(K);
  x := (K.1)[1,1]; y := (K.2)[1,1];
  gens := [Matrix(A,2,[0,1, 1,0]),Matrix(A,2,[x^n,0, 0,1]),Matrix(A,2,[x,0, 0,x^r]),
    Matrix(A,2,[y,0, 0,y^-1])];
  G := MatrixGroup< 2,A | gens >;
  assert #G eq 8*m*ell;
  return G, a2;
end function;

qDD := function(m,ell)
  error if ell le 1, "second argument must be at least 2";
  error if m notin [ell,2*ell], "invalid arguments";
  D<x,y> := qD(m);
  A := BaseRing(D);
  refs := [ Matrix(A,2,[0,1, 1,0]), Matrix(A,2,[y[1,1],0,0,1]) ];
  if m eq ell then
    Append(~refs, Matrix(A,2, [x[1,1],0, 0,1]) );
  else 
    Append(~refs, Matrix(A,2, [x[1,1],0, 0,x[1,1]^-1]) );
  end if;
  G := MatrixGroup< 2,A | refs >;
  assert #G eq 32*m*ell;
  return G;
end function;

qT := function()
  _<t> := PolynomialRing(Rationals());
  K<a2,tau> := NumberField([t^2-2,t^2-t-1]);
  A<i,j,k> := QuaternionAlgebra<K | -1,-1>;
  return MatrixGroup< 1,A | [Matrix(A,1,[x]) : x in [A.1,(1/2)*(-1+i+j+k)]] >;
end function;

qO := function()
  _<t> := PolynomialRing(Rationals());
  K<a2,tau> := NumberField([t^2-2,t^2-t-1]);
  A<i,j,k> := QuaternionAlgebra<K | -1,-1>;
  omega := (1/2)*(-1+i+j+k);
  gamma := a2^-1*(1+i);
  return MatrixGroup< 1,A | [Matrix(A,1,[x]) : x in [omega,gamma]] >;
end function;

qI := function()
  _<t> := PolynomialRing(Rationals());
  K<a2,tau> := NumberField([t^2-2,t^2-t-1]);
  A<i,j,k> := QuaternionAlgebra<K | -1,-1>;
  sigma := (1/2)*(tau^-1 + i + tau*j);
  return MatrixGroup< 1,A | [Matrix(A,1,[x]) : x in [A.1,sigma]] >;
end function;

qTT := function()
  T := qT();
  A := BaseRing(T);
  zi := (T.1)[1,1];
  omega := (T.2)[1,1];
  zj := omega^-1*zi*omega;
  r1 := Matrix(A,2,[zi,0,0,1]);
  r2 := Matrix(A,2,[1,0,0,omega]);
  r3 := Matrix(A,2,[0,zj,zj^-1,0]);
  return MatrixGroup< 2,A | r1, r2, r3 >, {@ x : x in T @};
end function;

qTD2d := function()
  T := qT();
  A<i,j,k> := BaseRing(T);
  K<a2,tau> := BaseRing(A);
  omega := 1/2*(-1+i+j+k);
  delta := 1/a2*(i-j);
  g1 := Matrix(A,2, [0,1, 1,0]);
  g2 := Matrix(A,2, [i,0, 0,1]);
  g3 := Matrix(A,2, [j,0, 0,1]);
  g4 := Matrix(A,2, [i,0, 0,-j]);
  g5 := Matrix(A,2, [omega,0, 0, delta*omega*delta^-1]);
  return MatrixGroup< 2,A | g1,g2,g3,g4,g5 >;
end function;

qTC2d := function()
  T := qT();
  A<i,j,k> := BaseRing(T);
  K<a2,tau> := BaseRing(A);
  delta := 1/a2*(i-j);
  g1 := (T.1)[1,1];
  g2 := (T.2)[1,1];
  r1 := Matrix(A,2, [0,1, 1,0]);
  r2 := Matrix(A,2, [-1,0, 0,1]);
  r3 := Matrix(A,2, [g1,0, 0,delta*g1*delta^-1]);
  r4 := Matrix(A,2, [g2,0, 0,delta*g2*delta^-1]);
  return MatrixGroup< 2,A | r1,r2,r3,r4 >;
end function;

qTC1d := function()
  T := qT();
  A<i,j,k> := BaseRing(T);
  K<a2,tau> := BaseRing(A);
  delta := 1/a2*(i-j);
  g1 := (T.1)[1,1];
  g2 := (T.2)[1,1];
  r1 := Matrix(A,2, [0,1, 1,0]);
  r2 := Matrix(A,2, [g1,0, 0,delta*g1*delta^-1]);
  r3 := Matrix(A,2, [g2,0, 0,delta*g2*delta^-1]);
  return MatrixGroup< 2,A | r1,r2,r3 >;
end function;

qOO := function()
  O := qO();
  A := BaseRing(O);
  r1 := Matrix(A,2,[(O.1)[1,1],0,0,1]);
  r2 := Matrix(A,2,[(O.2)[1,1],0,0,1]);
  r3 := Matrix(A,2,[0,1, 1,0]);
  return MatrixGroup< 2,A | r1, r2, r3 >;
end function;

qOT := function()
  O := qO();
  A := BaseRing(O);
  omega := (O.1)[1,1];
  gamma := (O.2)[1,1];
  r1 := Matrix(A,2,[omega,0, 0,1]);
  r2 := Matrix(A,2,[0,1, 1,0]);
  r3 := Matrix(A,2,[0,omega^-1*gamma*omega, omega^-1*gamma^-1*omega,0]);
  return MatrixGroup< 2,A | r1, r2, r3 >;
end function;

qOD2 := function()
  O := qO();
  A<i,j,k> := BaseRing(O);
  gamma := (O.2)[1,1];
  K<a2,tau> := BaseRing(A);
  r1 := Matrix(A,2,[gamma^2,0, 0,1]);
  r2 := Matrix(A,2,[0,1, 1,0]);
  r3 := Matrix(A,2,[0,x,x^-1,0]) where x is 1/a2*(j+k);
  r4 := Matrix(A,2,[0,x,x^-1,0]) where x is 1/a2*(1+j);
  D2 := sub< O | O.2^2, O.1^-1*O.2^2*O.1 >;
  return MatrixGroup< 2,A | r1,r2,r3,r4 >;
end function;

qOC2 := function()
  O := qO();
  A<i,j,k> := BaseRing(O);
  K<a2,tau> := BaseRing(A);
  g1 := (O.1)[1,1];
  g2 := (O.2)[1,1];
  r1 := Matrix(A,2, [0,1, 1,0]);
  r2 := Matrix(A,2, [-1,0, 0,1]);
  r3 := Matrix(A,2, [g1,0, 0,g1]);
  r4 := Matrix(A,2, [g2,0, 0,g2]);
  return MatrixGroup< 2,A | r1,r2,r3,r4 >;
end function;

qOC1b := function()
  O := qO();
  A<i,j,k> := BaseRing(O);
  K<a2,tau> := BaseRing(A);
  g1 := (O.1)[1,1];
  g2 := (O.2)[1,1];
  beta := hom< O -> O | O.1, -O.2 >; // central automorphism
  r1 := Matrix(A,2, [0,1, 1,0]);
  r2 := Matrix(A,2, [g1,0, 0,g1]);
  r3 := Matrix(A,2, [g2,0, 0,-g2]);
  return MatrixGroup< 2,A | r1,r2,r3 >, beta;
end function;

qOC1d := function()
  O := qO();
  A<i,j,k> := BaseRing(O);
  K<a2,tau> := BaseRing(A);
  delta := 1/a2*(i-j);
  g1 := (O.1)[1,1];
  g2 := (O.2)[1,1];
  r1 := Matrix(A,2, [0,1, 1,0]);
  r2 := Matrix(A,2, [g1,0, 0,delta*g1*delta^-1]);
  r3 := Matrix(A,2, [g2,0, 0,delta*g2*delta^-1]);
  return MatrixGroup< 2,A | r1,r2,r3 >;
end function;

qII := function()
  I := qI();
  A<i,j,k> := BaseRing(I);
  r1 := Matrix(A,2,[0,1,1,0]);
  r2 := Matrix(A,2,[1,0,0,i]);
  r3 := Matrix(A,2,[1,0,0,(I.2)[1,1]]);
  return MatrixGroup< 2,A | r1, r2, r3 >;
end function;

qIC2 := function()
  I := qI();
  A<i,j,k> := BaseRing(I);
  K<a2,tau> := BaseRing(A);
  g1 := (I.1)[1,1];
  g2 := (I.2)[1,1];
  r1 := Matrix(A,2, [0,1, 1,0]);
  r2 := Matrix(A,2, [-1,0, 0,1]);
  r3 := Matrix(A,2, [g1,0, 0,g1]);
  r4 := Matrix(A,2, [g2,0, 0,g2]);
  return MatrixGroup< 2,A | r1,r2,r3,r4 >;
end function;

qIC2t := function()
  I := qI();
  A<i,j,k> := BaseRing(I);
  K<a2,tau> := BaseRing(A);
  M := Matrix(A,1,[-1/2*tau - 1/2*i + 1/2*(tau - 1)*k]);
  theta := hom< I -> I | -I.1, M >;
  g1 := (I.1)[1,1];
  g2 := (I.2)[1,1];
  r1 := Matrix(A,2, [0,1, 1,0]);
  r2 := Matrix(A,2, [-1,0, 0,1]);
  r3 := Matrix(A,2, [g1,0, 0,-g1]);
  r4 := Matrix(A,2, [g2,0, 0,M[1,1]]);
  return MatrixGroup< 2,A | r1,r2,r3,r4 >, theta;
end function;

qIC1t := function()
  I := qI();
  A<i,j,k> := BaseRing(I);
  K<a2,tau> := BaseRing(A);
  delta := 1/a2*(i-j);
  M := Matrix(A,1,[-1/2*tau - 1/2*i + 1/2*(tau - 1)*k]);
  theta := hom< I -> I | -I.1, M >;
  g1 := (I.1)[1,1];
  g2 := (I.2)[1,1];
  r1 := Matrix(A,2, [0,1, 1,0]);
  r2 := Matrix(A,2, [g1,0, 0,-g1]);
  r3 := Matrix(A,2, [g2,0, 0,M[1,1]]);
  return MatrixGroup< 2,A | r1,r2,r3 >, theta;
end function;

qIC1j := function() 
  I := qI();
  A<i,j,k> := BaseRing(I);
  K<a2,tau> := BaseRing(A);
  delta := 1/a2*(i-j);
  g1 := (I.1)[1,1];
  g2 := (I.2)[1,1];
  r1 := Matrix(A,2, [0,1, 1,0]);
  r2 := Matrix(A,2, [g1,0, 0,j*g1*j^-1]);
  r3 := Matrix(A,2, [g2,0, 0,j*g2*j^-1]);
  return MatrixGroup< 2,A | r1,r2,r3 >;
end function;

qmuT := function(d)
  error if d mod 6 ne 0, "d must be divisible by 6";
  mu := qC(d);
  A<i,j,k> := BaseRing(mu);
  g0 := Matrix(A,2, [z,0,0,z]) where z is (mu.1)[1,1];
  g1 := Matrix(A,2, [i,0, 0,-i]);
  g2 := 1/2*Matrix(A,2, [-1+i, 1+i, -1+i, -1-i]);
  s := Matrix(A,2, [0,-j, j,0] );
  return MatrixGroup< 2, A | g0, g1, g2, s >;
end function;

qmuO := function(d)
  error if d mod 4 ne 0, "d must be divisible by 4";
  mu, a2 := a2C(d);
  A<i,j,k> := BaseRing(mu);
  g0 := Matrix(A,2, [z,0,0,z]) where z is (mu.1)[1,1];
  g1 := 1/2 * Matrix(A, 2, [-1+i, 1+i, -1+i, -1-i]);
  g2 := Matrix(A, 2, [kappa, 0, 0, kappa^-1]) where kappa is 1/a2*(1+i);
  s := Matrix(A,2, [0,-j, j,0] );
  return MatrixGroup< 2, A | g0, g1, g2, s >;
end function;

qmu2O := function(d)
  error if d mod 4 ne 0, "d must be a multiple of 4";
  error if d mod 16 eq 0, "d must not be a multiple of 16";
  mu, a2 := a2C(d);
  A<i,j,k> := BaseRing(mu);
  g0 := Matrix(A,2, [z,0, 0,z]) where z is (mu.1)[1,1];
  g1 := 1/2 * Matrix(A, 2, [-1+i, 1+i, -1+i, -1-i]);
  g2 := Matrix(A, 2, [kappa, 0, 0, kappa^-1]) where kappa is 1/a2*(1+i);
  s := Matrix(A,2, [0,-j, j,0] );
  return MatrixGroup< 2, A | g0^2, g1, g2^2, g0*g2, s >;
end function;

qmuI := function(d)
  error if (d mod 4) ne 0 and (d mod 6) ne 0 and (d mod 10) ne 0,
    "d must be divisible by 4, 6, or 10";
  mu, a5 := a5C(d);
  A<i,j,k> := BaseRing(mu);
  tau := 1/2*(1+a5);
  g0 := Matrix(A,2, [z,0, 0,z]) where z is (mu.1)[1,1];
  g1 := 1/2*Matrix(A,2, [tau^-1 + tau*i, 1, -1, tau^-1 - tau*i]);
  g2 := Matrix(A,2, [i,0, 0,-i]);
  s := Matrix(A,2, [0,-j, j,0] );
  return MatrixGroup< 2, A | g0, g1, g2, s >;
end function;

qDotProduct := function(u,v)
  error if v notin Generic(Parent(u)), "vectors must be in the same space";
  vbar := Vector([Conjugate(x) : x in Eltseq(v)]);
  return InnerProduct(u,vbar);
end function;

qReflection := function(a,xi)
  V := Generic(Parent(a));
  A := BaseRing(V);
  lambda := (qDotProduct(a,a))^-1*(One(A) - xi);
  return Matrix([ v - qDotProduct(v,a)*lambda*a : v in Basis(V)]);
end function;

reflections := func< G | 
  &join[Class(G,c[3]) : c in Classes(G) | Rank(Matrix(c[3])-G.0) eq 1 ] >;

qImprimitive := function(K,H,n)
  error if not DerivedGroup(K) subset H, "K/H must be abelian";
  A := BaseRing(K);
  gens := [ PermutationMatrix(A,g) : g in Generators(Sym(n)) ];
  M := RSpace(A,n);
  gens cat:= [ qReflection(M.1,xi[1,1]) : xi in Generators(H) ];
  gens cat:= [ qReflection(M.1 - xi[1,1]*M.2, -Id(A)) : xi in Generators(K) ];
  return MatrixGroup< n, A | gens >;
end function;

qqO1 := function()
  _<t> := PolynomialRing(Rationals());
  K<a,b> := NumberField([t^2-2,t^2-3]);
  A<ii,jj,kk> := QuaternionAlgebra< K | -1, -1 >;
  gamma := 1/2*(-1 + b*ii);
  delta := 1/2*(-1 + b*jj);
  V := RSpace(A,2);
  v1 := elt< V | 0, gamma - 1 >;
  v2 := elt< V | a*jj, gamma >;
  v3 := elt< V | -x*a^-1*delta*jj, x*gamma^2 > where x is (1-kk);
  rts := [v1, v2, v3];
  refs := [ qReflection(v,gamma) : v in rts ];
  return MatrixGroup< 2,A | refs >, rts;
end function;

qqO2 := function()
  _<t> := PolynomialRing(Rationals());
  K<a,b> := NumberField([t^2-2,t^2-3]);
  A<ii,jj,kk> := QuaternionAlgebra< K | -1, -1 >;
  gamma := 1/2*(-1 + b*ii);
  delta := 1/2*(-1 + b*jj);
  V := RSpace(A,2);
  v1 := elt< V | 0, gamma - 1 >;
  v2 := elt< V | a*jj, gamma >;
  v3 := elt< V | -x*a^-1*delta*jj, x*gamma^2 > where x is (1-kk);
  v4 := elt< V | a*b*(1-kk), 0 >;
  rts := [v1, v2, v3, v4];
  refs := [ qReflection(v,gamma) : v in rts ];
  return MatrixGroup< 2,A | refs >;
end function;

qqO3 := function()
  _<t> := PolynomialRing(Rationals());
  K<a,b,c> := NumberField([t^2-2,t^2-3,t^2-5]);
  A<ii,jj,kk> := QuaternionAlgebra< K | -1, -1 >;
  gamma := 1/2*(-1 + b*ii);
  tau := 1/2*(2 + b + ii + b*jj - kk);
  V := RSpace(A,2);
  v1 := elt< V | a*jj, 1 + b >;
  v2 := elt< V | gamma*a*jj, 1 + b >;
  v3 := elt< V | a*jj, gamma*(1 + b) >;
  v4 := elt< V | tau*gamma*(1+kk)/a, tau*gamma >;
  v5 := elt< V | x*gamma^2*(1-kk)/a, -x*gamma > where x is ii*tau;
  rts := [v1,v2,v3,v4,v5];
  refs := [ qReflection(v,-1) : v in rts ];
  return MatrixGroup< 2,A | refs >;
end function;

qqP1 := function()
  A<ii,jj,kk> := QuaternionAlgebra<Rationals() | -1,-1>;
  V := RSpace(A,2);
  v1 := elt< V | 1,0>;
  v2 := elt< V | x,x > where x is 1/2*(1-ii);
  rts := [v1,v2];
  refs := [ qReflection(v,jj) : v in rts ];
  return MatrixGroup< 2,A | refs >, rts;
end function;

qqP2 := function()
  A<ii,jj,kk> := QuaternionAlgebra<Rationals() | -1,-1>;
  V := RSpace(A,2);
  v1 := elt< V | 1,0>;
  v2 := elt< V | 0,1>;
  v3 := elt< V | x,x*jj > where x is 1/2*(1-ii);
  r1 := qReflection(v1, kk);
  r2 := qReflection(v2, jj);
  r3 := qReflection(v3, kk);
  return MatrixGroup< 2,A | [r1, r2, r3 ] >, [v1,v2,v3];
end function;

qqP3 := function()
  _<t> := PolynomialRing(Rationals());
  K<a> := NumberField(t^2-2);
  A<ii,jj,kk> := QuaternionAlgebra< K | -1, -1 >;
  V := RSpace(A,2);
  v1 := elt< V | 1+a,ii>;
  v2 := elt< V | 1+a,jj>;
  v3 := elt< V | 1+a,kk>;
  v4 := elt< V | x/a*(1+jj),x*jj> where x is 1 + a^-1*(1+jj);
  v5 := elt< V | x/a*(1-kk),x*kk> where x is 1 + a^-1*(1+kk);
  rts := [v1,v2,v3,v4,v5];
  refs := [ qReflection(v,-1) : v in rts ];
  return MatrixGroup< 2,A | refs >;
end function;

qqQ := function()
  _<t> := PolynomialRing(Rationals());
  K<c> := NumberField(t^2-5);
  A<ii,jj,kk> := QuaternionAlgebra< K | -1, -1 >;
  V := RSpace(A,3);
  alpha := 2^-1*(1-ii-jj-c*kk);
  v1 := elt< V | 2,0,0 >;
  v2 := elt< V | 0,2,0 >;
  v3 := elt< V | ii,ii,alpha>;
  v4 := elt< V | alpha*ii,1,ii>;
  rts := [v1,v2,v3,v4];
  refs := [ qReflection(v,-1) : v in rts ];
  return MatrixGroup< 3,A | refs >;
end function;

qqR := function()
  _<t> := PolynomialRing(Rationals());
  K<tau> := NumberField(t^2-t-1);
  A<ii,jj,kk> := QuaternionAlgebra< K | -1, -1 >;
  V := RSpace(A,3);
  c := 2*tau - 1;
  sigma := 1 - tau;
  alpha := 2^-1*(1-ii-jj-c*kk);
  omega := 1/2*(-1+ii+jj+kk);
  v1 := elt< V | 2,0,0 >;
  v2 := elt< V | 0,2,0 >;
  v3 := elt< V | ii,ii,alpha>;
  v4 := elt< V | tau*omega^2, sigma*omega*jj,1>;
  rts := [v1,v2,v3,v4];
  refs := [ qReflection(v,-1) : v in rts ];
  return MatrixGroup< 3,A | refs >;
end function;

qqS1 := function()
  A<ii,jj,kk> := QuaternionAlgebra<Rationals() | -1,-1>;
  V := RSpace(A,4);
  v1 := elt< V | x,x,0,0> where x is 1 - ii;
  v2 := elt< V | x,0,x,0> where x is 1 - ii;
  v3 := elt< V | x,0,0,x> where x is 1 - ii;
  v4 := elt< V | 1,ii,jj,kk >;
  rts := [v1,v2,v3,v4];
  refs := [ qReflection(v,-1) : v in rts ];
  return MatrixGroup< 4,A | refs >;
end function;

qqS2 := function()
  A<ii,jj,kk> := QuaternionAlgebra<Rationals() | -1,-1>;
  V := RSpace(A,4);
  v1 := elt< V | x,x,0,0> where x is 1 - ii;
  v2 := elt< V | x,0,0,x> where x is 1 - ii;
  v3 := elt< V | 1,ii,jj,kk >;
  v4 := elt< V | 0,0,0,2>;
  rts := [v1,v2,v3,v4];
  refs := [ qReflection(v,-1) : v in rts ];
  return MatrixGroup< 4,A | refs >;
end function;

qqS3 := function()
  A<ii,jj,kk> := QuaternionAlgebra<Rationals() | -1,-1>;
  V := RSpace(A,4);
  v1 := elt< V | x,x*ii,0,0> where x is 1 - ii;
  v2 := elt< V | x,x*jj,0,0> where x is 1 - ii;
  v3 := elt< V | x,0,x*kk,0> where x is 1 - ii;
  v4 := elt< V | x*ii,0,x,0> where x is 1 - ii;
  v5 := elt< V | 1,jj,kk,ii >;
  rts := [v1,v2,v3,v4,v5];
  refs := [ qReflection(v,-1) : v in rts ];
  return MatrixGroup< 4,A | refs >;
end function;

qqT := function()
  _<t> := PolynomialRing(Rationals());
  K<tau> := NumberField(t^2-t-1);
  A<ii,jj,kk> := QuaternionAlgebra< K | -1, -1 >;
  V := RSpace(A,4);
  sigma := 1 - tau;
  omega := 1/2*(-1+ii+jj+kk);
  p := 1/2*(tau + sigma*ii - jj);
  v1 := elt< V | 1/2*tau,1/2*sigma,-1/2,0 >;
  v2 := elt< V | 1/2*sigma,1/2*tau,0,-1/2 >;
  v3 := elt< V| 1,0,0,0 >;
  v4 := elt< V| x,x*p^-1*ii*p,x*p^-1*jj*p,x*p^-1*kk*p> where x is 1/2;
  rts := [v1,v2,v3,v4];
  refs := [ qReflection(v,-1) : v in rts ];
  return MatrixGroup< 4,A | refs >;
end function;

qqU := function()
  _<t> := PolynomialRing(Rationals());
  K<tau> := NumberField(t^2-t-1);
  A<ii,jj,kk> := QuaternionAlgebra< K | -1, -1 >;
  V := RSpace(A,5);
  omega := 1/2*(-1+ii+jj+kk);
  v1 := elt< V| 2,0,0,0,0>;
  v2 := elt< V| 0,0,0,2,0>;
  v3 := elt<V | 0,1,omega,omega,1>;
  v4 := elt<V | 1,0,1,omega,omega>;
  v5 := elt< V | 0,1,ii*omega*ii,jj*omega*jj,ii >;
  rts := [v1,v2,v3,v4,v5];
  refs := [ qReflection(v,-1) : v in rts ];
  return MatrixGroup< 5,A | refs >; 
end function;

