Vhakendata := function(length) 
local F, fgens,x,y,Dx,Dxi,Dy,Dyi,G,ggens,dx,dy,t,rels,H,hgens,
 Hauto, Kauto, words, goodwords, vect, countervect, Isreduced,
 spot, counterspot, mats, mat, J, Ind, K, P, kgens, word, exp1,  
 I, firstletter, twist1, twistdata, twist, exp2, originalexp1, originalexp2,
 originalword, i,j, underword, old_exp1, old_exp2, old_word, ind,n_three,
 n_four, Spot, w, v, data, perm, dup, list, conjugates, gw, ggoodwords,matt,
twistedslope1,twistedslope,Slope,
 pvar,qvar,svar,p,q,s,Newmat3,Mx3,My3, Newmat4,Mx4,My4;    
twistdata := [];



data := [];

G := FreeGroup(2); 
ggens := FreeGeneratorsOfFpGroup(G); 
dx := ggens[1];   
dy := ggens[2];   
t := (dx*dy)^3;  
rels := [t*(dx*dy*dx)^-2, t^2, t*dx*t^-1*dx^-1, t*dy*t^-1*dy^-1]; 
H := G/rels;     
hgens := GeneratorsOfGroup(H); 
dx := hgens[1];   
dy := hgens[2];   


## Alexander polynomial functions

pvar := Indeterminate(Rationals,"p");
p := LaurentPolynomialByCoefficients(FamilyObj(pvar),[pvar],0);
qvar := Indeterminate(Rationals,"q");
q := LaurentPolynomialByCoefficients(FamilyObj(qvar),[qvar],0);
svar := Indeterminate(Rationals,"s");
s :=  LaurentPolynomialByCoefficients(FamilyObj(pvar),[pvar^0],1);

Mx3 := [[s,s^-1,-s^-1,s-s], [s-s,s^0,s-s,s-s],
 [s-s,s-s^-1,s^-1,s-s],[s-s,s-s,s-s,s^0]];



My3 := [[s^0,s-s,s-s,s-s],[s-s,s^0,s-s,s-s],
[s^0+s+s^2,s-s,s^0,s-s],[s^0,s-s,s-s,s^0]];



Mx4 :=
  [[s,-1-s^-1,s-s,s-s, s^-1, s],
 [s-s,s^-1,s-s,s-s, s^0-s^-1, s-s],
 [s-s,s-s,s^0,s-s,s-s,s-s],
 [s-s,-s^-1,s-s,s^0,s^-1,s-s],
 [s-s,s-s,s-s,s-s,s^0,s-s],
 [s-s,s-s,s-s,s-s,s-s,s^0]];

My4 := 
 [[s^0,s-s,s-s,s-s, s-s,s-s],
 [1+s,s^0,s-s,s-s, s-s, s-s],
 [s^0,s-s,s^0,s-s,s-s,s-s],
 [s-s,s-s,s-s,s^0,s-s,s-s],
 [s-s,s-s,s-s,s-s,s^0,s-s],
 [s-s,s-s,s-s,s-s,s-s,s^0]];




Newmat3 := function(mat)
local x,w,z;
x := [[s^0,s-s,s-s],[s-s,s-s,s^0],[s-s,s-s,s-s]];
w := x - mat{[1,3,4]}{[1,2,3]};
z:= [[s-s,s-s,s-s,s-s],[s-s,s-s,s-s,s-s],[s^-1-s^0,s-s,s-s,s-s],[p,s-s,q,s-s]];
z{[1,2,3]}{[2,3,4]} := w;
return z;
end;

Newmat4 := function(mat)
local x,w,z;
x := [[s^0,s-s,s-s,s-s,s-s,s-s],
[s-s,s^0,s-s,s-s,s-s,s-s],
[s-s,s-s,s^0,s-s,s-s,s-s],
[s-s,s-s,s-s,s^0,s-s,s-s],
[s-s,s-s,s-s,s-s,s^0,s-s],
[s-s,s-s,s-s,s-s,s-s,s^0]];
w := x - mat;
z:= [[s-s,w[1][1], w[1][2],w[1][5],w[1][6]],
[s-s,w[2][1],w[2][2],w[2][5],w[2][6]],
[s^0-s,w[3][1],w[3][2],w[3][5],w[3][6]],
[p, p*w[4][1],p*w[4][2], p*w[4][5], p*w[4][6]+q],
[p,s-s,s-s,s-s,q,s-s]];
return z;
end;



