########################################################################
#
#  Software for "Roots of Dehn twists"
#
#     by Darryl McCullough and Kashyap Rajeevsarathy
#
#  Version of May 22, 2009
#
#  Contact:  dmccullough@math.ou.edu
#
#  Written for GAP 4.4.
#
########################################################################
#
#  The output appears on the console and is also written to a file 
#  named  roots.out . To change the name of the output file, alter 
#  the String constant in the next command. 

outputfile := String("roots.out");

#
########################################################################
#
#  To find the Margalit-Schleimer roots for a given value of g:
#  
#  gap> MSRoots( 17 );
#  For g = 17 there are 8 Margalit-Schleimer roots:
#  ( 35, 0, ( 2, 2 ); ( 31, 35 ) )
#  ( 35, 0, ( 3, 19 ); ( 13, 35 ) )
#  ( 35, 0, ( 4, 13 ); ( 18, 35 ) )
#  ( 35, 0, ( 9, 23 ); ( 3, 35 ) )
#  ( 35, 0, ( 12, 17 ); ( 6, 35 ) )
#  ( 35, 0, ( 18, 34 ); ( 18, 35 ) )
#  ( 35, 0, ( 24, 33 ); ( 13, 35 ) )
#  ( 35, 0, ( 27, 32 ); ( 11, 35 ) )
#  
#  To find the number of Margalit-Schleimer roots for a given value of g:
#
#  gap> NumMSRoots(1);   
#  1
#  gap> NumMSRoots(2);   
#  2
#  gap> NumMSRoots(1108);
#  369
#  
#  To find the number of Margalit-Schleimer roots for the values of
#  g in the interval [r, r+s], use NumMSRootsRange( r, s ):
#
#  gap> NumMSRootsRange(1,4);   
#  NumMSRoots( 1 ) = 1
#  NumMSRoots( 2 ) = 2
#  NumMSRoots( 3 ) = 3
#  NumMSRoots( 4 ) = 2
#  NumMSRoots( 5 ) = 5
#  gap> NumMSRootsRange(1105,3);
#  NumMSRoots( 1105 ) = 293
#  NumMSRoots( 1106 ) = 1106
#  NumMSRoots( 1107 ) = 662
#  NumMSRoots( 1108 ) = 369
#
########################################################################
#
#  For degrees and genera of specific ( d, e ) roots:
#
#  gap> N(3,5);
#  15
#  gap> G(3,5);
#  11
#  gap> N(3,35);
#  105
#  gap> N(15,7);
#  105
#  gap> G(3,35);
#  86
#  gap> G(15,7);
#  94
#
#  To find the degrees of all (d,e) roots for a given g (i. e. on the
#  surface of genus g+1):
#
#  gap> DERoots( 53135 ); 
#  [ 53141, 53361, 55671 ]
#
#  To find them in the range m <= g <= m + n:
#
#  gap> DERootsRange(101,19);
#  DERoots( 101 ) = [ 105, 115, 123 ]
#  DERoots( 102 ) = [ 103, 105, 123 ]
#  DERoots( 103 ) = [ 105, 115 ]
#  DERoots( 104 ) = [ 105 ]
#  DERoots( 105 ) = [  ]
#  DERoots( 106 ) = [ 107, 117, 129 ]
#  DERoots( 107 ) = [ 119, 129 ]
#  DERoots( 108 ) = [ 109 ]
#  DERoots( 109 ) = [ 111, 117 ]
#  DERoots( 110 ) = [ 111, 117, 119 ]
#  DERoots( 111 ) = [  ]
#  DERoots( 112 ) = [ 113, 115, 117, 125, 135 ]
#  DERoots( 113 ) = [  ]
#  DERoots( 114 ) = [ 115 ]
#  DERoots( 115 ) = [ 117, 119, 121 ]
#  DERoots( 116 ) = [ 117, 141 ]
#  DERoots( 117 ) = [ 141 ]
#  DERoots( 118 ) = [ 119 ]
#  DERoots( 119 ) = [ 135 ]
#  DERoots( 120 ) = [ 121, 133 ]
#
#  To find the genera such that there is a (d,e) root of degree n:
#
#  gap> DERootGenera(105);
#  [ 86, 87, 92, 94, 97, 99, 100, 101, 102, 103, 104 ]
#
#  and its range version:
#
#  gap> DERootGeneraRange(105,34);
#  DERootGenera( 105 ) = [ 86, 87, 92, 94, 97, 99, 100, 101, 102, 103, 104 ]
#  DERootGenera( 107 ) = [ 106 ]
#  DERootGenera( 109 ) = [ 108 ]
#  DERootGenera( 111 ) = [ 91, 92, 109, 110 ]
#  DERootGenera( 113 ) = [ 112 ]
#  DERootGenera( 115 ) = [ 101, 103, 112, 114 ]
#  DERootGenera( 117 ) = [ 97, 106, 109, 110, 112, 115, 116 ]
#  DERootGenera( 119 ) = [ 107, 110, 115, 118 ]
#  DERootGenera( 121 ) = [ 115, 120 ]
#  DERootGenera( 123 ) = [ 101, 102, 121, 122 ]
#  DERootGenera( 125 ) = [ 112, 122, 124 ]
#  DERootGenera( 127 ) = [ 126 ]
#  DERootGenera( 129 ) = [ 106, 107, 127, 128 ]
#  DERootGenera( 131 ) = [ 130 ]
#  DERootGenera( 133 ) = [ 120, 123, 129, 132 ]
#  DERootGenera( 135 ) = [ 112, 119, 121, 127, 128, 130, 131, 132, 133, 134 ]
#  DERootGenera( 137 ) = [ 136 ]
#  DERootGenera( 139 ) = [ 138 ]
#
########################################################################

