###########################################################################
 #                                                                         #
 # This file defines a number of GAP functions for computing with PLRs.    #
 # Please refer to "Primitive Lambda-Roots" for definitions.               #
 # Note one small difference - invariant factors and elementary divisors   #
 #    are ordered differently here! (This may change.)                     #
 # Note that GAP already defines such functions as:                        #
 # - Lambda(n)          (Carmichael's lambda-function)                     #
 # - Phi(n)             (Euler's phi-function)                             #
 # - MoebiusMu(n)       (Moebius' mu-function                              #
 # - OrderMod(x,n)      (The order of x mod n)                             #
 # - PowerMod(x,d,n)    (x^d mod n)                                        #
 #                                                                         #
 # The new functions are the following. The first few are calculated by    #
 #    formulae given in "Primitive Lambda-Roots".                          #
 # - Xi(n)              (Phi(n)/Lambda(n))                                 #
 # - LambdaStar(m)      (The maximum n for which Lambda(n) divides m)      #
 # - ElemDivs(n)        (The list of elementary divisors of U(n),          #
 #                       ordered by prime power and decreasing exponent)   #
 # - InvFacts(n)        (The list of invariant factors of U(n),            #
 #                       ordered by reverse divisibility)                  #
 # - NrPLRs(n)          (The number of PLRs of n)                          #
 # - NrPCs(n)           (The number of power classes of PLRs of n)         #
 # - NrNegatingPLRs(n)  (The number of negating PLRs of n)                 #
 # - NrNegatingPCs(n)                                                      #
 # - NrFraternities(n)  (The number of fraternities of PLRs of n)          #
 #                                                                         #
 # From here on we don't know formulae, so we list!                        #
 #    This is likely to be much slower.                                    #
 # - PowerClass(x,n)    (The power class of x mod n, where x is a PLR)     #
 # - EPowerClass(x,n)   (All powers of x mod n)                            #
 # - IsPLR(x,n), PLRs(n)  (The test for and list of PLRs of n)             #
 # - SmallestPLR(n)     (the smallest PLR of n)                            #
 # - IsNegatingPLR(x,n), NegatingPLRs(n)                                   #
 # - IsInwardPLR(x,n), InwardPLRs(n)                                       #
 # - NrInwardPLRs(n) and NrInwardPCs(n)                                    #
 # - IsStrongPLR(x,n), StrongPLRs(n)  (strong = inward and non-negating)   #
 # - NrStrongPLRs(n) and NrStrongPCs(n)                                    #
 # - IsSelfSeekingPLR(x,n), SelfSeekingPLRs(n)                             #
 # - SubofUnits(l,n)     (The subgroup of U(n) generated by the list l)    #
 # - ListSubofUnits(l,n), NrSubofUnits(l,n)                                #
 #                                                                         #
 ###########################################################################
 #
 # Xi(n) = Phi(n)/Lambda(n), measure of how far PLRs from primitive roots
 #
Xi:=function(n)
  return QuoInt(Phi(n),Lambda(n));
end;
 #
 # The function "LambdaStar" is computed using the algorithm
 # given in "Primitive Lambda-Roots"
 #
LambdaStar:=function(m)
  local l,x,y,p;
  if RemInt(m,2)=1 then return 2; fi;
  l:=PrimePowersInt(m);
  x:=2^(l[2]+2);
  for p in [1..QuoInt(m,2)] do
    if RemInt(m,2*p)=0 then
      if IsPrimeInt(2*p+1) then
        y:=1;
        while RemInt(m,y)=0 do y:=y*(2*p+1); od;
        x:=x*y;
      fi;
    fi;
  od;
  return x;
end;
 #
 # The function "ppfi" computes the prime power factors of its integer
 #    argument, based on the GAP function PrimePowersInt (which gives
 #    primes and exponents).
 #
ppfi:=function(n)
  local l,ll,m;
  l:=PrimePowersInt(n);
  ll:=[];
  for m in [1..QuoInt(Size(l),2)] do
    Add(ll, l[2*m-1]^l[2*m]);
  od;
  return ll;
  end;
 #
 # The function "orderpp" is used to put the result in order: we order
 #    prime powers first by the prime involved and then by the exponent.
 #
