// freeze;

/***********************************************************************************
 *                                                                                 *
 * M. Kirschmer, M. Ruether,                                                       *
 * The explicit membership problem for discrete two-generator subgroups of SL(2,R) *
 *                                                                                 *
 * Some complimentary Magma code, M. Kirschmer & M. Ruether, July 2016             *
 *                                                                                 *
 **********************************************************************************/

declare type SL2R2Gen;

declare attributes SL2R2Gen:
  A, B,		// user given generators
  U, V,		// normalized witness pair
  WordsUV,	// how to express \bar{U},\bar{V} in {\bar{A},\bar{B}}
  MinusOne,	// 0 or a word w such that w(A,B)=-1.
  S,		// how to conjugate U, V into "nice generators".
  K,		// base field for the original generators
  plc,		// embedding K -> R which yields G -> SL2(R)
  PSL,		// should we interpret G as a subgroup of PSL2(R)?
  finite,	// data for finite groups
  type;		
/* 
type can be one of the following
"nd" -> not discrete,
<0,0,0> -> free group & discrete
<a,0,0> -> Z/aZ * Z free product & discrete
<a,b,0> -> Z/aZ * Z/bZ free product & discrete
<a,b,c> -> Triangle group <a,b,c> & discrete
"finite" -> finite (and thus discrete)
"ed" -> elementary & discrete but not finite
"comm" -> commutator of finite order & discrete
*/

oo:= Infinity();
ZZ:= Integers();
MAXTRIES:= 100;


myNF:= func< plc | Type(plc) eq Infty select Rationals() else NumberField(plc) >;
function myEval(x, plc : ForceReal:= false, Precision:= 300)
  if Type(plc) eq Infty then
    return ForceReal select RealField(Precision) ! x else x;
  end if;
  return Evaluate(x, plc : Precision:= Precision);
end function;
comm:= func< g,h | g*h*g^-1*h^-1 >;
mySign:= func< x, plc | Sign(myEval(x, plc)) >;
l2norm:= func< x | x[1]^2 + x[2]^2 >;
// Test if P is in the circle centered at (Q1, 0) with radius r^(1/2)
InCircle:= func< P, Q1, r, plc | myEval(r - ((P[1] - Q1)^2 + P[2]^2) , plc) gt 0 >;
OnCircle:= func< P, Q1, r | r eq 0 select P[1] eq Q1 else (P[1] - Q1)^2 + P[2]^2 eq r >;
// Cosh of the hyperbolic distance
HypDistance:= func< P, Q | 1 + ((P[1]-Q[1])^2 + (P[2]-Q[2])^2)/(2*P[2]*Q[2]) >;

function mysqrt(x, plc)
  K:= myNF(plc);
  p:= Polynomial([K | -x,0, 1]);
  r:= Roots(p);
  if #r ne 0 then return myEval(r[1,1], plc) ge 0 select r[1,1] else r[2,1], plc; end if;
  L:= ext< K | p >;
  D:= [ d[1] : d in Decomposition(L, plc) ];
  return L.1, Evaluate(L.1, D[1]) ge 0 select D[1] else D[2];
end function;

apply:= function(X, x)
  if x cmpeq oo then
    return X[2,1] eq 0 select oo else X[1,1]/X[2,1];
  end if;
  d:= X[2,1] * x + X[2,2];
  return d eq 0 select oo else (X[1,1] * x + X[1,2])/d;
end function;

pi:= Pi(RealField(300));
function GetExp(z, A, t, plc);
  k:= 0;
  if t ge 2 then	// t is the projective order.
    u2:= -A[1,2]/A[2,1];
    uR:= Sqrt(myEval(u2, plc : ForceReal));
    s:= u2 eq 1 select 1 else -1;
    if z[1] eq 0 then
      b:= (u2 eq 1 and myEval(z[2], plc) lt 1) or (u2 ne 1 and myEval(z[2], plc) gt uR) select pi else 0;
    else
      M:= myEval((z[1]^2 + z[2]^2 - u2)/(2*z[1]), plc: ForceReal);
      b:= -Arctan(uR/M);
      if s*mySign(z[1], plc)*Sign(M) eq -1 then b +:= pi; end if;
    end if;
    kk:= -b*t / (2*pi);
    k:= Round(kk) mod t;
    if k ne 0 then z:= apply(A^k, z); end if;
    if Abs(kk - Floor(kk) - 1/2) le 10^-6 then	 // close to 1/2, so we might sit on the boundary!
      B:= A^(s*mySign(z[1], plc));
      if OnCircle(z, -B[2,2]/B[2,1], B[2,1]^-2) then return 0,0; end if;
    end if;
  elif t eq 1 and A[2,1] eq 0 then	// A is a translation 
    k:= Round(myEval(-z[1]/A[1,2], plc));
    z +:= k*A[1,2];
    if 2*z[1] eq A[1,2] or 2*z[1] eq -A[1,2] then return 0,0; end if;
  else
    for i in [1,-1] do
      if i eq -1 then A:= A^-1; end if;
      m:= -A[2,2]/A[2,1]; 
      r2:= A[2,1]^-2;
      
      while true do
        l:= (z[1]-m)^2 + z[2]^2;
	if l eq r2 then return 0,0; end if;
	if myEval(l-r2, plc) gt 0 then break; end if;
	k +:= i;
	z:= apply(A, z);
      end while;
      if k ne 0 then break; end if;
    end for;
  end if;
  return k, z;
end function;

procedure swap(~X, ~Y)
  temp:= X; X:= Y; Y:= temp;
end procedure;

function RationalReconstruction(x,y)
  if y lt x then swap(~x, ~y); end if;
  d:= y-x;
  assert d gt 0;
  l:= Max(0, Ceiling(Log(10, d))) + 2;
  r:= Round( (x+y)/2 * 10^l ) / 10^l;
  assert x lt r and r lt y;
  return r;
end function;

cross:= function(a,b,c,d)
try
  if a cmpeq oo then return (b-d)/(b-c); end if;
  if b cmpeq oo then return (a-c)/(a-d); end if;
  if c cmpeq oo then return (d-b)/(d-a); end if;
  if d cmpeq oo then return (c-a)/(c-b); end if;
  return (a-c)*(b-d)/(b-c)/(a-d);
catch e
  return oo;
end try;
end function;

fixed_exact:= function(X, plc)
  a,b,c,d:= Explode(Eltseq(X));
  disc:= (a+d)^2-4;
  if disc eq 0 then
    return (c eq 0 select [oo] else [(a-d) / (2*c)]), plc;
  elif myEval(myNF(plc) ! disc, plc) lt 0 then 
    D, plc:= mysqrt(-disc, plc);
    if myEval(c, plc) lt 0 then D:= -D; end if;
    return [ (a-d)/(2*c), D/(2*c) ], plc;
  end if;
  if c eq 0 then return [* -b/(a-d), oo *], plc; end if;
  D, plc:= mysqrt(disc, plc); c:= 2*c;
  x:= a-d;
  return Sort([(x-D)/c, (x+D)/c], func< r,s | myEval(r-s, plc) >), plc;