Slope := function(mat)
local det,coeffs,z;
det := Determinant(mat);
coeffs := CoefficientsOfLaurentPolynomial(det);
if Length(coeffs[1]) > 0 then
z := coeffs[1][Length(coeffs[1])];
fi;
if Length(coeffs[1]) = 0 then
z := 0;
fi;
return z;
end;



# Generate monodromy words   

words := []; 
goodwords := [];  
ggoodwords := [];

vect := [];  
for i in [1 .. length] do 
Add(vect, [1,2]);  
od; 

countervect := [];     
for i in [1 .. length] do   
Add(countervect, 0); 
od;   

Spot := function(v)local i,z;
 z := length;    
for i in [1 .. length] do  
if v[i] = 0 then  
 return i;   
fi; 
 od;
for i in [0 .. length-1] do  
if v[length-i] = 1 then 
return length-i; 
fi;
od; 
return 0;  
end;   

counterspot := 1;     
while counterspot > 0 do 
countervect[counterspot] := countervect[counterspot]+1; 
if counterspot < length then
for i in [counterspot+1 .. length] do
countervect[i] := 0;  
 od; 
 fi;    


Add(words, ShallowCopy(countervect)); 

counterspot := ShallowCopy(Spot(countervect));
od; 

mats := [];
for v in words do
mat := [[1,0],[0,1]]; 
for i in [1 .. length] do 
if v[i] = 1 then
mat := mat *[[1,1],[0,1]];;
fi;
if v[i] = 2 then
mat := mat *[[1,0],[1,1]];;
fi;
od;


w := dx*dx^-1;
 for i in [1 .. Length(v)] do
 if v[i] = 1 then w := w*dx^-1;
fi; 
if v[i] = 2 then w := w*dy;
fi;
od;

gw := ggens[1]*ggens[1]^-1;
 for i in [1 .. Length(v)] do
 if v[i] = 1 then gw := gw*ggens[1]^-1;
fi; 
if v[i] = 2 then gw := gw*ggens[2];
fi;
od;


conjugates := [];
for i in [1 .. Length(gw)] do
Add(conjugates, Subword(gw, 1, i)^-1 * gw * Subword(gw, 1, i));
od;
dup := 0;
for i in [1 .. Length(gw)] do
if conjugates[i] in ggoodwords then dup := 1;
fi;
od;


if Trace(mat) > 2  
 and not mat in mats
 and dup = 0 
 then Add(goodwords, w); 
 for i in [1 .. Int(length/Length(gw))+1] do
 Add(ggoodwords, gw^i);
 od;  
 Add(mats, mat);
 fi; 


od; 


for w in goodwords do 
Add(twistdata, rec(old_word := dx*dx^-1,
  n_three := 0,
  n_four := 0,
  degree_of_cover := 0,
  word_for_f_to_the_n_three := [],
  word_for_f_to_the_n_four := [],  
  dx_exponent_for_H_three := 0, 
  dx_exponent_for_H_four := 0, 
  dy_cubed_exponent := 0,
  dy_to_the_fourth_exponent := 0,
  twist_for_beta_three := 0,
  twist_for_beta_four := 0,
  beta_three_one := [0,0],
  beta_three_two := [0,0],
  beta_four_one := [0,0],
  beta_four_two := [0,0],
 success := 0
  ));
od;





# determine n st f^n \in <dx, dy^3> or <dx, dy^4>


for j in [1 .. Length(goodwords)] do 


w := goodwords[j];


twistdata[j].old_word := w;

J := Subgroup(H, [dx, dy^3]);

Ind := function(u) 
local i; 
 for i in [1 .. 6] do 
if u^i in J then return i; 
fi; 
od; 
end;

twistdata[j].n_three := Ind(w);