orderpp:=function(q1,q2) # both prime powers
                         # order by smallest prime then largest exponent
  local p1,p2;
  p1:=FactorsInt(q1)[1]; p2:=FactorsInt(q2)[1];
  if p1<>p2 then
    return p1<p2;
  else
    return q1>q2;
  fi;
  end;
 #
 # The function "ElemDivs" computes the elementary divisors of U(n),
 #    in the order given by the function "orderpp" above.
 # It computes the elementary divisors of U(q) for the
 #    prime power factors q of n and concatenate.
 # If q is even then ElemDivs = [], [2], or [q/4,2] according as q=2, 4, or >4.
 # If q is odd then ElemDivs is the list of prime power factors of Phi(q).
 # Finally sort the list using the ordering given by "orderpp".
 #
ElemDivs:=function(n)
  local l,ll,x;
  l:=ppfi(n);
  ll:=[];
  for x in l do
    if RemInt(x,2)=0 then
      if x>2 then Add(ll,2); fi;
      if x>4 then Add(ll,QuoInt(x,4)); fi;
    else Append(ll,ppfi(Phi(x)));
    fi;
  od;
  Sort(ll,orderpp); return ll;
  end;
 #
 # The function "InvFacts" computes the invariant factors of U(n),
 #    in decreasing order by divisibility.
 # It starts with the elementary divisors and combines them in the usual way.
 # Note that the largest invariant factor is Lambda(n).
 #
isntnull:=function(x)
  return x<>[];
end;
 #
InvFacts:=function(n)
  local l,ll,lll,m,x;
  l:=ElemDivs(n); lll:=[];
  while l<>[] do
    ll:=List(l, x->PrimePowersInt(x));
    m:=Size(l); x:=1;
    while m>=1 do
      if (m=1) or (ll[m][1]<>ll[m-1][1]) then
        x:=x*l[m]; l[m]:=[]; 
      fi;
    m:=m-1;
    od;
    Add(lll,x); l:=Filtered(l, isntnull);
  od;
  return lll;
end;
 #
 # The function "eltsoforder" finds the number of elements of order m
 #    in an abelian group where l is the list of cyclic factors in any
 #    direct sum decomposition. See "Primitive Lambda-Roots".
 #
eltsoforder:=function(m,l)
  local t,x,y,z;
  x:=0;
  for t in [1..m] do
    if RemInt(m,t)=0 then
       y:=MoebiusMu(QuoInt(m,t));
       for z in l do y:=y*GcdInt(t,z); od;
       x:=x+y;
    fi;
  od;
  return x;
end;
 #
 # Now the number of PLRs is obtained when l is the list of invariant factors
 #    and m is its first element
 #
NrPLRs:=function(n)
  local l;
  l:=InvFacts(n);
  if n<=2 then return 1; fi;
  return eltsoforder(l[1],l);
end;
 #
 # The size of a power class is Phi(Lambda(n)). Dividing by this number gives
 # the number of power classes of PLRs
 #
NrPCs:=function(n)
  return QuoInt(NrPLRs(n),Phi(Lambda(n)));
end;
 #
 # The number of negating PLRs is zero if the Sylow $2$-subgroup is not
 #    homocyclic, and is 1/(2^s-1) of the number of PLRs if it is
 #    homocyclic of rank s. See "Primitive Lambda-Roots".
 #
NrNegatingPLRs:=function(n)
  local l,s;
  if n<=2 then return 1; fi; # otherwise l[1] is even
  l:=ElemDivs(n);
  s:=0;
  while Size(l)>=s+1 and RemInt(l[s+1],2)=0 do s:=s+1; od;
  if l[1]=l[s] then return QuoInt(NrPLRs(n),2^s-1);
    else return 0;
  fi;
end;
 #
NrNegatingPCs:=function(n)
  return QuoInt(NrNegatingPLRs(n),Phi(Lambda(n)));
end;
 #
 # NrFraternities(n) returns the number of fraternities of n
 #
NrFraternities:=function(n)
  local l,s,e;
  if n<=2 then return 1; fi;
  l:=ElemDivs(n);
  s:=0;
  while Size(l)>=s+1 and RemInt(l[s+1],2)=0 do s:=s+1; od;
  if l[1]>2 then return QuoInt(NrPLRs(n), 2^s);
    else return QuoInt(NrPLRs(n), 2^s-1);
  fi;
