# remember to first load Farey Symbols:
# LoadPackage("congruence");  

# it will be useful to find the maximum value of the labels
# though if space is not a problem, this could just return
# the Length of the labels.

max_label:= function(L)
  local s, i;
  s:=1;
  for i in [1..Length(L)] do
    if (not L[i] = "even") and (not L[i] = "odd")  and  L[i] > s then
       s := L[i];
    fi;
  od;
  return s;
end;;


# For a list of labels L such as
# [1,3,4,7,4,7,3,1,"odd","even"], for reference, indices are:
#  1 2 3 4 5 6 7 8  9     10
# want to produce a list:
# [[9],[10],[1,8],[],[2,7],[3,5],...]
# this is the list of the form:
# [[all indices with L[x] = "odd"],[all indices with L[x] = "even"],
# [all indices with L[x] = 1], ....]

# assume L is a list of integers, or "odd" or "even".

edgepairs := function(L)
  local max, pairs, i;
  pairs:=[];
  max:=max_label(L);
  for i in [1..max+2] do
      pairs[i] := [];
  od;
  for i in [1..Length(L)] do
      if L[i]="odd" then
        Add(pairs[1],i);
      elif L[i]="even" then
        Add(pairs[2],i);
      else
        Add(pairs[L[i]+2],i);
      fi;
  od;
  return pairs;
end;;

# for each edge of a Farey Symbol, we compute the generator
# which maps that edge to another edge.
# (this is done at the same time as the fundamental
# domain is computed, but the data may not have been stored,
# and has to be recomputed; suggest change for a future version)

# this function gives "edge gluing matrices" as a number in the
# list of generators (gens); negative entries mean the inverse matrix,
# e.g., -5 would mean (5th generator)^(-1)
# (note, the list of labels in a Farey sequence says which edge is
# glued to which; -2 and -3 means there is an elliptic point order
# 2 or 3)
#
#
# the input is assumed to be a FareySymbol;
# another version of this function could take input to be the group
#
# Note, if the output of this function was
# stored as an attribute of the FareySymbol,
# then it would not have to be recomputed
#



gluing_matrices := function(FS)
   local cusps, gens, label_list, glue_list, l, i, index, gfs, labels, matrix;
   # the following is a list of the cusps of the sequence,
   # and other data extracted from the FareySymbol
   gfs := GeneralizedFareySequence(FS);
   labels := LabelsOfFareySymbol(FS);
   gens := GeneratorsByFareySymbol( FS );
   # make a list of which edges have a given label:
   label_list := edgepairs(labels);
   # the following list will be what is finally returned, 
   # a list of integers as described above.
   glue_list := [];
   # make list of which generator joins two edges,
   # in the non elliptic case
   for i in [3..Length(label_list)] do
      l := label_list[i];
      matrix := MatrixByFreePairOfIntervals( gfs, l[1], l[2] );
      index := PositionNthOccurrence( gens ,matrix,1);
      if index = "fail" then
        index := -PositionNthOccurrence(gens,matrix^(-1),1);
      fi;
      glue_list[l[1]] := index;
      glue_list[l[2]] := -index;
   od;
   # Now deal with elliptic elements:
   for i in label_list[1] do
      matrix := MatrixByOddInterval( gfs, i );
      index := PositionNthOccurrence(gens,matrix,1);
      if index = "fail" then
         index := -PositionNthOccurrence(gens,matrix^(-1),1);
         glue_list[i] := -index;
      else
         glue_list[i] := -index;
      fi;
   od;
   for i in label_list[2] do
      matrix := MatrixByEvenInterval( gfs, i );
      index := PositionNthOccurrence(gens,matrix,1);
      if index = "fail" then
         index := -PositionNthOccurrence(gens,matrix^(-1),1);
         glue_list[i] := -index;
      else
         glue_list[i] := -index;
      fi;
   od;
   return glue_list;     
end;;
   