J := Subgroup(H, [dx, dy^4]);


Ind := function(u) 
local i; 
 for i in [1 .. 6] do 
if u^i in J then return i; 
fi; 
od; 
end;

twistdata[j].n_four := Ind(w);

od;


# determine first slope \beta_{3,1}

for j in [1 .. Length(goodwords)]  do
old_word := twistdata[j].old_word;
n_three :=  twistdata[j].n_three;
n_four := twistdata[j].n_four;
ind := Lcm(n_four,n_three);
twistdata[j].degree_of_cover := ind;


K := Subgroup(H, [dx, dy^3,old_word^n_three]);;
P := PresentationSubgroupMtc(H,K,0);;   
TzInitGeneratorImages(P);   
TzOptions(P).protected := 2;
TzOptions(P).printLevel := 0;
 TzGoGo(P);;
kgens := GeneratorsOfPresentation(P);; 
word := TzImagesOldGens(P)[3];

exp1 := (ind/n_three)*ExponentSumWord(word, kgens[1]);
exp2 := (ind/n_three)*ExponentSumWord(word, kgens[2]);

twistdata[j].word_for_f_to_the_n_three := word;
twistdata[j].dx_exponent_for_H_three :=  exp1; 
twistdata[j].dy_cubed_exponent :=  exp2; 
 
underword :=  UnderlyingElement(old_word);
old_exp1 := ExponentSumWord(underword, ggens[1]);
old_exp2 := ExponentSumWord(underword, ggens[2]);

twist := (ind*(old_exp1+old_exp2) - (exp1+3*exp2))/12;
twistdata[j].twist_for_beta_three := twist;

twistdata[j].beta_three_one := [3,exp1 -3*twist]/Gcd(3, exp1);

## Determine second slope \beta_{3,2}



matt := MappedWord(word,[kgens[1],kgens[2]],[Mx3,My3]);
twistedslope1 :=[Value(Slope(Newmat3(matt)),[pvar,qvar],[1,0]),
Value(Slope(Newmat3(matt)),[pvar,qvar],[0,1])];
twistedslope := [-twistedslope1[2],twistedslope1[1]]/
                              Gcd(twistedslope1[1], twistedslope1[2]);
twistdata[j].beta_three_two := [twistedslope[1],-twistedslope[1]*twist+twistedslope[2]];
 
#twistedslope1 := Determinant(Newmat3(matt));
twistdata[j].beta_three_two :=  twistedslope1;


## Determine \beta_4^1


K := Subgroup(H, [dx, dy^4,old_word^n_four]);
 P := PresentationSubgroupMtc(H,K,0);   
TzInitGeneratorImages(P);   
TzOptions(P).protected := 2;
TzOptions(P).printLevel := 0;
 TzGoGo(P);;   
kgens := GeneratorsOfPresentation(P); 
word := TzImagesOldGens(P)[3];  

exp1 := (ind/n_four)*ExponentSumWord(word, kgens[1]);
exp2 := (ind/n_four)*ExponentSumWord(word, kgens[2]);

twistdata[j].word_for_f_to_the_n_four := word;
twistdata[j].dx_exponent_for_H_four :=  exp1; 
twistdata[j].dy_to_the_fourth_exponent :=  exp2; 
 

underword :=  UnderlyingElement(old_word);
old_exp1 := ExponentSumWord(underword, ggens[1]);
old_exp2 := ExponentSumWord(underword, ggens[2]);


twist := (ind*(old_exp1+old_exp2) - (exp1+4*exp2))/12;
twistdata[j].twist_for_beta_four := twist;
twistdata[j].beta_four_one := [2, exp1 - 2*twist]/Gcd(2,exp1);


## Determine \beta_4^2

matt := MappedWord(word,[kgens[1],kgens[2]],[Mx4,My4]);
#if Slope(Newmat4(matt)) <> 0 then
#twistedslope1 :=[-Value(Slope(Newmat4(matt)),[pvar,qvar],[0,1]),
#Value(Slope(Newmat4(matt)),[pvar,qvar],[1,0])];
#fi;
#if Slope(Newmat4(matt)) = 0 then
#twistedslope := [0,133];
#fi;
#twistdata[j].beta_four_two := [twistedslope[1],-twistedslope[1]*twist+twistedslope[2]];