end function;

// Diagonalize an hyperbolic element. The larger eigenvalue comes first.
function diagonalize(A, plc)
  if IsDiagonal(A) then
    S:= Parent(A) ! 1;
  elif A[2,1] eq 0 then
    S:= Matrix(2, [(A[1,1] - A[2,2]) / A[1,2], 1, 0, 1]);
    A[1,2]:= 0;
  else
    f, plc:= fixed_exact(A, plc);
    assert #f eq 2;
//    S:= Matrix(2, [1, -f[1], 1, -f[2]]);
    S:= Matrix(2, [1, -f[1], -1, f[2]]);  // Det > 0
    A:= S*A*S^-1;
    assert IsDiagonal(A);
  end if;
  if myEval(A[1,1] - A[2,2], plc) ge 0 then return A, S, plc; end if;
  return DiagonalMatrix([A[2,2], A[1,1]]), Matrix(2, [BaseRing(A) | 0, 1, -1, 0]) * S, plc;
end function;

eqn:= function(z,w)
  if z cmpeq oo then
    return w cmpeq oo select [0, 0, 1, 0] else  [1, 0, -w, 0];
  end if;
  return w cmpeq oo select [0,0,z,1] else [z, 1, -w*z, -w];
end function;


function MyTransform(F1, F2)
  T:= Type(F1);
  if T eq Infty or ISA(T, RngElt) then
    F1:= [* F1 *];
  end if;
  T:= Type(F2);
  if T eq Infty or ISA(T, RngElt) then
    F2:= [* F2 *];
  end if;

  assert #F1 eq #F2 and #F1 le 3;
  while #F1 lt 3 do
    ok:= exists(i){i: i in [-2..2] | Index(F1, i) eq 0 and Index(F2, i) eq 0 };
    assert ok;
    Append(~F1, i);
    Append(~F2, i);
  end while;

  M:= Matrix([ eqn(F1[i], F2[i]) : i in [1..3] ]);
  K:= Kernel(Transpose(M));
  assert Dimension(K) eq 1;
  return Matrix(2, Eltseq(K.1));
end function;

procedure MarkElementary(G)
  // case 1: G is finite
  if IsFinite(G) then
    G`type:= "finite";
    return;
  end if;
  
  A:= G`A; B:= G`B; plc:= G`plc;
  // case 2: all elements in G are parabolic.
  T:= Trace(comm(A,B));
  if {Trace(A), Trace(B)} subset {-2,2} and T eq 2 then 
    f:= fixed_exact(A, plc)[1]; 		// does not change plc!
    assert f eq apply(B,f);
    T:= MyTransform(f, oo);
    A:= T*A*T^-1;
    B:= T*B*T^-1;
    if {A[1,1], B[1,1]} subset {-1,1} then 
      G`U:= A; G`V:= B;
      G`WordsUV:= [ F.1, F.2 ] where F:= FreeGroup(2);
      G`type:= "ed";
      G`S:= T;
    else
      G`type:= "nd";
    end if;
    return;
  end if;

  // case 3a: G contains only hyperbolic elements.
  if Abs(myEval(Trace(A), plc)) gt 2 and Abs(myEval(Trace(B), plc)) gt 2 and T eq 2 then
    fA, plc:= fixed_exact(A, plc);
    T:= MyTransform(fA, [*0, oo*]); TI:= T^-1;
    A:= T*A*TI; B:= T*B*TI;
    a:= Abs(myEval(A[1,1], plc));
    b:= Abs(myEval(B[1,1], plc));
    l:= Log(a, b);
    L:= BestApproximation(l, 10^20);
    n:= Numerator(L);
    m:= Denominator(L);
    if a^n eq b^m then
      G`U:= A; G`V:= B;
      G`type:= "ed";
      G`WordsUV:= [ F.1, F.2 ] where F:= FreeGroup(2);
      G`S:= T;
      G`plc:= plc;
    else
      G`type:= "nd";
    end if;
    return;
  end if;

  // case 3b: extension of Z by a subgroup of C4.
  // {A,B,AB} contains elliptic and hyperbolic elements -> find them.
  hyp:= false;
  ell:= false;
  for i in [<1,0>, <0,1>, <1,1>] do
     X:= A^i[1] * B^i[2];
     TX:= Trace(X);
     if IsOne(X^4) then 
       if ell cmpeq false then ell:= i; end if;
     //elif TX in {-2,2} or Abs(myEval(TX, plc)) lt 2 then 
     //  G`type:= "nd"; return;
     //else
     elif Abs(myEval(TX, plc)) gt 2 then
       if hyp cmpeq false then hyp:= i; end if;
     end if;
  end for;
  if hyp cmpne false and ell cmpne false then
    U:= A^hyp[1] * B^hyp[2];
    V:= A^ell[1] * B^ell[2];
    f, plc:= fixed_exact(U, plc);
    T:= MyTransform(f, [*0, oo*]); TI:= T^-1;
    U:= T*U*TI;
    V:= T*V*TI;
    if IsOne(V^2) or {V[1,1], V[2,2]} eq {0} then
      G`U:= U; G`V:= V;
      G`type:= "ed";
      G`S:= T;
      G`WordsUV:= [ F.1^hyp[1] * F.2^hyp[2], F.1^ell[1]*F.2^ell[2] ] where F:= FreeGroup(2);
      G`plc:= plc;
    end if;
  end if;

  // the remaining case
  if T eq 2 then 
    G`type:= "nd";
  end if;
  // we let all other cases slip through. They will be discovered later.
end procedure;

// Returns the order of g in PSL.
// 0 means parabolic / hyperbolic, -1 means elliptic of infinite order.
// Otherwise it returns the order n of g in PSL.
// The 2nd return value gives m such that Tr(g) = cos(m*pi/n)
function MyOrder(g, plc)
  T:= myEval(Trace(g), plc);
  if Abs(T) ge 2 then return 0,0,0; end if;
  if not HasFiniteOrder(g) then return -1, 0 ,0; end if;
  o:= Order(g); n:= o;
  t:= myEval(Trace(g)/2, plc : ForceReal);
  m:= Round(Arccos(t)*n/Pi(Parent(t)));
  if IsEven(n) and IsEven(m) then n div:= 2; m div:= 2; end if;
  assert GCD(n,m) eq 1 and m ge 1 and m lt n;
  return n, m, o;
end function;



// The user interface
intrinsic Print(G::SL2R2Gen)
{internal}
  printf "Two generator subgroup of %oSL(2,R) generated by\n%o,\n%o", G`PSL select "P" else "", G`A, G`B;
end intrinsic;

intrinsic IsCoercible(G::SL2R2Gen, A::. ) -> BoolElt, .
{internal}
  return false, "invalid coercion";
