;  This set of functions was used in November 1990 by Robert H. Lewis to 
; answer an unsolved question in group theory. A certain property of some 
; groups called the Gkj groups was checked by counting the number of 
; homomorphisms from Gkj -> SL(2, Z/4). It turns out that for certain k and 
; j, this number can vary from the expected number. Gilbert Baumslag 
; presented the question in this form.  Gkj is a one relator group with two 
; generators; the relation involves the two generators raised to powers k 
; and j. To find a homomorphism from Gkj into a group it is necessary and 
; sufficient to find 2 elements of that group that satisfy the same relation. 
; SL(2, Z/4) is picked because it is finite, 2x2 matrices are easy to 
; compute with, and it has an appropriate lower central series.

; Fermat note: Fermat reads the first six statements and does them. Then 
; it reads the function definitions. Then it executes the last 11 statements. 
; There is no reason to have this particular order of statements in this file. 
; In particular, the function definitions can be in any order.

&(t = 1);
&(m = 0);     ;; Turn off the mouse interrupt capability. 

 
;; Find [p] among the list of 2x2 matrices representing SL(2,Z/4).
 Function Find(i,size,j,ans) =  
    ans := 0; 
    for j = i, size do 
       if p[1] <> e[j] then &]
       fi; 
       if p[2] <> f[j] then &]
       fi; 
       if p[3] <> g[j] then &]
       fi; 
       if p[4] = h[j] then 
          ans := &_m(j); 
          & > 
       fi
    od; 
    ans.;

;; Find where each element [p]'s inverse transpose is in SL(2,Z/4).
 Function SetInvTr(i) =  
    for i = 1, invs do 
       [p] := [(e[i],f[i],g[i],h[i])]; 
       { ith element of SL(2,Z/4)  }
       [p] := Transpose[p]; 
       [p] := [(p[4],-p[2],-p[3],p[1])]; 
       { fast inverse  }
       it[i] := Find(1, &_m(invs))
    od.;

;; Insert is used by DelConj.
 Function Insert(x1,i,k) =  
    i := 1; 
    while i <= j do 
       if x1 > conjs[i] then 
          for k = j + 1, i + 1, - 1 do 
             conjs[k] := conjs[k-1]
          od; 
          conjs[i] := x1; 
          &_m(j :+ ); 
          &}
       else
          if x1 = conjs[i] then 
             &}
          else
             &_m(i :+ )
          fi
       fi
    od; 
    conjs[i] := x1; 
    &_m(j :+ ).;

; Compute the conjugacy classes in SL(2,Z/4). Look at each element of 
; SL(2,Z/4). Delete from the list of SL(2,Z/4) all of its conjugates.  
; When finished, one representative is left in each class. 
 Function DelConj(i,j,k,n,spot) =  
    Array conjs[invs\2]; 
    Array a1[2,2]; 
    Array count[2]; 
    [count] := 0; 
    &_m(size := invs); 
    i := 1; 
    while i <= size do 
       [a1] := [(e[i],f[i],g[i],h[i])]; 
       j := 0; 
       { Fill the array conjs to tell where the conjugates of [a1] are. }
       for k = 1, invs do 
          [p] := [(e1[k],f1[k],g1[k],h1[k])]; 
          [q] := [(p[4],-p[2],-p[3],p[1])]; 
          [p] := [q]*[a1]*[p]; 
          spot := Find(i, size); 
          if spot = 0 then 
             !!'spot is 0'
          fi; 
          if spot <> i then 
             Insert(spot)
          fi
       od; 
       !!('number of conjugates is ', &_m(j + 1)); 
       [count] := [count] \_ &_m(j + 1); 
       { Delete the conjugates of [a1].  Change [it] accordingly.   }
       for k = 1, j do 
          &_m(size:-); 
          for n = 1, conjs[k] - 1 do 
             if it[n] = conjs[k] then 
                it[n] := i
             fi; 
             if it[n] > conjs[k] then 
                &_m(it[n] := it[n] - 1)
             fi
          od; 
          for n = conjs[k], size do 
             e[n] := e[n+1]; 
             f[n] := f[n+1]; 
             g[n] := g[n+1]; 
             h[n] := h[n+1]; 
             it[n] := it[n+1]; 
             if it[n] = conjs[k] then 
                it[n] := i
             fi; 
             if it[n] > conjs[k] then 
                &_m(it[n]:-)
             fi
          od
       od; 
       &_m(i :+ )
    od; 
    [count] := [count[3~]]; 
    @([conjs], [a1]).;
 