N := function( d, e)
   return (d * e)/GcdInt( d, e );
   end;  

G := function( d, e )
   return (d * e)/GcdInt( d, e ) - ( d + e )/ ( 2 * GcdInt( d, e ) );
   end;

FactorPairs := function( n )
   local
      factorList, prime, power, maxPower, divisorPair, 
         divisorsList, count, addDivisorsList;

   factorList := PrimePowersInt( n );
   addDivisorsList := [ ];
   divisorsList := [ [1,1] ];
   for count in [1..Length(factorList)/2] do
      prime := factorList[2*count-1];
      maxPower := factorList[2*count];        
      for power in [1..maxPower] do
         for divisorPair in divisorsList do
            Append( addDivisorsList, 
               [[ prime^power * divisorPair[1], divisorPair[2] ]] );
            Append( addDivisorsList, 
               [[ divisorPair[1], prime^power * divisorPair[2] ]] );
         od;
      od;
      Append( divisorsList, addDivisorsList );
      addDivisorsList := [ ];
   od;
   return 
      Filtered( divisorsList, divisorPair -> ( divisorPair[1] < n) and 
         ( divisorPair[2] < n ) );
   end;
   
DERootGenera := function( n )
   # negative input value signals do not write output to file
   # (OK, I admit it's an ugly hack!)
   local
      generaList, divisorPair, output;

   if IsEvenInt(n) then
      Print("n must be odd.\n");
      AppendTo(outputfile, "n must be odd.\n" );
      return;
   fi;

   # check for no output
   if n < 0 then
      output := false;
      n := -n;
   else
      output := true;
   fi;

   generaList := [ ];
   for divisorPair in FactorPairs( n ) do
      Append( generaList, [n - (divisorPair[1] + divisorPair[2])/2 ] );
   od;
   Sort( generaList );
   if output then
      AppendTo( outputfile, "DERootGenera( ", n, " ) = ",
                DuplicateFreeList( generaList ),"\n");
   fi;
   return DuplicateFreeList( generaList );
   end;