end intrinsic;

function create(A,B,plc,PSL)
  T:= Type(plc);
  if T ne Infty and not (ISA(T, PlcNumElt) and IsReal(plc)) then
    return false, "The third argument must be a real infinite place";
  elif Determinant(A) ne 1 or Determinant(B) ne 1 then
    return false, "The matrices must have determinant 1.";
  end if;
  K:= myNF(plc);
  okA, A:= CanChangeRing(A, K);
  okB, B:= CanChangeRing(B, K);
  if not (okA and okB) then return false, "The base rings of A and B must be the base field of the place."; end if;
  G:= New(SL2R2Gen);
  G`A:= A; G`B:= B; G`K:= K; G`plc:= plc; G`PSL:= PSL;
  return true, G;
end function;

intrinsic TwoGeneratorSubgroupSL2R(A::Mtrx, B::Mtrx, plc::.) -> SL2R2Gen
{The subgroup of SL(2,R) generated by A, B using the embedding plc}
  ok, G:= create(A,B,plc,false);
  require ok : G;
  return G;
end intrinsic;

intrinsic TwoGeneratorSubgroupPSL2R(A::Mtrx, B::Mtrx, plc::.) -> SL2R2Gen
{The subgroup of SL(2,R) generated by A, B using the embedding plc}
  ok, G:= create(A,B,plc,true);
  require ok : G;
  return G;
end intrinsic;

intrinsic '.'(G::SL2R2Gen, i::RngIntElt) -> Mtrx
{The i-th generator of the group G}
  require i in [-2..2] : "Index out of range";
  return (Abs(i) eq 2 select G`B else G`A)^Sign(i);
end intrinsic;

// first return value: -1 -> ell (infinite order), 0 -> para/hyp., > 0 -> ell.
normalize:= function(A, plc)
  n, m, o:= MyOrder(A, plc);
  if n eq -1 then return -1, _,_,_;
  elif n eq 0 then 
    k:= 1;
  else
    if m gt n/2 then m:= n-m; end if;
    k:= InverseMod(m, n);
    A:= A^k;
    if IsOdd(o) then o:= 0; end if;
  end if;
  if myEval(Trace(A), plc) lt 0 then A:= -A; end if;
  return n, A, k, o div 2;
end function;

// set up G`type and U, V if necessary
procedure prepare(G)
  if assigned G`type then return; end if;
  MarkElementary(G);			// test is G is elementary and discrete
  if assigned G`type then return; end if;

  C:= comm(G`A,G`B); T:= Trace(C);			// T \ne 2
  plc:= G`plc;
  F:= SLPGroup(2);
  Minus:= 0;

  n, U, k, minus:= normalize(G`A, plc);
  if n eq -1 then G`type:= "nd"; return; end if;
  if minus ne 0 then Minus:= F.1^minus; end if;
  WU:= F.1^k;
  n, V, k, minus:= normalize(G`B, plc);
  if n eq -1 then G`type:= "nd"; return; end if;
  if (minus ne 0) and (Minus cmpeq 0) then Minus:= F.2^minus; end if;
  WV:= F.2^k;

  if myEval(T, plc) le -2 or T eq -2 then 
    G`type:= <0,0,0>; 				// G is free
  elif myEval(T, plc) lt 2 then
    G`type:= "nd";
    n, m, o:= MyOrder( -C, plc );
    if n ne 0 then
      if IsEven(m*n) then Minus:= comm(F.1, F.2)^n; end if;
      if m in {1,2} then G`type:= "comm";
      elif m in {3,4,6} or (m eq 4 and n eq 7) then		// find (R,S)
	TU:= Trace(U); TV:= Trace(V);
        repeat
          if myEval(TV-TU, plc) lt 0 then swap(~U, ~V); swap(~WU, ~WV); end if;
          S:= [ U*V, U^-1*V ];
          T:= [ Trace(s): s in S ];
          E:= [ myEval(t, plc): t in T ];
	  mm, i:= Min(E);
	  done:= mm ge myEval(TV, plc);
	  if not done then V:= U^(i eq 1 select 1 else -1) * V; TV:= T[i]; end if;
        until done;
	if i eq 2 then V:= V^-1; WV:= WV^-1; end if;
	// now we find some funny words in R,S
	R:= U; S:= V; WR:= WU; WS:= WV;
	C:= comm(R,S); WC:= comm(WR, WS);
	s:= InverseMod(m,n);
        if m eq 6 then
	  if #{ Trace(R), Trace(S), Trace(R*S) } eq 1 then 
	    G`type:= <2,3,n>;
	    V:= C^-(2*s)*R; U:= V^-2 * S * V^2;
	    WV:= WC^-(2*s)*WR; WU:= WV^-2 * WS * WV^2;
	  end if;
        elif m eq 4 then
	  if Trace(R*S) eq Trace(R)^2/2 then
	    G`type:= <2,4,n>;
	    V:= (-1)^(1+s mod 2) * C^-s * R; U:= V^2 * R^-1;
	    WV:= WC^-s * WR; WU:= WV^2 * WR^-1;
	  elif n eq 7 and #{Trace(R)+1, Trace(S), Trace(R*S)} eq 1 then
	    G`type:= <2,3,7>;
	    T:= (R*S)^-1; WT:= (WR*WS)^-1;
            U:= -comm(S^-1,R) * comm(S^-1*R*S, R) * comm(R, R*S^-1) * comm(R*S^-1, T) * comm(S^-1*R, T) * comm(S^-1, T) * T;
            WU:= comm(WS^-1,WR) * comm(WS^-1*WR*WS, WR) * comm(WR, WR*WS^-1) * comm(WR*WS^-1, WT) * comm(WS^-1*WR, WT) * comm(WS^-1, WT) * WT;
	    V:= -( comm(S^-1,R)^3 * R*S * comm(T,S^-1*R) * comm(T, R*S^-1) * comm(R*S^-1, R)^2 * comm(R,S^-1)^2 )^2;
	    WV:= ( comm(WS^-1,WR)^3 * WR*WS * comm(WT,WS^-1*WR) * comm(WT, WR*WS^-1) * comm(WR*WS^-1, WR)^2 * comm(WR,WS^-1)^2 )^2;
	  end if;
	elif m eq 3 then
	  if #{ Trace(R), Trace(S), Trace(R*S) } eq 1 then
	    G`type:= <3,3,n>;
	    if IsOdd(n) then		// check if -1 in G. 
	      h:= hom< F -> GL(2, G`K) | G`A, G`B >;
              RR:= h(WR); SS:= h(WS);
	      if IsOne(-(C^s * SS)^3) then
	        Minus:= (WC^s * WS)^3;
              elif IsOne(- (C^-s * RR)^3) then
	        Minus:= (WC^-s * WR)^3;
	      end if;
	    end if;

	    U:= -(C^s*S)^-1; V:= -C^-s * R;
	    WU:= -(WC^s*WS)^-1; WV:= -WC^-s * WR;
	  end if;
        end if;
      end if;
    end if;

  else		// Trace([U,V])> 2:
    while true do
      if myEval(Trace(V) - Trace(U), plc) lt 0 then swap(~U, ~V); swap(~WU, ~WV); end if;
      n1, UV, k1, minus1:= normalize(U*V, plc);
      n2, IV, k2, minus2:= normalize(U^-1*V, plc);
      if n1 eq -1 or n2 eq -1 then G`type:= "nd"; return; end if;

      T:= [ myEval(Trace(UV), plc), myEval(Trace(IV), plc) ];
      m,i:= Minimum(T);
      if k1 eq 1 and k2 eq 1 and m ge myEval(Trace(V), plc) then break; end if;
      if k1 ne 1 then i:= 1; elif k2 ne 1 then i:= 2; end if;		// reducing tr[U,V] takes precedence.
      V:= i eq 1 select UV else IV;
      k:= i eq 1 select k1 else k2;
      n:= i eq 1 select n1 else n2;
      W:= WU^(i eq 1 select 1 else -1) * WV;
      if IsEven(k) and Minus cmpeq 0 and IsOne(-Evaluate(W, [ GL(2,G`K) | G`A, G`B ])^n) then
	Minus:= W^n; 
      end if;
      WV:= W^k;
    end while;
    if myEval(Trace(U*V), plc) gt 0 then V:= V^-1; WV:= WV^-1; end if;
    assert myEval(Trace(U*V), plc) lt 0;
    G`type:= < MyOrder(U, plc), MyOrder(V, plc), MyOrder(-U*V, plc) >;
  end if;
  G`U:= U; G`V:= V; G`WordsUV:= [WU, WV];
  if Minus cmpne 0 then 
    assert IsOne(-Evaluate(Minus, [GL(2, G`K) | G`A, G`B]));
    G`MinusOne:= Minus;
  end if;
end procedure;


/* Predicates */

intrinsic IsDiscrete(G::SL2R2Gen) -> BoolElt
{Returns true if and only if G is discrete}
  prepare(G);
  return G`type cmpne "nd";
end intrinsic;

intrinsic IsFinite(G::SL2R2Gen) -> BoolElt, RngIntElt
{Returns true if and only if G is discrete and free}
  if not assigned G`finite then
    GG:= sub< GL(2, G`K) | G`A, G`B >;
    if IsFinite(GG) then
      G`finite:= <true, GG>;
    else
      G`finite:= <false>;
    end if;
  end if;
  if G`finite[1] then return true, #(G`finite[2]); end if;
  return false, _;