# following function determines which "gap" an image ImL of
# a domain, given by a list of cusps L, belongs to.
# (cusps means rational or infinity)
#
# Assume that L always has the form ["infinity",..."infinity"]
# if not this program will not work correctly.
# L is a generalized Farey sequence
#
# The function either returns a index of a "gap",
# which is a number between 1 and #L-1,
# or it returns 0, meaning that the domains ImL and L are equal
# or it returns "overlap" meaning that there is overlap, but not equality.

which_gap := function(L,ImL)
   local finite_cusps_ImL,finite_cusps_L,i,maxL, maxImL,minL,minImL, start, Max, Min, index, found;
   # need a list of cusps not including infinity:   
   finite_cusps_ImL := [];
   finite_cusps_L := [];
   for i in [1..Length(ImL)] do
      if not ImL[i] = infinity then
        Add(finite_cusps_ImL,ImL[i]);
      fi; 
      if not ImL[i] = infinity then
        Add(finite_cusps_L,L[i]);
      fi; 
   od;
   maxImL := Maximum(finite_cusps_ImL);
   minImL := Minimum(finite_cusps_ImL);
   maxL := Maximum(finite_cusps_L);
   minL := Minimum(finite_cusps_L);
   if minImL >= maxL then
      return Length(L)-1;
   elif maxImL <= minL then
      return 1;
   elif minImL = minL and maxImL = maxL then
      if ImL = L then
         return 0;
      else 
         return "overlap";
      fi;
   else
      i := 1;
      found := false;      
      while not found and i < Length(finite_cusps_ImL) do
         i := i + 1; 
         if maxImL <= finite_cusps_L[i] then
            found := true;
            index := i;
         fi;
      od;
   fi;
   if minImL >= finite_cusps_L[i-1] then
      return index-1;
   else
      return "overlap";
   fi;
end;;

# Need to be able to apply action of matrices to cusps

fractionallineartransformation:= function(g,c)
  local den, num;
  if c = infinity then
    if g[2][1] = 0 then
      return infinity;
    else    
      return g[1][1]/g[2][1];
    fi;  
  else
    num:=g[1][1]*c + g[1][2];
    den:=g[2][1]*c + g[2][2];
    if den = 0 then
      return infinity;
    else
      return num/den;
    fi;
  fi;
end;;

PSL2multiply := function(g,L)
  local imL, i;
  imL := [];
  for i in [1..Length(L)] do
    Add(imL,fractionallineartransformation(g,L[i]));
  od;
  return imL;
end;;


# finally, we can give an algorithm to determine a word for
# a given matrix g in G in terms of the generators:
#
# example of using this function:
#
#  Read("membershiptest1.g");
#  G16:=Gamma0(16);
#  FS16:=FareySymbol(G16);
#  glue_list:=gluing_matrices(FS16);
#  g:=Random(G16);
#  find_word(FS16,glue_list,g);

find_word := function(FS,glue_list,g)
   local gens, L, ImL, done, word,letter, gap,i;
   gens := GeneratorsByFareySymbol( FS );
   L := GeneralizedFareySequence( FS );
   ImL := PSL2multiply(g,L);
   word:=[];
   done := false;
   while not done do;      
      gap := which_gap(L,ImL);
      if gap = "overlap" then
         Print("matrix not in group\n");
         return word;
      elif gap = 0 then
         done := true;
      else
         # get next "letter" in the word for the matrix:
         letter := glue_list[gap];
         Add(word,-letter);
         ImL:=PSL2multiply(gens[AbsoluteValue(letter)]^(SignInt(letter)),ImL);
      fi;
   od;
   return word;
end;;


# the following function is for testing purposes

# gens is a list of generators,
# "word" a sequence of integers, none
# of which is bigger than the size of the list
# of generators.  
# a words [4,6,-3] will return the product
# gens[4]*gens[6]*gens[3]^(-1);

multiply_out_word:=function(gens,word)
   local g, i;
   g := [[1,0],[0,1]];
   for i in word do
      g := g*gens[AbsoluteValue(i)]^SignInt(i);
   od;
   return g;
end;;