;; OneDet tests every possible 2x2 matrix over Z/4 and saves the 
;; invertible ones.  j1 is owned by CreateData.
 Function OneDet(i,j,k,l) =  
    for i = 0, dim - 1 do 
       for j = 0, dim - 1 do 
          for k = 0, dim - 1 do 
             for l = 0, dim - 1 do 
                if i*l - k*j = 1 then 
                   &_m(j1 := j1 + 1); 
                   e[j1] := i; 
                   f[j1] := j; 
                   g[j1] := k; 
                   h[j1] := l
                fi
             od
          od
       od
    od.;

;; Count the number of all 2x2 matrices over Z/4 that have determinant 1.
 Function CountInv(i,j,k,l,j1) =  
    for i = 0, dim - 1 do 
       for j = 0, dim - 1 do 
          for k = 0, dim - 1 do 
             for l = 0, dim - 1 do 
                if i*l - k*j = 1 then 
                   &_m(j1 :+ )
                fi
             od
          od
       od
    od; 
    j1.;

 Function CreateData(j1) =  
    Array e[invs]; 
    Array f[invs]; 
    Array g[invs]; 
    Array h[invs]; 
    OneDet; 
    [e1] := [e]; 
    [f1] := [f]; 
    [g1] := [g]; 
    [h1] := [h].;
            
;    Main tries one pair of exponents, exp1 and exp2.  It counts the number 
; of homomorphisms from the group they define into the 2x2 matrix group.
; The basic word that defines the group has been massaged to minimize
; the amount of work within the innermost loop.  That is also the reason 
; for using the polynomials -- part of the word can be evaluated in advance 
; "once and for all".
 Function Main(exp1,exp2,i,j,k,l,yes) =  
    time := &T; 
    hits := 0; 
    for i = 1, size do 
       [c] := [(e[i],f[i],g[i],h[i])]; 
       [ci] := [c]^exp1; 
       [cj] := [c]^exp2; 
       [prod1] := [cj]*[b]; 
       [prod2] := [b]*[cj]; 
       { If [c] is in the center, so that count[i] = 1,
         the answer is easy. Increment and cycle.   }
       if count[i] = 1 then 
          &_m(hits := hits + invs); 
          &]
       fi; 
       for j = 1, invs do 
          &(m = 1); 
          { Allow mouse interrupts.   }
          [a] := [(e1[j],f1[j],g1[j],h1[j])]; 
          [p] := [ci]*[a]; 
          &(m = 0); 
          { Turn off mouse interrupts.   }
          [a1] := [(p[4],-p[2],-p[3],p[1])]; 
          [test] := [a1]*[a]*[p]; 
          [test] := [prod1] - [prod2]*[test]; 
          for k = 1, invs do 
             yes := 1; 
             for l = 1, 4 do 
                if test[l]#(h1[k], g1[k], f1[k], e1[k]) <> 0 then 
                   yes := 0; 
                   & > 
                fi
             od; 
             if yes = 1 then 
                &_m(hits := hits + count[i])
             fi
          od
       od
    od; 
    !!('time is ', &_m((&T - time)/60), 'seconds.'); 
    !!('i,  j,  hits are ', exp1:3, exp2:3, hits).;