end intrinsic;

intrinsic IsTriangleGroup(G::SL2R2Gen) -> BoolElt, Tup
{Returns true if and only if G is a triangle group}
  prepare(G);
  if Type(G`type) ne Tup or G`type[3] eq 0 then return false, _; end if;
  return true, G`type;
end intrinsic;

intrinsic IsDiscreteAndElementary(G::SL2R2Gen) -> BoolElt
{Returns true if and only if G is elementary and discrete}
  prepare(G);
  return G`type cmpeq "ed" or G`type cmpeq "finite";
end intrinsic;

intrinsic IsDiscreteAndFree(G::SL2R2Gen) -> BoolElt
{Returns true if and only if G is discrete and free of rank 2}
  prepare(G);
  return G`type cmpeq <0,0,0>;
end intrinsic;

intrinsic IsDiscreteAndFreeProduct(G::SL2R2Gen) -> BoolElt, RngIntElt, RngIntElt
{Returns true if and only if G is discrete and a free product}
  prepare(G);
  if Type(G`type) ne Tup then return false, _, _; end if;
  return true, G`type[1], G`type[2];
end intrinsic;



/* Membership */

// Centre and radius^2 of the geodesic through z and w in H.
// radius 0 indicates  vertical line.
function Geodesic(z,w)
  if z[1] eq w[1] then return z[1], 0; end if;
  c:= (w[2]^2 - z[2]^2 + w[1]^2 - z[1]^2) / (2*(w[1] - z[1]));
  r2:= (c-z[1])^2 + z[2]^2;
  return c, r2;
end function;

function ConjugateEllEll(U,V,plc)
  fU, plc:= fixed_exact(U, plc);
  fV, plc:= fixed_exact(V, plc);
  if fU[1] eq fV[1] then
    P1:= [* fU[1] , oo *];
  else
    //x:= (fV[2]^2 - fU[2]^2 + fV[1]^2 - fU[1]^2) / (2*(fV[1] - fU[1]));
    //r, plc:= mysqrt( (x-fU[1])^2 + fU[2]^2, plc);
    x, r2:= Geodesic(fU, fV);
    r, plc:= mysqrt(r2, plc);
    P1:= [x-r, x+r];
  end if;
  S:= MyTransform(P1, [*0,oo*]);
  L:= ext< myNF(plc) | Polynomial([1,0,1]) >;
  fU:= Eltseq(apply(S, L ! fU));
  assert fU[1] eq 0;
  S:= DiagonalMatrix([1,fU[2]]) * S;
  U:= S*U*S^-1;
  if myEval(U[1,2], plc) lt 0 then S:= Matrix(myNF(plc), 2, [0,1,1,0]) * S; end if;
  V:= S*V*S^-1;
  if myEval(-V[1,2]/V[2,1], plc) lt 1 then S:= Matrix(myNF(plc), 2, [0,1,-1,0]) * S; end if;
  return S, plc;
end function;

function ConjugateEllPar(U,V,plc)
  fU, plc:= fixed_exact(U, plc);
  fV     := fixed_exact(V, plc)[1];
  if fV eq oo then P:= [*fU[1], oo*];
  else
    P:= [ fV-(tmp^2+ fU[2]^2)/tmp, fV ] where tmp:= fV - fU[1];
  end if;
  S:= MyTransform(P, [*0,oo*]);
  L:= ext< myNF(plc) | Polynomial([1,0,1]) >;
  z:= apply(S, L ! fU);
  assert z[1] eq 0;
  S:= DiagonalMatrix([1,z[2]]) * S;
  U:= S*U*S^-1;
  if myEval(U[1,2], plc) lt 0 then S := Matrix(myNF(plc), 2, [1,0,0,-1]) * S; end if;
  return S, plc;
end function;

function ConjugateEllHyp(U,V,plc)
  fU, plc:= fixed_exact(U, plc);
  fV, plc:= fixed_exact(V, plc);
  if fV[1] cmpeq oo or fV[2] cmpeq oo then
    t:= fV[1] cmpeq oo select fV[2] else fV[1];
    w:= fU;
    w[1]:= 2*t - w[1];
    m, r2:= Geodesic(fU, w);
    r, plc:= mysqrt(r2, plc);
  else
    M:= (fV[1] + fV[2])/2; R2:= (fV[1]-fV[2])/2;
    m:= (M^2 - R2 - l2norm(fU)) / (2*(M-fU[1]));
    r, plc:= mysqrt((M-m)^2 - R2, plc);
  end if;
  S:= MyTransform([* fV[1], fV[2], m-r *], [-1,1,0]);
  L:= ext< myNF(plc) | Polynomial([1,0,1]) >;
  z:= apply(S, L ! fU);
  assert z[1] eq 0;
  S:= DiagonalMatrix([1,z[2]]) * S;
  U:= S*U*S^-1;
  if myEval(U[1,2], plc) lt 0 then S:= Matrix(myNF(plc), 2, [0,1,1,0]) * S; end if;
  V:= S*V*S^-1;
  if myEval(V[1,2], plc) lt 0 then S:= Matrix(myNF(plc), 2, [0,1,-1,0]) * S; end if;
  return S, plc;
end function;

function ConjugateParPar(U,V,plc)
  fU:= fixed_exact(U, plc)[1];
  fV:= fixed_exact(V, plc)[1];
  S:= MyTransform([* fU, fV *], [* oo, 0 *]);
  U:= S*U*S^-1;
  if myEval(U[1,2], plc) lt 0 then S:= DiagonalMatrix([myNF(plc) | -1,1]) * S; end if;
  return S, plc;
end function;

function ConjugateParHyp(U,V,plc)
  fU:= fixed_exact(U, plc)[1];
  fV, plc:= fixed_exact(V, plc);
  S:= MyTransform([* fU, fV[1], fV[2] *], [* oo, -1, 1 *]);
  U:= S*U*S^-1;
  if myEval(U[1,2], plc) lt 0 then S:= DiagonalMatrix([myNF(plc) | -1,1]) * S; end if;
  return S, plc;
end function;

function ConjugateHypHyp(U,V,plc)
  fU, plc:= fixed_exact(U, plc);
  fV, plc:= fixed_exact(V, plc);
  s, plc:= mysqrt(cross(fU[1], fU[2], fV[1], fV[2]), plc);
  k:= (1+s)/(1-s);
  S:= MyTransform([* fU[1], fU[2], fV[1] *], [* -1, 1, -k *]);
  U:= S*U*S^-1;
  if myEval(U[1,2], plc) lt 0 then S:= DiagonalMatrix([myNF(plc) | -1,1]) * S; end if;
  return S, plc;
end function;

/* Membership for free products with tr[A,B] > -2 */
function MembershipFP(G,M)
  t:= G`type;
  tU:= Trace(G`U) eq 2 select 1 else t[1];
  tV:= Trace(G`V) eq 2 select 1 else t[2];
  if not assigned G`S then
    U:= G`U; V:= G`V; plc:= G`plc;
    if t[1] ne 0 and t[2] ne 0 then
      S, plc:= ConjugateEllEll(U,V,plc);
    elif t[1] ne 0 and tV eq 1 then
      S, plc:= ConjugateEllPar(U,V,plc);
    elif t[1] ne 0 then
      S, plc:= ConjugateEllHyp(U,V,plc);
    elif tU eq 1 and tV eq 1 then
      S, plc:= ConjugateParPar(U,V,plc);
    elif tU eq 1 then
      S, plc:= ConjugateParHyp(U,V,plc);
    else
      S, plc:= ConjugateHypHyp(U,V,plc);
    end if;
    G`S:= S; G`plc:= plc; 
  end if;
  plc:= G`plc; S:= G`S; SI:= S^-1; U:= S*G`U*SI; V:= S*G`V*SI; M:= S*M*SI;

  L:= ext< myNF(plc) | Polynomial([1,0,1]) >;
  if t[1] ne 0 and t[2] ne 0 then
    u:= mysqrt(-V[1,2]/V[2,1], plc);
    z0:= (1+u)/2 * L.1;
//    u:= Sqrt(myEval(-V[1,2]/V[2,1], plc));
//    z0:= RationalReconstruction(1, u)/2 * L.1;
  elif t[1] ne 0 then
    z0:= 2*L.1;
  else
    z0:= L.1;
  end if;

  z:= apply(M, z0);
  F:= FreeGroup(2);
  w:= F ! 1;

  repeat
    k, z:= GetExp(z, U, tU, plc);
    if z eq 0 then break; end if;		// we hit the boundary
    w := w * F.1^-k;
    k, z:= GetExp(z, V, tV, plc);
    if z eq 0 then break; end if;
    w := w * F.2^-k;
  until k eq 0;

  return z eq z0, w;
end function;

function get_exponent(z, c, xy, r, plc: S:= false, Tries:= 3)
  if S cmpne false then z:= apply(S[1], z); end if;
  l:= l2norm(z);
  czl:= c*z[1]/l;
  i:= 0; done:= false;
  while not done and i lt Tries do
    prec:= 300 + 100*i;
    try
      tmp:= myEval(czl, plc: ForceReal, Precision:= prec) + Sqrt( myEval(czl^2 - xy/l, plc: ForceReal, Precision:= prec) );
      tmp:= Log(tmp) / Log(myEval(r, plc : ForceReal, Precision:= prec));
      k:= Ceiling(tmp);
      done:= true;
    catch e
      "rounding error, trying to recover";
      k:= 0;
    end try;
    i+:= 1;
  end while;
  if k ne 0 then
    z *:= r^k;
    if S cmpne false then z:= apply(S[2], z); end if;
  end if;
  return k, z;
end function;

// The case that tr[A,B] <= -2
function Membership1(G, M)
  if not assigned G`S then
    _, S, G`plc:= diagonalize(G`U, G`plc);
    G`S:= DiagonalMatrix([(S*G`V*S^-1)[2,1], 1]) * S;
  end if;
  plc:= G`plc; S:= G`S; SI:= S^-1;
  U:= S * G`U * SI; V:= S * G`V * SI; VI:= V^-1;

  // construct domain and init
  a,b,_,d:= Explode(Eltseq(V)); r:= U[1,1]^2;
  R2P:= a/d; R2M:= d/a;
  rR:= myEval(r, plc: ForceReal);
  s, plc1:= mysqrt(R2P, plc);

  if Trace(G`A*G`B*G`A^-1*G`B^-1) eq -2 then
    x:= a-s; y:= -d+1/s;
    assert plc eq plc1;
  else // We can get away without using s. This usually safes one relative extension.
    sR:= myEval(s, plc1);
    x:= RationalReconstruction( (aR + sR)/rR, aR - sR ) where aR:= myEval(a, plc);
    y:= RationalReconstruction( dR + 1/sR, (dR - 1/sR)/rR ) where dR:= myEval(-d, plc);
  end if;
  xy:= x*y;

  VV, SV, plc1:= diagonalize(V, plc1);
  SVI:= SV^-1;
  kk:= VV[1,1]^2;
  xV:= apply(SV, -d+1/s); yV:=apply(SV, -d-1/s); cV:= (xV+yV)/2; xVyV:= xV*yV;

  // the inner circle is centered at (c,  0) with radius r1
  // the outer circle is centered at (r*c,0) with radius r*r1
  c:= (x+y)/2; r1:= (x-y)/2;

  L := ext< myNF(plc ) | Polynomial([1,0,1]) >;           // plc  is real, so this works.
  L1:= ext< myNF(plc1) | Polynomial([1,0,1]) >;           // plc1 is real, so this works.
  t:= L.1 * RationalReconstruction(s, rR*s) where s:= Sqrt(myEval(-xy, plc));
  assert InCircle(t, r*c, (r*r1)^2, plc) and not InCircle(t, c, r1^2, plc);

  // start algorithm
  z:= apply(S*M*SI, t);
  F:= FreeGroup(2);
  w:= F ! 1;

  repeat