#twistedslope1 := Determinant(Newmat4(matt));

twistedslope1 := IsConstantRationalFunction(Slope(Newmat4(matt))/(pvar*qvar));

twistdata[j].beta_four_two := twistedslope1;

 




if
(
not  twistdata[j].beta_four_one
       in  [twistdata[j].beta_three_one, twistdata[j].beta_three_two]
 and  twistdata[j].beta_four_two = (1=1)

## so surfaces from 3 and 4-fold cover suffice
)


 then  twistdata[j].success := 1;
fi;
Add(data, twistdata[j]);
od; 


return twistdata;
end; 










#failures := [];
#for i in [1 .. Length(data)] do
# if data[i].success = 0
# then Add(failures, i);
#fi;
#od;






























































































































Vhakencheck := function(v)
local F, fgens,x,y,Dx,Dxi,Dy,Dyi,G,ggens,dx,dy,t,rels,H,hgens,
 Hauto, Kauto, words, goodwords, vect, countervect, Isreduced,
 spot, counterspot, mats, mat, J, Ind, K, P, kgens, word, exp1,  
 I, firstletter, twist1, twistdata, twist, exp2, originalexp1, originalexp2,
 originalword, i,j, underword, old_exp1, old_exp2, old_word, ind,n_three,
 n_four, Spot, w, data, perm, dup, list, conjugates, gw, ggoodwords;  

twistdata := [];
  
G := FreeGroup(2); 
ggens := FreeGeneratorsOfFpGroup(G); 
dx := ggens[1];   
dy := ggens[2];   
t := (dx*dy)^3;  
rels := [t*(dx*dy*dx)^-2, t^2, t*dx*t^-1*dx^-1, t*dy*t^-1*dy^-1]; 
H := G/rels;     
 hgens := GeneratorsOfGroup(H); 
dx := hgens[1];   
 dy := hgens[2];   


# Generate monodromy words   

w := dx*dx^-1;
 for i in [1 .. Length(v)] do
 if v[i] = 1 then w := w*dx;
fi; 
if v[i] = 2 then w := w*dy;
fi;
 if v[i] = 3 then w := w*dx^-1;
fi; 
if v[i] = 4 then w := w*dy^-1;
fi;


od;

goodwords := [w];

for w in goodwords do 
Add(twistdata, rec(old_word := dx*dx^-1,
  n_three := 0,
  n_four := 0,
  degree_of_cover := 0,
  word_for_f_to_the_n_three := [],
  word_for_f_to_the_n_four := [],  
  dx_exponent_for_H_three := 0, 
  dx_exponent_for_H_four := 0, 
  dy_cubed_exponent := 0,
  dy_to_the_fourth_exponent := 0,
  twist_for_beta_three := 0,
  twist_for_beta_four := 0,
  beta_three_one := [0,0],
  beta_four_one := [0,0],
  beta_three_two := [0,0],
  beta_four_two := [0,0],
 success := 0
  ));
od;



# determine n st f^n \in <dx, dy^3> or <dx, dy^4>


for j in [1 .. Length(goodwords)] do 


w := goodwords[j];


twistdata[j].old_word := w;

J := Subgroup(H, [dx, dy^3]);

Ind := function(u) 
local i; 
 for i in [1 .. 6] do 
if u^i in J then return i; 
fi; 
od; 
end;

twistdata[j].n_three := Ind(w);


J := Subgroup(H, [dx, dy^4]);


Ind := function(u) 
local i; 
 for i in [1 .. 6] do 
if u^i in J then return i; 
fi; 
od; 
end;

twistdata[j].n_four := Ind(w);

od;


# determine exponent sum of dx and twist

for j in [1 .. Length(goodwords)]  do
old_word := twistdata[j].old_word;
n_three :=  twistdata[j].n_three;
n_four := twistdata[j].n_four;
ind := Lcm(n_four,n_three);
twistdata[j].degree_of_cover := ind;