DERoots := function( g )
   local
      rootList, n, startN, endN;
   rootList := [ ];

   if IsOddInt( g ) then 
      startN := g;
   else 
      startN := g+1; 
   fi;

   if IsOddInt( Int( (6*(g+2))/5 ) ) then
      endN := Int( (6*(g+2))/5 );
   else
      endN := Int( (6*(g+2))/5 ) - 1;
   fi;

   for n in [startN,startN+2..endN] do
      if g in DERootGenera( -n ) then
         Append( rootList, [n] );
      fi;
      # -n was used to suppress writing output to file by DERootGenera()
   od;
   AppendTo( outputfile, "DERoots( ", g, " ) = ", rootList, "\n");
   return rootList;
   end;
         
DERootGeneraRange := function( n, k )
   local
      startN, endN, degree;

   if IsEvenInt(n) then
      startN := n+1;
   else
      startN := n;
   fi;

   if IsEvenInt(n+k) then
      endN := n+k-1;
   else
      endN := n+k;
   fi;

   for degree in [startN, startN+2..endN] do
      Print("DERootGenera( ",degree," ) = ",DERootGenera( degree ),"\n");
   od;
   end;

DERootsRange := function( g, k )
   local
      genus;

   for genus in [g..g+k] do
      Print("DERoots( ",genus," ) = ",DERoots( genus ),"\n");
   od;
   end;

# functions for finding Margalit-Schleimer roots

FindInverse := function( x, n )
   local
      y;
   if x < 0 then
      for y in [1..n-1] do
         if RemInt( x*y, -n ) = 1-n then
            return y;
         fi;
      od;
   else
      for y in [1..n-1] do
         if RemInt( x*y, n ) = 1 then
            return y;
         fi;
      od;
   fi;
   end;
   
PairOrderFunction := function( x, y )
   local
      a1, a2;
   a1 := x[1];
   a2 := y[1];
   if a1 < a2 then
      return true;
   else
      return false;
   fi;
   end;

SumProductSolutions := function( n )
   local
      solutionsList, aInverse, aSolution, a, b;
   
   solutionsList := [ ];
   for aInverse in [2..(n+1)/2] do;
      if GcdInt( aInverse, n ) = 1 and GcdInt( 1 - aInverse, n ) = 1 then
         a := FindInverse( aInverse, n );
         b := FindInverse( 1 - aInverse, n );
         aSolution := [a, b];
         Sort( aSolution );
         Append( solutionsList, [ aSolution ] );
      fi;
   od;
   Sort( solutionsList, PairOrderFunction );
   return solutionsList;
   end;

MSRoots := function( g )
   local
      n, aSolution, a, b, c;
   if g < 1 then
      AppendTo(outputfile, "g must be odd at least 1");
      Print("g must be odd at least 1");
      return;
   fi;

   n := 2*g + 1;
   
   if Length(SumProductSolutions( n )) = 1 then
      AppendTo(outputfile, "For g = ",g," there is one Margalit-Schleimer root:\n");
      Print("For g = ",g," there is one Margalit-Schleimer root:\n");
   else
      AppendTo(outputfile, "For g = ",g," there are ",Length(SumProductSolutions( n )),
             " Margalit-Schleimer roots:\n");
      Print("For g = ",g," there are ",Length(SumProductSolutions( n )),
             " Margalit-Schleimer roots:\n");
   fi;

   for aSolution in SumProductSolutions( n ) do
      a := aSolution[1];
      b := aSolution[2];
      c := n - a - b;
      if c < 0 then
         c := n + c;
      fi;
      AppendTo(outputfile, "( ", n, ", 0, ( ", a, ", ", b , " ); ( ",
               c , ", ", n , " ) )\n");
      Print("( ", n, ", 0, ( ", a, ", ", b , " ); ( ", c , ", ", n , " ) )\n");
   od;
   end;

NumMSRoots := function( g )
   local
      output;

   AppendTo(outputfile, "NumMSRoots( ", g, " ) = ",
            Length( SumProductSolutions( 2 * g + 1 ) ),"\n");

   return Length( SumProductSolutions( 2 * g + 1 ) );
   end;

NumMSRootsRange := function( g, k )
   local
      genus;

   for genus in [g..g+k] do
      Print("NumMSRoots( ", genus," ) = ", NumMSRoots( genus ),"\n");
   od;
   end;