// TODO: boundary tests
    // get power of U
    k, z:= get_exponent(z, c, xy, r, plc: Tries:= oo);
    w:= w*F.1^-k;
//    "applied U^", -k;

    k:= 0;
    while InCircle(z, -d, R2M, plc) and  k le MAXTRIES do k+:= 1; z:= apply(V,  z); end while;
    while InCircle(z,  a, R2P, plc) and -k le MAXTRIES do k-:= 1; z:= apply(VI, z); end while;
    if Abs(k) gt MAXTRIES then
      kV, zz:= get_exponent(L1 ! Eltseq(z), cV, xVyV, kk, plc1: S:= [SV, SVI]);
      z:= L ! Eltseq(zz);
      k +:= kV;
    end if;
    w:= w*F.2^-k;
//    "applied V^", -k;
  until k eq 0;

  return z eq t, w;
end function;

// membership for elementary discrete groups
function MembershipED(G,M)
  assert G`type eq "ed";
  U:= G`U; V:= G`V; M:= G`S*M*G`S^-1; plc:= G`plc;
  wU, wV:= Explode(G`WordsUV);
  F:= FreeGroup(2);
  if Trace(U) in {-2,2} then		// all parabolic
    if M[1,1] notin {U[1,1],V[1,1]} then return false, _; end if;
    // make M[1,1] = +1
    if M[1,1] eq 1 then w:= F!1;
    elif U[1,1] eq -1 then
      w:= F.1; M:= U^-1 * M;
    else
      w:= F.2; M:= V^-1 * M;
    end if;

    // compute the subgroup of G with entries +1
    if U[1,1] eq -1 and V[1,1] eq -1 then
      V:= U*V; wV:= wU*wV;
      U:= U^2; wU:= wU^2;
    elif U[1,1] eq -1 then
      U:= U^2; wU:= wU^2;
    elif V[1,1] eq -1 then
      V:= V^2; wV:= wV^2;
    end if;

    Mat:= Matrix( [ Eltseq(U[1,2]), Eltseq(V[1,2]) ]  );
    d:= Denominator(Mat);
    Mat:= Matrix(ZZ, d*Mat);
    ok, sol:= IsConsistent( Mat, Vector(ZZ, Eltseq(d*M[1,2])) );
    if not ok then return false, _; end if;
    return true, w * F.1^sol[1] * F.2^sol[2];

  elif IsOne(V^4) then		// V has order 1,2,4
    a:= U[1,1]; plc:= G`plc;
    for i in [0..3] do
      X:= V^-i * M;
      if {X[1,2], X[2,1]} ne {0} then continue; end if;
      b:= X[1,1];
      l:= Round(Log(Abs(myEval(a, plc)), Abs(myEval(b, plc))));
      if a^l eq b then return true, F.2^i * F.1^l; end if;
    end for;
    return false, _;

  else				// all hyperbolic
    sM:= mySign(M[1,1], plc);
    sU:= mySign(U[1,1], plc);
    sV:= mySign(V[1,1], plc);
    if not IsDiagonal(M) or sM notin {sU, sV} then return false, _; end if;
    if sM eq 1 then w:= F!1;
    elif sU eq -1 then
      w:= F.1; M:= U^-1 * M;
    else
      w:= F.2; M:= V^-1 * M;
    end if;

    if sU eq -1 and sV eq -1 then
      V:= U*V; wV:= wU*wV;
      U:= U^2; wU:= wU^2;
    elif sU eq -1 then
      U:= U^2; wU:= wU^2;
    elif sV eq -1 then
      V:= V^2; wV:= wV^2;
    end if;

    a:= U[1,1]; b:= V[1,1];
    r:= BestApproximation(Log(myEval(b, plc), myEval(a, plc)), 10^20);
    n:= Numerator(r); m:= Denominator(r);
    _, k,l:= XGCD(n,m);
    lambda:= a^k * b^l;
    e:= Round(Log(myEval(lambda, plc), myEval(M[1,1], plc))); 
    if M eq (U^k * V^l)^e then return true, w * (F.1^(k*e) * F.2^(l*e)); end if;
    return false, _;
  end if;
end function;

// membership for triangle groups
function MembershipTriangle(G, M)
  if not assigned G`S then
    S, plc:= ConjugateEllEll(G`U, G`V, G`plc);
    G`S:= S; G`plc:= plc;
  end if;
  S:= G`S; plc:= G`plc; SI:= S^-1;
  U:= S*G`U*SI; V:= S*G`V*SI; M:= S*M*SI;
  u:= mysqrt(-V[1,2]/V[2,1], plc);
  L:= ext< myNF(plc) | Polynomial([1,0,1]) >;
  z0:= (1+u)/2 * L.1;
  z:= apply(M, z0);
  F:= FreeGroup(2);
  w:= F ! 1;

  Gens:= [ U,V,U^-1,V^-1 ];
  Words:= [ F.1^-1, F.2^-1, F.1, F.2 ];
  d:= myEval(HypDistance( z0, z ), plc);
  while true do
    P:= [ apply(g, z) : g in Gens ];
    D:= [ myEval(HypDistance( z0, p ), plc) : p in P ];
    m, i:= Min(D);
    if m gt d then break; end if;
    d:= m;
    w:= w * Words[i];
    z:= P[i];
  end while;

  return z eq z0, w;