K := Subgroup(H, [dx, dy^3,old_word^n_three]);;
P := PresentationSubgroupMtc(H,K,0);;   
TzInitGeneratorImages(P);   
TzOptions(P).protected := 2;
TzOptions(P).printLevel := 0;
 TzGoGo(P);;
kgens := GeneratorsOfPresentation(P);; 
word := TzImagesOldGens(P)[3];

exp1 := (ind/n_three)*ExponentSumWord(word, kgens[1]);
exp2 := (ind/n_three)*ExponentSumWord(word, kgens[2]);

twistdata[j].word_for_f_to_the_n_three := word;
twistdata[j].dx_exponent_for_H_three :=  exp1; 
twistdata[j].dy_cubed_exponent :=  exp2; 
 
underword :=  UnderlyingElement(old_word);
old_exp1 := ExponentSumWord(underword, ggens[1]);
old_exp2 := ExponentSumWord(underword, ggens[2]);

twist := (ind*(old_exp1+old_exp2) - (exp1+3*exp2))/12;
twistdata[j].twist_for_beta_three := twist;

twistdata[j].beta_three_one := [3,exp1 -3*twist]/Gcd(3, exp1);
twistdata[j].beta_three_two := [1,-twist];
 

K := Subgroup(H, [dx, dy^4,old_word^n_four]);
 P := PresentationSubgroupMtc(H,K,0);   
TzInitGeneratorImages(P);   
TzOptions(P).protected := 2;
TzOptions(P).printLevel := 0;
 TzGoGo(P);;   
kgens := GeneratorsOfPresentation(P); 
word := TzImagesOldGens(P)[3];  

exp1 := (ind/n_four)*ExponentSumWord(word, kgens[1]);
exp2 := (ind/n_four)*ExponentSumWord(word, kgens[2]);

twistdata[j].word_for_f_to_the_n_four := word;
twistdata[j].dx_exponent_for_H_four :=  exp1; 
twistdata[j].dy_to_the_fourth_exponent :=  exp2; 
 

underword :=  UnderlyingElement(old_word);
old_exp1 := ExponentSumWord(underword, ggens[1]);
old_exp2 := ExponentSumWord(underword, ggens[2]);


twist := (ind*(old_exp1+old_exp2) - (exp1+4*exp2))/12;
twistdata[j].twist_for_beta_four := twist;
twistdata[j].beta_four_one := [2, exp1 - 2*twist]/Gcd(2,exp1);
twistdata[j].beta_four_two := [1, - twist];

if
(
not  twistdata[j].beta_three_one
       in  [twistdata[j].beta_four_one, twistdata[j].beta_four_two]
 and  not twistdata[j].beta_three_two
      in  [twistdata[j].beta_four_one, twistdata[j].beta_four_two]

## so surfaces from 3 and 4-fold cover suffice


)


 then  twistdata[j].success := 1;
fi;

od; 


return twistdata;
end; 











































































#F := FreeGroup(2);
#fgens := FreeGeneratorsOfFpGroup(F);
#x := fgens[1];  
#y := fgens[2]; 
#Dx := GroupHomomorphismByImages(F, F, fgens, [x,y*x^-1]); 
#Dxi := GroupHomomorphismByImages(F, F, fgens, [x,y*x]);
#Dy := GroupHomomorphismByImages(F, F, fgens, [y*x,y]);   
#Dyi := GroupHomomorphismByImages(F, F, fgens, [y^-1*x,y]);

# Hauto := function(w) 
#local i,y,z;  
# y := UnderlyingElement(w);
#z := Dx*Dxi;    
# for i in [1 .. Length(y)] do 
# if Subword(y, i, i) = ggens[1] then  
#z := z*Dx;   
#fi; 
# if Subword(y, i, i) = ggens[1]^-1 then    
#z := z*Dxi;          
#fi;   
# if Subword(y, i, i) = ggens[2] then   
#z := z*Dy; 
# fi;    
#if Subword(y, i, i) = ggens[2]^-1 then     
#z := z*Dyi;     
# fi; 
# od;  
#return z;
# end;