end;
 #
 # IsPLR(x,n) checks if x is a PLR of n
 #
IsPLR:=function(x,n)
  return OrderMod(x,n)=Lambda(n);
end;
 #
 # SmallestPLR(n) does what it says
 #
SmallestPLR:=function(n)
  local x,m;
  x:=0; m:=Lambda(n);
  while OrderMod(x,n)<>m do x:=x+1; od;
  return x;
end;
 #
 # PowerClass(x,n) gives the power class of x mod n
 #
PowerClass:=function(x,n) # we assume x is a PLR
  local l,i;
  l:=[];
  for i in [1..Lambda(n)] do
    if GcdInt(i,Lambda(n))=1 then Add(l,PowerMod(x,i,n)); fi;
  od;
  return Set(l);
end;
 #
EPowerClass:=function(x,n) # all powers!
  local l,i;
  l:=[];
  for i in [1..Lambda(n)] do Add(l,PowerMod(x,i,n)); od;
  return Set(l);
end;
 #
 # Produce a list of PLRs by direct checking
 #
PLRs:=function(n)
  local x,l,m;
  l:=[]; m:=Lambda(n);
  for x in [0..n-1] do
    if OrderMod(x,n)=m then Add(l,x); fi;
  od;
  return l;
end;
 #
 # Code for listing negating PLRs
 #
IsNegatingPLR:=function(x,n)
  local e;
  e:=QuoInt(Lambda(n),2);
  return PowerMod(x,e,n)=n-1;
end;
 #
NegatingPLRs:=function(n)
  local e;
  e:=QuoInt(Lambda(n),2);
  return Filtered(PLRs(n), x->PowerMod(x,e,n)=n-1);
end;
 # 
 # Code for listing inward PLRs
 #
IsInwardPLR:=function(x,n)
  return GcdInt(x-1,n)=1;
end;
 #
InwardPLRs:=function(n)
  return Filtered(PLRs(n), x->GcdInt(x-1,n)=1);
end;
 #
NrInwardPLRs:=function(n);
  return Size(InwardPLRs(n));
end;
 #
NrInwardPCs:=function(n);
  return QuoInt(NrInwardPLRs(n),Phi(Lambda(n)));
end;
 #
 # Code for strong PLRs
 #
IsStrongPLR:=function(x,n)
  return IsInwardPLR(x,n) and not IsNegatingPLR(x,n);
end;
 #
StrongPLRs:=function(n)
  local e;
  e:=QuoInt(Lambda(n),2);
  return Filtered(InwardPLRs(n),x->PowerMod(x,e,n)<>n-1);
end;
 #
NrStrongPLRs:=function(n)
  return Size(StrongPLRs(n));
end;
 #
NrStrongPCs:=function(n)
  return QuoInt(NrStrongPLRs(n),Phi(Lambda(n)));
end;
 # 
 # Code for self-seeking PLRs
 #
IsSelfSeekingPLR:=function(x,n)
  local l;
  l:=EPowerClass(x,n);
  return (IsStrongPLR(x,n) and ((x-1 in l) or (n+1-x in l)));
end;
 #
SelfSeekingPLRs:=function(n)
  return Filtered(StrongPLRs(n),x->IsSelfSeekingPLR(x,n));
end;
 #
 # "SubofUnits" returns the subgroup of U(n) generated by l
 #
SubofUnits:=function(l,n)
  local R,U,x;
  for x in l do
    if GcdInt(x,n)<>1 then return fail; fi;
  od;
  R:=ZmodnZ(n); U:=Units(R);
  return Subgroup(U,List(l,x->ZmodnZObj(x,n)));
end;
 #
ListSubofUnits:=function(l,n)
  local ll;
  ll:=SubofUnits(l,n); if ll=fail then return fail; fi;
  return Set(List(Set(ll), x-> Int(x)));
end;
 #
NrSubofUnits:=function(l,n)
  local ll;
  ll:=SubofUnits(l,n); if ll=fail then return fail; fi;
  return Size(ll);
end;