end function;

//Does the geodesic (x,r) intersect the segment [z1,z2]?
function DoesIntersect(x, r, z1, z2, plc)
  if r eq 0 then
    x1:= myEval(z1[1], plc);
    x2:= myEval(z2[1], plc);
    x:= myEval(x, plc);
    return (x1 le x) xor (x2 le x);
  else
    return InCircle(z1, x, r, plc) xor InCircle(z2, x, r, plc);
  end if;
end function;

// Map the geodesic through z and zz to i*R
function GeodesicToiR(z, zz, plc)
  x,r:= Geodesic( z, zz );
  if r eq 0 then
    return Matrix( 2, [ 1, -x, 0, 1]), plc;
  end if;
  s, plc:= mysqrt(r, plc);
  P:= [* x-s, x+s *];
  S:= MyTransform(P, [* 0, oo *]);
  if myEval(Determinant(S),  plc) le 0 then S[1] *:= -1; end if;
  return S, plc;
end function;

function Test(fp, z, plc)
  for i in [1..3] do
    j,k:= Explode([j: j in [1..3] | j ne i ]);
    x,r:= Geodesic(fp[j], fp[k]);
    if InCircle(fp[i], x,r, plc) ne InCircle(z,x,r,plc) then return false, i; end if;
  end for;
  return true, 0;