; Driver is the 'main program'. It first counts the number of invertible 
; 2x2 matrices over Z/dim, using CountInv. (dim = 4 is the right choice, 
; for which invs is 48). CreateData then actually creates the 48 matrices. 
; Matrices are stored in 4 arrays [e], [f], [g], [h],one for the upper 
; left coordinates, one for the lower left, etc. (This is inelegant and 
; could be done better with Fermat's array of arrays.)  These are 
; duplicated into [e1], [f1], etc. Then SetInvTr finds where each matrix's 
; inverse transpose is. This data is stored in array [it]. DelConj breaks 
; the group into conjugacy classes. (This is done because the basic defining 
; word is invariant under conjugation, so we save some time this way.)  
; When it's finished, [e], [f], [g], [h] store 1 representative from each 
; class, and [it] says where the inverse transpose of each of these 
; representatives is.  But it turns out that it[j] = j for all j.  Had 
; this not been true, a further speed up would have been possible. Then a 
; loop indexed on i tests the exponents from the arrays [try1] and [try2] 
; thereby looking for maps from the group Gkj,where k = try1[i] and j = 
; try2[i]. Main counts the number of homomorphisms from the Gkj group into 
; SL(2,Z/4).  
 Function Driver(size,time,invs) =  
    time := &T; 
    Array a[2,2]; 
    Array b[2,2]; 
    Array c[2,2]; 
    Array ci[2,2]; 
    Array cj[2,2]; 
    Array p[2,2]; 
    Array a1[2,2]; 
    Array b1[2,2]; 
    Array q[2,2]; 
    Array test[2,2]; 
    [b] := [(x,z,y,w)]; 
    invs := CountInv; 
    !!('order of group is ', invs); 
    Array it[invs]; 
    CreateData; 
    SetInvTr; 
    !!('initial [it]'); 
    ![it;   !; 
    DelConj; 
    !!('time is ', &_m((&T - time)/60), ' seconds.'); 
    !!('size is ', size); 
    { size = number of conjugacy classes.   }
    !!('Now [it] is '); 
    ![it[1~size];  !; 
    !!'[count] is '; 
    ![count;  !; 
    for i = 1, 100 do 
       Main(try1[i], try2[i])
    od.;
 
?dim;                      ; Ask user for dim.
&(p = dim);                ; Go to modular mode mod dim. In converting, dim
                           ;   becomes 0.  (of course)
dim := -1;                 
&_m(dim := dim + 1);       ;  dim now holds dim, not 0.
Array try1[100];           ; Create exponents to try for k and j.
Array try2[100];
[try1[1~7]] := 1;
[try1[8~14]] := 2;
&(a = 0);
[try2] := [<i=0,99> ( &_m (i|7 + 2) )];   
&(a = 1);

&x;

 Function Degree(n,i,j,k,diff,answer) =  
    { Create matrix of degrees of entries in [x]. n = dim. of [x]. }
    Array y[n,n]; 
    for i = 1, n*n do 
       if x[i] = 0 then  y[i] := -10^10   { " -infinity" }
          else y[i] := Deg(x[i]) fi
    od; 
    answer := 0; 
    { This loop simulates row and column manipulations in [x] by working
       with the degrees in [y].  }
    for i = 1, n do 
       { Try to factor out as many powers as you can.  }
       Rowscan(n, i); 
       Colscan(n, i); 
       { Find spot with lowest degree >= 0.  }
       j := Det_ + ([y[i~,i~]] + 1); 
       if j = -1 then  & > fi; 
       { Convert linear to (row, column) coordinates in [y].  }
       k := (j - 1)\(n - i + 1) + i; 
       j := (j - 1)|(n - i + 1) + i; 
       { If the remaining matrix is all constants, stop now.  }
       if Maxdegree(n, i) = 0 then & >  fi; 
       { Simulate row and column manipulations.  }
       Switchrow([y[,i~]], i, j); 
       Switchcol([y[i~,]], i, k); 
       for j = i + 1, n do 
          diff := y[j,i] - y[i,i]; 
          for k = i + 1, n do 
             [y[j,k]] := Max(y[j,k], diff + y[i,k])
          od
       od; 
       answer := answer + y[i,i]
    od; 
    answer.;
					
 Function Rowscan(n,i,r,j,deg) =  
    for r = i, n do 
       { Find position and value of smallest nonnegative entry in rth row
              in columns >= i.  }
       j := Det_ + ([y[r,i~]] + 1); 
       deg := y[r,j+i-1]; 
       { If deg is positive, factor it out, changing answer accordingly.  }
       if deg > 0 then 
          [y[r,i~]] := [y[r,i~]] - deg; 
          answer := answer + deg
       fi
    od.;
						
 Function Colscan(n,i,r,j,deg) =  
    for r = i, n do 
       { Find position and value of smallest nonnegative entry in rth column
              in rows >= i.  }
       j := Det_ + ([y[i~,r]] + 1); 
       deg := y[j+i-1,r]; 
       { If deg is positive, factor it out, changing answer accordingly.  }
       if deg > 0 then 
          [y[i~,r]] := [y[i~,r]] - deg; 
          answer := answer + deg
       fi
    od.;
					
 Function Maxdegree(n,i,r,s,deg) =  
    deg := 0; 
    for r = i, n do 
       for s = i, n do 
          if y[r,s] > deg then 
             deg := y[r,s]; 
             &> 
          fi
       od; 
       if deg > 0 then &>  fi
    od; 
    deg.;
		 
 Function Max(a,b) =  if a > b then a else b fi.;		 

{ Assume that t is the polynomial variable. }   

 Function CharPoly(a,n,i,pm) =  
    { Assume t is a polynomial variable.  }
    n := Cols[a]; 
    Array c[n]; 
    [q] := [a]; 
    c[1] := -Trace[a]; 
    for i = 2, n do 
       [q] := [q] + c[i-1]*[1]; 
       [q] := [a]*[q]; 
       c[i] := -Trace[q]/i
    od; 
    { Correct for odd n.  }
    if n|2 then 
       pm :=  - 1
    else
       pm := 1
    fi.;

 Function Invert(a,q,n,i) =  
    { Set [q] = the inverse of [a]  }
    n := Cols[a]; 
    Array c[n]; 
    [p] := [a]; 
    c[1] := -Trace[a]; 
    for i = 2, n do 
       [q] := [p] + c[i-1]*[1]; 
       [p] := [a]*[q]; 
       c[i] := -Trace[p]/i
    od; 
    [q] := [q]/(-c[n]); 
    @([p], [c]).;