end function;


// This function is much slower than the others, but then this case is quite rare.
// It would be much faster to use inexact arithmetic ...
function MembershipComm(G, M)
  Mold:= M;
  plc:= G`plc;
  C:= -comm(G`A, G`B);
  p,n:= MyOrder(C, plc);
  assert n in {1,2};
  assert n eq 1 or IsOdd(p);
  f, plc:= fixed_exact(G`A, plc);
  g, plc:= fixed_exact(G`B, plc);
  assert #f eq 2 and #g eq 2;
  A:= Matrix(myNF(plc), G`A); B:= Matrix(myNF(plc), G`B);
  U:= MyTransform( [* f[1], f[2], g[1] *], [* f[2], f[1], g[2] *] );
  V:= U*G`A;
  W:= G`B*U;

  // let F be the triangle with vertices fp[i] (i=1,2,3).
  fp:= [];
  fp[1], plc:= fixed_exact(W*U*V, plc);
  fp[2], plc:= fixed_exact(U*V*W, plc);
  fp[3], plc:= fixed_exact(V*W*U, plc);

  // find a point in F^o (without enlarging the field).
  i:= [ i: i in [2..3] | fp[1,1] ne fp[i,1] ][1];
  x:= RationalReconstruction( myEval(fp[1,1], plc), myEval(fp[i,1], plc));
  j:= Rep({1..3} diff {1,i});
  T:= DoesIntersect(x,0, fp[1], fp[j], plc)  select < 1,j > else < i,j >;
  Y:= [];
  for t in [T, <1,i>] do
    xx,r2:= Geodesic( fp[t[1]], fp[t[2]]); assert r2 ne 0;
    Append(~Y, Sqrt(myEval(r2 - (x-xx)^2, plc: ForceReal)));
  end for;
  y:= RationalReconstruction(Y[1], Y[2]);
  // so z0=ix+iy in F^o.

  L:= ext< myNF(plc) | Polynomial([1,0,1]) >;
  z0:= L ! [x,y];
  assert Test(fp, z0, plc);
  z:= apply(M, z0);
  S, plc:= GeodesicToiR(z0, z, plc);
  L:= ext< myNF(plc) | Polynomial([1,0,1]) >;
  z0:= apply(S, L ! [ myNF(plc) | z0[1], z0[2]]);
  z := apply(S, L ! [ myNF(plc) | z[1] , z[2] ]);
  if myEval(z0[2], plc) gt myEval(z[2], plc) then
    T:= Matrix(BaseRing(S), 2, [0,1,-1,0]);
    z0:= apply(T, z0);
    z:= apply(T, z);
    S:= T*S;
  end if;

  SI:= S^-1; U:= S*U*SI; V:= S*V*SI; W:= S*W*SI; M:= S*M*SI;

  // The triangle with vertices fp[i] (i=1,2,3) is the fundamental domain for <U,V,W>
  fp:= [ L | ];
  fp[1], plc:= fixed_exact(W*U*V, plc);
  fp[2], plc:= fixed_exact(U*V*W, plc);
  fp[3], plc:= fixed_exact(V*W*U, plc);
  
  F:= FreeGroup(3);
  fold:= fp;
  h:= hom< F -> GL(2, myNF(plc)) | [U,V,W] >;
  w:= F.0;

  t0:= myEval(z0[2], plc);
  t := myEval(z[2] , plc);
  assert t0 le t;
  N:= M^0;

  // We move the triangle F to the point z.
  // This is slightly different from Alg. 6, were we mobe p to F.
  // But this way, we can keep the geodesic [z,z0] as as subset of i*R
  while true do
    // check which sides J of the triangle F intersect i*R
    // start with a vertex on i*R
    I:= [i: i in [1..3] | fp[i,1] eq 0 ];
    if #I ne 0 then
      E:= [ myEval( fp[i,2], plc) : i in I ];
      m, i:= Max(E); i:= I[i];
      if m ge t then break; end if;	// done yet?
      J:= Remove([1,2,3], i);
      if #I eq 1 then
        y:= Sqrt(myEval(r2-x^2, plc: ForceReal)) where x,r2:= Geodesic(fp[J[1]], fp[J[2]]) ;
        // If the vertex is not on top, we use the generic routine.
        if E[1] le y then i:= 0; end if;
      end if;
      if i ne 0 then
        if IsEven(p) then
          g:= (F.J[1] * F.i * F.J[2])^(p div 2); 
        else
          g:= F.0;
          // There are 6 candidates of gF \cap F = vertex on iR and gF \cap [z,z0] \ne {}.
          for gg in [ F.0, F.J[1], F.J[2] ], l in [(p-1) div 2, (p+1) div 2] do
            ggg:= gg * (F.J[1] * F.i * F.J[2])^l;
            f:= [ apply(N*h(ggg) , f): f in fold ];
            II:= [ i: i in [1..3] | f[i] ne fp[i] ];
            S:= { Sign(myEval(f[j, 1], plc)): j in II };
            if #S eq 2 then g:= ggg; end if;
          end for;
          assert g ne F.0;
        end if;
        w:= w*g;
        N:= N*h(g);
        fp:= [ apply(N, f): f in fold ];
        continue;
      end if;
    end if;
    // the generic case, no vertex on iR or the vertex is not on top of the triangle.
    S:= [ Sign(myEval(f[1], plc)): f in fp ];
    assert #Set(S) eq 2;
    s:= &*S;
    i:= Index(S, s);
    J:= Indices(S, -s);
    Y:= [ Sqrt(myEval(r2-x^2, plc: ForceReal)) where x,r2:= Geodesic(fp[i], fp[j]) : j in J ];
    y, m:= Max(Y); j:= J[m];
    // are we done yet?
    if y ge t then break; end if;
    k:= Rep({1..3} diff {i,j});
    w:= w * F.k;
    N:= N * h(F.k);
    fp:= [ apply(N, f): f in fold ];
  end while;

  if N notin {M, -M} then return false, _; end if;

  // Now rewrite w as a word in A and B.
  // set U=f.3, V=U.A = F.3*F.1 and W=B*U=F.2*F.3
  w:= f(w) where f:= hom< F -> F | [ F.3, F.3*F.1, F.2*F.3 ] >;
  FF:= FreeGroup(2);
  k:= 0;
  v:= FF.0;
  for e in Eltseq(w) do
    if Abs(e) eq 3 then k := (k+1) mod 2;
    elif k eq 0 then
      v *:= FF.e;
    else
      v *:= FF.(-e);
    end if;
  end for;
  if k eq 1 then
    if n eq 1 then return false, _; end if;
    v *:= comm(FF.1, FF.2)^(p div 2) * FF.1 * FF.2;
  end if;
//  h:= hom< FF -> GL(2, G`K) | [ G`A, G`B] >;
//  assert h(v) in {Mold, -Mold};

  return true, v;
end function;

// user interface to the membership functions
intrinsic Membership(G::SL2R2Gen, M::Mtrx) -> BoolElt, GrpFPElt
{Returns true if M is in G. If this is the case, a word expressing M as the generators of G is also returned}
  require Ncols(M) eq Nrows(M) and Determinant(M) eq 1 : "The matrix must have determinant 1";
  ok, M:= CanChangeRing(M, G`K);
  require ok : "The group and the matrix must be over the same field.";
  // Do NOT remove the following line. It fills the data structure if necessary!
  require IsDiscrete(G) : "The group is not discrete.";

  if IsFinite(G) then
    GG:= G`finite[2];
    ok, M:= IsCoercible(GG, M);
    if not ok and G`PSL then 
      ok, M:= IsCoercible(GG, -M);
    end if;
    if not ok then return false, _; end if;
    _, h:= WordGroup(GG);
    return true, Evaluate(M @@ h, FreeGroup(2));
  elif IsDiscreteAndElementary(G) then
    ok, w:= MembershipED(G, M);
    if not ok and G`PSL then
      ok, w:= MembershipED(G, -M);
    end if;
    if not ok then return false, _; end if;
    // a test
    N:= w @ h where h:= hom < Parent(w) -> GL(2, G`K) | G`A, G`B >;
    assert N eq M or (G`PSL and N eq -M);
    return true, w;
  end if;

  // In the remaining cases, we first solve membership in \PSL.
  if G`type cmpeq "comm" then
    ok, w:= MembershipComm(G,M);
  elif IsTriangleGroup(G) then
    ok, w:= MembershipTriangle(G,M);
  else
    t:= Trace(G`A*G`B*G`A^-1*G`B^-1);
    if myEval(t+2, G`plc) le 0 or t eq -2 then
      ok, w:= Membership1(G, M);
    else
      ok, w:= MembershipFP(G, M);
    end if;
  end if;
  if ok then 
    // w is the only possible candidate (in U,V !)
    // Rewrite w as a word in A,B
    F:= Parent(w);
    if assigned G`WordsUV then
      w:= w @ h where h:= hom< F -> F | Evaluate(G`WordsUV[1], F), Evaluate(G`WordsUV[2], F) >;
    end if;
    // Test if w does the trick
    h:= hom < F -> GL(2, G`K) | G`A, G`B >;
    N:= w @ h;
    if N eq M or (G`PSL and N eq -M) then return true, w; end if;
    if (N eq -M) and assigned G`MinusOne then
      w:= w * Evaluate(G`MinusOne, F);
      N:= w @ h;
      assert N eq M;
      return true, w;
    end if;
  end if;
  return false, _;
end intrinsic;

intrinsic WitnessPair(G::SL2R2Gen)-> Mtrx, Mtrx, []
{Returns a witness pair for G as well as SLPs expressing the pair in terms of the generators}
  prepare(G);
  require assigned G`U : "The group has no witness pair";
  return G`U, G`V, G`WordsUV;
end intrinsic;


// These functions were used to produce the tables in Section 10.
//
// Simply copy the functions from the beginning of the file
// (mysqrt etc.) and the two functions below into a magma shell.

// Times the membership test for G with 10 words of length l.
function timemember(G, l)
  F:= FreeGroup(2);
  h:= hom< F -> GL(2, G`K) | G`A, G`B >;
  W:= [];
  for i in [1..10] do
    w:= F ! 1;
    while #w lt l do
      w*:= F.Random([-2,-1,1,2]);
    end while;
    Append(~W, w);
  end for;
  //[ #w: w in W ];
  W:= [ w @ h : w in W ];
  
  for w in W  do
    time ok, ww:= Membership(G, Matrix(w));
    assert ok;
  end for;
  return 0;
end function;

// The standard cos & sin matrix with trace +2*cos(pi/n).
function ConstructMatrix(n)
  assert n ge 3;

  C:= CyclotomicField(2*n); 
  K:= sub< C | (C.1 + C.1^-1)/2 >;
  cc:= Cos(Pi(RealField())/n);
  c:= K ! (C.1 + C.1^-1)/2;
  plc:= [ i: i in RealPlaces(K) | Abs(Evaluate(c, i) - cc) lt 10^-10 ][1];
  assert Abs(Evaluate(c, plc) - cc) lt 10^-10;
  s, plc:= mysqrt(1-c^2, plc);
  assert AbsoluteDegree(NumberField(plc)) eq EulerPhi(n);
  M:= Matrix(2, [c,s,-s,c]);
  assert IsOne(-M^n);
  return M,plc;
end function;

/*  
 // These lines were used to produce Figure 5.
 K:= NumberField(Polynomial([-1,0,-2,0,1,0,1]));
 a:= K.1;
 A:= Matrix(K, 2, [0,1,-1,0]);
 B:= 1/2 * Matrix(2, [1, -a^5+a^4+a^2+2*a-1, -a^5-a^4-a^2+2*a+1, 1 ]);
 cc:= Sqrt(2*Cos(2*Pi(RealField())/7));
 plc:= [ i: i in RealPlaces(K) | Abs(Evaluate(a, i) - cc) lt 10^-10 ][1];
 G:= TwoGeneratorSubgroupSL2R(A,B,plc);
 timemember(G, 100);		// replace 100 by different values to test larger words.


 // These lines were used to produce Figure 6.
 n:= 5;
 X, plc:= ConstructMatrix(n);
 Y:= Matrix(2, [1, Ceiling(2*Cot(Pi(RealField())/(2*n)) + 1), 0, 1]);
 G:= TwoGeneratorSubgroupSL2R(X,Y,plc);
 timemember(G, 100);

 // An example with commutator [A,B] of order 5.
 C:= CyclotomicField(5);
 K:= sub< C  | C.1 + C.1^-1 >;
 plc:= RealPlaces(K)[2];
 alpha:= -K.1/(5-1/5);
 a:= 1;
 d:= alpha^2+1;
 X:= Matrix(K, 2, [0,1,-1,0]);
 Y:= Matrix(K, 2, [0,-5,1/5,0]);
 Z:= Matrix(K, 2, [alpha,-d,a,-alpha]);
 A:= X*Y;
 B:= Z*X;
 G:= TwoGeneratorSubgroupSL2R(A,B,plc);
*/
