#############################################################################
##
##  This file is part of GAP, a system for computational discrete algebra.
##  This file's authors include Thomas Breuer.
##
##  Copyright of GAP belongs to its developers, whose names are too numerous
##  to list here. Please refer to the COPYRIGHT file for details.
##
##  SPDX-License-Identifier: GPL-2.0-or-later
##
##  This file is being maintained by Thomas Breuer.
##  Please do not make any changes without consulting him.
##  (This holds also for minor changes such as the removal of whitespace or
##  the correction of typos.)
##
##  This file contains methods for cyclotomics.
##


#############################################################################
##
#V  Cyclotomics . . . . . . . . . . . . . . . . . .  field of all cyclotomics
##
BindGlobal( "Cyclotomics", Objectify( NewType(
    CollectionsFamily( CyclotomicsFamily ),
    IsField and IsAttributeStoringRep ),
    rec() ) );
SetName( Cyclotomics, "Cyclotomics" );
SetLeftActingDomain( Cyclotomics, Rationals );
SetIsFiniteDimensional( Cyclotomics, false );
SetIsFinite( Cyclotomics, false );
SetIsWholeFamily( Cyclotomics, true );
SetDegreeOverPrimeField( Cyclotomics, infinity );
SetDimension( Cyclotomics, infinity );
SetRepresentative(Cyclotomics, 0);


#############################################################################
##
#M  Conductor( <list> ) . . . . . . . . . . . . . . . . . . . . .  for a list
##
##  (This works not only for lists of cyclotomics but also for lists of
##  lists of cyclotomics etc.)
##
InstallOtherMethod( Conductor,
    "for a list",
    [ IsList ],
    function( list )
    local n, entry;
    n:= 1;
    for entry in list do
      n:= LcmInt( n, Conductor( entry ) );
    od;
    return n;
    end );


#############################################################################
##
#M  RoundCyc( <cyc> ) . . . . . . . . . . cyclotomic integer near to <cyc>
##
InstallMethod( RoundCyc, "general cyclotomic", [ IsCyclotomic ],
    function ( x )
    local n, cfs, int, i, c;
    n:= Conductor( x );
    cfs:= ExtRepOfObj( x );
    int:= EmptyPlist( n );
    for i in [ 1 .. n ] do
      c:= cfs[i];
      if IsInt( c ) then
        int[i]:= c;
      elif c < 0 then
        int[i]:= Int( c - 1/2 );
      else
        int[i]:= Int( c + 1/2 );
      fi;
    od;
    return CycList( int );
end );


#############################################################################
##
#M  RoundCycDown( <cyc> ) . . . . . . . . .  cyclotomic integer near to <cyc>
##                                                       rounding halves down
##
InstallMethod( RoundCycDown, "general cyclotomic", [ IsCyclotomic ],
    x -> CycList( List( ExtRepOfObj( x ), RoundCycDown ) ) );


#############################################################################
##
#M  ComplexConjugate( <cyc> )
#M  ComplexConjugate( <list> )
##
InstallMethod( ComplexConjugate,
    "for a cyclotomic",
    [ IsCyc ],
    cyc -> GaloisCyc( cyc, -1 ) );

InstallMethod( ComplexConjugate,
    "for a list",
    [ IsList ],
    function( list )
    local result, i;

    result:= [];
    for i in [ 1 .. Length( list ) ] do
      if IsBound( list[i] ) then
        result[i]:= ComplexConjugate( list[i] );
      fi;
    od;
    return result;
    end );


#############################################################################
##
#M  RealPart( <z> )
#M  RealPart( <list> )
##
InstallMethod( RealPart,
    "for a scalar",
    [ IsScalar ],
    z -> ( z + ComplexConjugate( z ) ) / 2 );

InstallMethod( RealPart,
    "for a list",
    [ IsList ],
    function( list )
    local result, i;

    result:= [];
    for i in [ 1 .. Length( list ) ] do
      if IsBound( list[i] ) then
        result[i]:= RealPart( list[i] );
      fi;
    od;
    return result;
    end );


#############################################################################
##
#M  ImaginaryPart( <z> )
#M  ImaginaryPart( <list> )
##
InstallMethod( ImaginaryPart,
    "for a cyclotomic",
    [ IsCyc ],
    z -> E(4) * ( ComplexConjugate( z ) - z ) / 2 );

InstallMethod( ImaginaryPart,
    "for a list",
    [ IsList ],
    function( list )
    local result, i;

    result:= [];
    for i in [ 1 .. Length( list ) ] do
      if IsBound( list[i] ) then
        result[i]:= ImaginaryPart( list[i] );
      fi;
    od;
    return result;
    end );


#############################################################################
##
#M  ExtRepOfObj( <cyc> )
##
##  <#GAPDoc Label="ExtRepOfObj:cyclotomics">
##  <ManSection>
##  <Meth Name="ExtRepOfObj" Arg='cyc' Label="for a cyclotomic"/>
##
##  <Description>
##  The external representation of a cyclotomic <A>cyc</A> with conductor
##  <M>n</M> (see <Ref Attr="Conductor" Label="for a cyclotomic"/> is
##  the list returned by <Ref Func="CoeffsCyc"/>,
##  called with <A>cyc</A> and <M>n</M>.
##  <P/>
##  <Example><![CDATA[
##  gap> ExtRepOfObj( E(5) ); CoeffsCyc( E(5), 5 );
##  [ 0, 1, 0, 0, 0 ]
##  [ 0, 1, 0, 0, 0 ]
##  gap> CoeffsCyc( E(5), 15 );
##  [ 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, -1, 0 ]
##  ]]></Example>
##  </Description>
##  </ManSection>
##  <#/GAPDoc>
##
InstallMethod( ExtRepOfObj,
    "for an internal cyclotomic",
    [ IsCyc and IsInternalRep ],
    COEFFS_CYC );


#############################################################################
##
#F  CoeffsCyc( <z>, <N> )
##
InstallGlobalFunction( CoeffsCyc, function( z, N )
    local coeffs,      # coefficients list (also intermediately)
          n,           # length of `coeffs'
          quo,         # factor by that we have to blow up
          s,           # denominator of `quo' (we have to reduce)
          factors,     #
          second,      # product of primes in second reduction step
          first,
          val,
          kk,
          pos,
          j,
          k,
          p,
          nn,
          newcoeffs;

    if IsCyc( z ) then

      # `z' is an internal cyclotomic, and therefore it is represented
      # in the smallest possible cyclotomic field.

      coeffs:= ExtRepOfObj( z ); # returns `CoeffsCyc( z, Conductor( z ) )'
      n:= Length( coeffs );
      quo:= N / n;
      if not IsInt( quo ) then
        return fail;
      fi;

    else

      # `z' is already a coefficients list
      coeffs:= z;
      n:= Length( coeffs );
      quo:= N / n;
      if not IsInt( quo ) then

        # Maybe `z' was not represented with respect to
        # the smallest possible cyclotomic field.
        # Try to reduce `z' until the denominator $s$ disappears.

        s:= DenominatorRat( quo );

        # step 1.
        # First we get rid of
        # - the $p$-parts for those primes $p$ that divide both $s$
        #   and $n / s$, and hence remain in the reduced expression,
        # - all primes but the squarefree part of the rest.

        factors:= PrimeDivisors( s );
        second:= 1;
        for p in factors do
          if s mod p <> 0 or ( n / s ) mod p <> 0 then
            second:= second * p;
          fi;
        od;
        if second mod 2 = 0 then
          second:= second / 2;
        fi;
        first:= s / second;
        if first > 1 then
          newcoeffs:= ListWithIdenticalEntries( n / first, 0 );
          for k in [ 1 .. n ] do
            if coeffs[k] <> 0 then
              pos:= ( k - 1 ) / first + 1;
              if not IsInt( pos ) then
                return fail;
              fi;
              newcoeffs[ pos ]:= coeffs[k];
            fi;
          od;
          n:= n / first;
          coeffs:= newcoeffs;
        fi;

        # step 2.
        # Now we process those primes that shall disappear
        # in the reduced form.
        # Note that $p-1$ of the coefficients congruent modulo $n/p$
        # must be equal, and the negative of this value is put at the
        # position of the $p$-th element of this congruence class.
        if second > 1 then
          for p in Factors(Integers, second ) do
            nn:= n / p;
            newcoeffs:= ListWithIdenticalEntries( nn, 0 );
            for k in [ 1 .. n ] do
              pos:= 0;
              if coeffs[k] <> 0 then
                val:= coeffs[k];
                for j in [ 1 .. p-1 ] do
                  kk:= ( ( k - 1 + j*nn ) mod n ) + 1;
                  if coeffs[ kk ] = val then
                    coeffs[ kk ]:= 0;
                  elif pos <> 0 or coeffs[kk] <> 0 or kk mod p <> 1 then
                    return fail;
                  else
                    pos:= ( kk - 1 ) / p + 1;
                  fi;
                od;
                newcoeffs[ pos ]:= - val;
              fi;
            od;
            n:= nn;
            coeffs:= newcoeffs;
          od;
        fi;
        quo:= NumeratorRat( quo );

      fi;

    fi;

    # If necessary then blow up the representation in two steps.

    # step 1.
    # For each prime `p' not dividing `n' we replace
    # `E(n)^k' by $ - \sum_{j=1}^{p-1} `E(n*p)^(p*k+j*n)'$.
    if quo <> 1 then

      for p in PrimeDivisors( quo ) do
        if p <> 2 and n mod p <> 0 then
          nn  := n * p;
          quo := quo / p;
          newcoeffs:= ListWithIdenticalEntries( nn, 0 );
          for k in [ 1 .. n ] do
            if coeffs[k] <> 0 then
              for j in [ 1 .. p-1 ] do
                newcoeffs[ ( ( (k-1)*p + j*n  ) mod nn ) + 1 ]:= -coeffs[k];
              od;
            fi;
          od;
          coeffs:= newcoeffs;
          n:= nn;
        fi;
      od;

    fi;

    # step 2.
    # For the remaining divisors of `quo' we have
    # `E(n*p)^(k*p)' in the basis for each basis element `E(n)^k'.
    if quo <> 1 then

      n:= Length( coeffs );
      newcoeffs:= ListWithIdenticalEntries( quo*n, 0 );
      for k in [ 1 .. Length( coeffs ) ] do
        if coeffs[k] <> 0 then
          newcoeffs[ (k-1)*quo + 1 ]:= coeffs[k];
        fi;
      od;
      coeffs:= newcoeffs;

    fi;

    # Return the coefficients list.
    return coeffs;
end );


#############################################################################
##
#F  IsGaussInt(<x>) . . . . . . . . . test if an object is a Gaussian integer
##
InstallGlobalFunction( IsGaussInt,
    x -> IsCycInt( x ) and (Conductor( x ) = 1 or Conductor( x ) = 4) );


#############################################################################
##
#F  IsGaussRat( <x> ) . . . . . . .  test if an object is a Gaussian rational
##
InstallGlobalFunction( IsGaussRat,
    x -> IsCyc( x ) and (Conductor( x ) = 1 or Conductor( x ) = 4) );


##############################################################################
##
#F  DescriptionOfRootOfUnity( <root> )
##
##  Let $\zeta$ denote the primitive $n$-th root of unity $`E'(n)$,
##  and suppose that we know the coefficients of $\pm\zeta^i$ w.r.t.
##  the $n$-th Zumbroich basis (see~"ZumbroichBasis").
##  These coefficients are either all $1$ or all $-1$.
##  More precisely, they arise in the base conversion from (formally)
##  successively multiplying $\pm\zeta^i$ by
##  $1 = - \sum_{j=1}^{p-1} \zeta^{jn/p}$,
##  for suitable prime diviors $p$ of $n$,
##  and then treating the summands $\pm\zeta^{i + jn/p}$ in the same way
##  until roots in the basis are reached.
##  It should be noted that all roots obtained this way are distinct.
##
##  Suppose the above procedure must be applied for the primes
##  $p_1$, $p_2$, \ldots, $p_r$.
##  Then $\zeta^i$ is equal to
##  $(-1)^r \sum_{j_1=1}^{p_1-1} \cdots \sum_{j_r=1}^{p_r-1}
##  \zeta^{i + \sum_{k=1}^r j_k n/p_k}$.
##  The number of summands is $m = \prod_{k=1}^r (p_k-1)$.
##  Since these roots are all distinct, we can compute the sum $s$ of their
##  exponents modulo $n$ from the known coefficients of $\zeta^i$,
##  and we get $s \equiv m ( i + r n/2 ) \pmod{n}$.
##  Either $m = 1$ or $m$ is even,
##  hence this congruence determines $\zeta^i$ at most up to its sign.
##
##  Now suppose that $g = \gcd( m, n )$ is nontrivial.
##  Then $i$ is determined only modulo $n/g$, and we have to decide which
##  of the $g$ possible values $i$ is.
##  This could be done by computing the values $j_{0,p}$ for one
##  candidate $i$, where $p$ runs over the prime divisors $p$ of $n$
##  for which $m$ is the product of $(p-1)$.
##
##  (Recall that each $n$-th root of unity is of the form
##  $\prod_{p\in P} \prod_{k_p=1}^{\nu_p-1} \zeta^{j_{k,p}n/p^{k_p+1}}$,
##  with $j_{0,p} \in \{ 0, 1, \ldots, p-1 \}$, $j_{k,2} \in \{ 0, 1 \}$
##  for $k > 0$, and $j_{k,p} \in \{ -(p-1)/2, \ldots, (p-1)/2 \}$ for
##  $p > 2$.
##  The root is in the $n$-th Zumbroich basis if and only if $j_{0,2} = 0$
##  and $j_{0,p} \not= 0$ for $p > 2$.)
##
##  But note that we do not have an easy access to the decomposition
##  of $m$ into factors $(p-1)$,
##  although in fact this decomposition is unique for $n \leq 65000$.
##
##  So the exponent is identified by dividing off one of the candidates
##  and then identifying the quotient, which is a $g$-th or $2 g$-th
##  root of unity.
##  (Note that $g$ is small compared to $n$.
##  $m$ divides the product of $(p-1)$ for prime divisors $p$
##  that divide $n$ at least with exponent $2$.
##  The maximum of $g$ for $n \leq 65000$ is $40$, so $2$ steps suffice
##  for these values of $n$.)
##
InstallGlobalFunction( DescriptionOfRootOfUnity, function( root )
    local coeffs,   # Zumbroich basis coefficients of `root'
          n,        # conductor of `n'
          sum,      # sum of exponents with nonzero coefficients
          num,      # number of nonzero coefficients
          i,        # loop variable
          val,      # one coefficient
          coeff,    # one nonzero coefficient (either `1' or `-1')
          exp,      # candidate for the exponent
          g,        # `Gcd( n, num )'
          groot;    # root in recursion

    # Handle the trivial cases that `root' is an integer.
    if root = 1 then
      return [ 1, 1 ];
    elif root = -1 then
      return [ 2, 1 ];
    fi;
    Assert( 1, not IsRat( root ) );

    # Compute the Zumbroich basis coefficients,
    # and number and sum of exponents with nonzero coefficient (mod `n').
    coeffs:= ExtRepOfObj( root );
    n:= Length( coeffs );
    sum:= 0;
    num:= 0;
    for i in [ 1 .. n ] do
      val:= coeffs[i];
      if val <> 0 then
        sum:= sum + i;
        num:= num + 1;
        coeff:= val;
      fi;
    od;
    sum:= sum - num;

    # `num' is equal to `1' if and only if `root' or its negative
    # belongs to the basis.
    # (The coefficient is equal to `-1' if and only if either
    # `n' is a multiple of $4$ or
    # `n' is odd and `root' is a primitive $2 `n'$-th root of unity.)
    if num = 1 then
      if coeff < 0 then
        if n mod 2 = 0 then
          sum:= sum + n/2;
        else
          sum:= 2*sum + n;
          n:= 2*n;
        fi;
      fi;
      Assert( 1, root = E(n)^sum );
      return [ n, sum mod n ];
    fi;

    # Let $N$ be `n' if `n' is even, and equal to $2 `n'$ otherwise.
    # The exponent is determined modulo $N / \gcd( N, `num' )$.
    g:= GcdInt( n, num );
    if g = 1 then

      # If `n' and `num' are coprime then `root' is identified up to its sign.
      exp:= ( sum / num ) mod n;
      if root <> E(n)^exp then
        exp:= 2*exp + n;
        n:= 2*n;
      fi;

    elif g = 2 then

      # `n' is even, and again `root' is determined up to its sign.
      exp:= ( sum / num ) mod ( n / 2 );
      if root <> E(n)^exp then
        exp:= exp + n / 2;
      fi;

    else

      # Divide off one of the candidates.
      # The quotient is a `g'-th or $2 `g'$-th root of unity,
      # which can be identified by recursion.
      exp:= ( sum / num ) mod ( n / g );
      groot:= DescriptionOfRootOfUnity( root * E(n)^(n-exp) );
      if n mod 2 = 1 and groot[1] mod 2 = 0 then
        exp:= 2*exp;
        n:= 2*n;
      fi;
      exp:= exp + groot[2] * n / groot[1];

    fi;

    # Return the result.
    Assert( 1, root = E(n)^exp );
    return [ n, exp mod n ];
end );


#############################################################################
##
#F  Atlas1( <n>, <i> )  . . . . . . . . . . . . . . . utility for EB, ..., EH
##
##  is the value $\frac{1}{i} \sum{j=1}^{n-1} z_n^{j^i}$
##  for $2 \leq i \leq 8$ and $<n> \equiv 1 \pmod{i}$;
##  if $i > 2$, <n> should be a prime to get sure that the result is
##  well-defined;
##  `Atlas1' returns the value given above if it is a cyclotomic integer.
##  (see: Conway et al, ATLAS of finite groups, Oxford University Press 1985,
##        Chapter 7, Section 10)
##
BindGlobal( "Atlas1", function( n, i )
    local atlas, k, pos;

    if not IsInt( n ) or n < 1 then
      Error( "usage: EB(<n>), EC(<n>), ..., EH(<n>) with appropriate ",
             "integer <n>" );
    elif n mod i <> 1 then
      Error( "<n> not congruent 1 mod ", i );
    fi;

    if n = 1 then
      return 0;
    fi;

    atlas:= [ 1 .. n ] * 0;

    if i mod 2 = 0 then
      # Note that `n' is odd in this case.
      for k in [ 1 .. QuoInt( n-1, 2 ) ] do
        # summand `2 * E(n)^(k^i)'
        pos:= ( (k^i) mod n ) + 1;
        atlas[ pos ]:= atlas[ pos ] + 2;
      od;
    else
      for k in [ 1 .. QuoInt( n-1, 2 ) ] do
        # summand `E(n)^(k^i) + E(n)^(n - k^i)'
        pos:= ( (k^i) mod n ) + 1;
        atlas[ pos ]:= atlas[ pos ] + 1;
        pos:= ( n - pos + 1 ) mod n + 1;
        atlas[ pos ]:= atlas[ pos ] + 1;
      od;
      if n mod 2 = 0 then
        # summand `E(n)^( (n/2)^i )'
        pos:= ( ( (n/2)^i ) mod n ) + 1;
        atlas[ pos ]:= atlas[ pos ] + 1;
      fi;
    fi;
    atlas:= CycList( atlas / i );
    if not IsCycInt( atlas ) then
      Error( "result divided by ", i, " is not a cyclotomic integer" );
    fi;
    return atlas;
end );


#############################################################################
##
#F  EB( <n> ), EC( <n> ), \ldots, EH( <n> ) . . .  some ATLAS irrationalities
##
InstallGlobalFunction( EB, n -> Atlas1( n, 2 ) );
InstallGlobalFunction( EC, n -> Atlas1( n, 3 ) );
InstallGlobalFunction( ED, n -> Atlas1( n, 4 ) );
InstallGlobalFunction( EE, n -> Atlas1( n, 5 ) );
InstallGlobalFunction( EF, n -> Atlas1( n, 6 ) );
InstallGlobalFunction( EG, n -> Atlas1( n, 7 ) );
InstallGlobalFunction( EH, n -> Atlas1( n, 8 ) );


#############################################################################
##
#F  NK( <n>, <k>, <d> ) . . . . . . . . . . utility for ATLAS irrationalities
##
##  Find the (<d>+1)-st automorphism *nk of order <k> on primitive <n>-th
##  roots of unity.
##
##  Optimization steps:
##  - Since <k> is very small, `PowerModInt' should be avoided, compared to
##    powering and then reducing mod <n>.
##  - Do not call `Phi' but look at the prime factors of <n> directly.
##  - Use multiplication instead of powering for computing squares.
##
InstallGlobalFunction( NK, function( n, k, deriv )
    local nk,
          nkn,   # nk mod n
          n1,    # n-1
          pow, pow2, pow3;

    if n <= 2 then
      return fail;
    fi;

    nk:= 1;
    nkn:= 1;
    n1:= n-1;
    if k = 2 then
      # We have always `Phi( n ) mod k = 0'.
      while true do
        if nkn <> 1 then
          # `*nk' has order larger than 1.
          pow:= (nkn * nkn) mod n;
          if pow = 1 then
            # `*nk' has order 2.
            if deriv = 0 then return nk; fi;
            deriv:= deriv - 1;
            if nkn <> n1 then
              # `**nk' has order larger than 1 and thus equal to 2.
              if deriv = 0 then return -nk; fi;
              deriv:= deriv - 1;
            fi;
          fi;
        elif nkn <> n1 then
          # `**nk' has order larger than 1.
          pow:= (nkn * nkn) mod n;
          if pow = 1 then
            # `**nk' has order dividing 2.
            if deriv = 0 then return -nk; fi;
            deriv:= deriv - 1;
          fi;
        fi;
        nk:= nk + 1;
        nkn:= nk mod n;
      od;
    elif k in [ 3, 5, 7 ] then   # for odd primes
      if ( n mod ( k*k ) <> 0 ) and
         ForAll( PrimeDivisors( n ), p -> (p-1) mod k <> 0 ) then
        return fail;
      fi;
      while true do
        if nkn <> 1 then
          # `*nk' has order larger than 1.
          pow:= (nkn ^ k) mod n;
          if pow = 1 then
            # `*nk' has order dividing `k'.
            if deriv = 0 then return nk; fi;
            deriv:= deriv - 1;
          fi;
          if nkn <> n1 and pow = n1 then
            # `**nk' has order larger than 1 and dividing `k'.
            if deriv = 0 then return -nk; fi;
            deriv:= deriv - 1;
          fi;
        elif nkn <> n1 then
          # `**nk' has order larger than 1.
          pow:= (nkn ^ k) mod n;
          if pow = n1 then
            # `**nk' has order dividing `k'.
            if deriv = 0 then return -nk; fi;
            deriv:= deriv - 1;
          fi;
        fi;
        nk:= nk + 1;
        nkn:= nk mod n;
      od;
    elif k = 4 then
      # An automorphism of order 4 exists if 4 divides $p-1$ for an odd
      # prime divisor $p$ of `n', or if 16 divides `n'.
      if ForAll( PrimeDivisors( n ), p -> (p-1) mod k <> 0 )
         and n mod 16 <> 0 then
        return fail;
      fi;
      while true do
        if nkn <> 1 and nkn <> n1 then
          # `*nk' and `**nk' have order at least 2.
          pow:= (nkn * nkn) mod n;
          if pow <> 1 and ( pow * pow ) mod n = 1 then
            # `*nk' (and thus also `**nk') has order 4.
            if deriv = 0 then
              return nk;
            elif deriv = 1 then
              return -nk;
            fi;
            deriv:= deriv - 2;
          fi;
        fi;
        nk:= nk + 1;
        nkn:= nk mod n;
      od;
    elif k = 6 then
      # An automorphism of order 6 exists if automorphisms of the orders
      # 2 and 3 exist; the former is always true.
      if ( n mod 9 <> 0 ) and
         ForAll( PrimeDivisors( n ), p -> (p-1) mod 3 <> 0 ) then
        return fail;
      fi;
      while true do
        if nkn <> 1 and nkn <> n1 then
          # `*nk' and `**nk' have order at least 2.
          pow2:= (nkn * nkn) mod n;
          if pow2 <> 1 then
            # `*nk' (and thus `**nk') has order larger than 2.
            pow:= (pow2 ^ 3) mod n;
            if pow = 1 then
              # `*nk' (and thus `**nk') has order dividing 6.
              pow3:= (nkn * pow2) mod n;
              if pow3 <> 1 then
                if deriv = 0 then return nk; fi;
                deriv:= deriv - 1;
              fi;
              if pow3 <> n1 then
                if deriv = 0 then return -nk; fi;
                deriv:= deriv - 1;
              fi;
            fi;
          fi;
        fi;
        nk:= nk + 1;
        nkn:= nk mod n;
      od;
    elif k = 8 then
      # An automorphism of order 8 exists if 8 divides $p-1$ for an odd
      # prime divisor $p$ of `n', or if 32 divides `n'.
      if ForAll( PrimeDivisors( n ), p -> (p-1) mod k <> 0 )
         and n mod 32 <> 0 then
        return fail;
      fi;
      while true do
        if nkn <> 1 and nkn <> n1 then
          # `*nk' and `**nk' have order at least 2.
          pow:= (nkn ^ 4) mod n;
          if pow <> 1 and ( pow * pow ) mod n = 1 then
            # `*nk' (and thus also `**nk') has order 8.
            if deriv = 0 then
              return nk;
            elif deriv = 1 then
              return -nk;
            fi;
            deriv:= deriv - 2;
          fi;

        fi;
        nk:= nk + 1;
        nkn:= nk mod n;
      od;
    fi;

    # We have `k < 2' or `k > 8'.
    return fail;
end );


#############################################################################
##
#F  Atlas2( <n>, <k>, <deriv> ) . . . . . . utility for ATLAS irrationalities
##
BindGlobal( "Atlas2", function( n, k, deriv )
    local nk, result, i, pos;

    if not ( IsInt( n ) and IsInt( k ) and IsInt( deriv ) ) then
      Error( "usage: ATLAS irrationalities require integral arguments" );
    fi;

    nk:= NK( n, k, deriv );
    if nk = fail then
      return fail;
    fi;

    result:= [ 1 .. n ] * 0;
    for i in [ 0 .. k-1 ] do
      # summand `E(n)^(nk^i);
      pos:= ( (nk^i) mod n ) + 1;
      result[ pos ]:= result[ pos ] + 1;
    od;
    return CycList( result );
end );


#############################################################################
##
#F  EY(<n>), EY(<n>,<deriv>) . . . . . . .  ATLAS irrationalities $y_n$ resp.
#F                                          $y_n^{<deriv>}$
#F  ... ES(<n>), ES(<n>,<deriv>)              ... $s_n$ resp. $s_n^{<deriv>}$
##
InstallGlobalFunction( EY, function(arg)
    if   Length(arg)=1 then return Atlas2(arg[1],2,0);
    elif Length(arg)=2 then return Atlas2(arg[1],2,arg[2]);
    else Error( "usage: EY(n) resp. EY(n,deriv)" ); fi;
end );

InstallGlobalFunction( EX, function(arg)
    if   Length(arg)=1 then return Atlas2(arg[1],3,0);
    elif Length(arg)=2 then return Atlas2(arg[1],3,arg[2]);
    else Error( "usage: EX(n) resp. EX(n,deriv)" ); fi;
end );

InstallGlobalFunction( EW, function(arg)
    if   Length(arg)=1 then return Atlas2(arg[1],4,0);
    elif Length(arg)=2 then return Atlas2(arg[1],4,arg[2]);
    else Error( "usage: EW(n) resp. EW(n,deriv)" ); fi;
end );

InstallGlobalFunction( EV, function(arg)
    if   Length(arg)=1 then return Atlas2(arg[1],5,0);
    elif Length(arg)=2 then return Atlas2(arg[1],5,arg[2]);
    else Error( "usage: EV(n) resp. EV(n,deriv)" ); fi;
end );

InstallGlobalFunction( EU, function(arg)
    if   Length(arg)=1 then return Atlas2(arg[1],6,0);
    elif Length(arg)=2 then return Atlas2(arg[1],6,arg[2]);
    else Error( "usage: EU(n) resp. EU(n,deriv)" ); fi;
end );

InstallGlobalFunction( ET, function(arg)
    if   Length(arg)=1 then return Atlas2(arg[1],7,0);
    elif Length(arg)=2 then return Atlas2(arg[1],7,arg[2]);
    else Error( "usage: ET(n) resp. ET(n,deriv)" ); fi;
end );

InstallGlobalFunction( ES, function(arg)
    if   Length(arg)=1 then return Atlas2(arg[1],8,0);
    elif Length(arg)=2 then return Atlas2(arg[1],8,arg[2]);
    else Error( "usage: ES(n) resp. ES(n,deriv)" ); fi;
end );


#############################################################################
##
#F  Atlas3( <n>, <k>, <deriv> ) . . . . . . utility for ATLAS irrationalities
##
BindGlobal( "Atlas3", function( n, k, deriv )
    local nk, l, pow, val, i;

    nk:= NK( n, k, deriv );
    if nk = fail then
      return fail;
    fi;
    l:= 0 * [ 1 .. n ];
    pow:= 1;
    val:= 1;
    for i in [ 1 .. k ] do
      l[ pow + 1 ]:= val;
      pow:= (nk * pow) mod n;
      val:= -val;
    od;

    return CycList( l );
end );


#############################################################################
##
#F  EM( <n>[, <deriv>] )  . . . . . . . . ATLAS irrationality $m_{<n>}$ resp.
##                                                        $m_{<n>}^{<deriv>}$
InstallGlobalFunction( EM, function( arg )
    local n;

    n:= arg[1];
    if 2 < Length( arg ) or not IsInt( n ) or n < 3 then
      Error( "usage: EM( <n>[, <deriv>] )" );
    elif Length( arg ) = 1 then
      return Atlas3( n, 2, 0 );
    else
      return Atlas3( n, 2, arg[2] );
    fi;
end );


#############################################################################
##
#F  EL( <n>[, <deriv>] )  . . . . . . . . ATLAS irrationality $l_{<n>}$ resp.
##                                                        $l_{<n>}^{<deriv>}$
InstallGlobalFunction( EL, function( arg )
    local n;

    n:= arg[1];
    if 2 < Length( arg ) or not IsInt( n ) or n < 3 then
      Error( "usage: EL( <n>[, <deriv>] )" );
    elif Length( arg ) = 1 then
      return Atlas3( n, 4, 0 );
    else
      return Atlas3( n, 4, arg[2] );
    fi;
end );


#############################################################################
##
#F  EK( <n>[, <deriv>] )  . . . . . . . . ATLAS irrationality $k_{<n>}$ resp.
##                                                        $k_{<n>}^{<deriv>}$
InstallGlobalFunction( EK, function( arg )
    local n;

    n:= arg[1];
    if 2 < Length( arg ) or not IsInt( n ) or n < 3 then
      Error( "usage: EK( <n>[, <deriv>] )" );
    elif Length( arg ) = 1 then
      return Atlas3( n, 6, 0 );
    else
      return Atlas3( n, 6, arg[2] );
    fi;
end );


#############################################################################
##
#F  EJ( <n>[, <deriv>] )  . . . . . . . . ATLAS irrationality $j_{<n>}$ resp.
##                                                        $j_{<n>}^{<deriv>}$
InstallGlobalFunction( EJ, function( arg )
    local n;

    n:= arg[1];
    if 2 < Length( arg ) or not IsInt( n ) or n < 3 then
      Error( "usage: EJ( <n>[, <deriv>] )" );
    elif Length( arg ) = 1 then
      return Atlas3( n, 8, 0 );
    else
      return Atlas3( n, 8, arg[2] );
    fi;
end );


#############################################################################
##
#F  ER( <n> ) . . . . ATLAS irrationality $r_{<n>}$ (pos. square root of <n>)
##
##  Different from other {\ATLAS} irrationalities,
##  `ER' (and its synonym `Sqrt') can be applied also to noninteger
##  rationals.
##
InstallGlobalFunction( ER, function( n )
    local factor, factors, pair;

    if IsInt( n ) then

      if n = 0 then
        return 0;
      elif n < 0 then
        factor:= E(4);
        n:= -n;
      else
        factor:= 1;
      fi;

      # Split `n' into the product of a square and a squarefree number.
      factors:= Collected( Factors( n ) );
      n:= 1;
      for pair in factors do
        if pair[2] mod 2 = 0 then
          factor:= factor * pair[1]^(pair[2]/2);
        else
          n:= n * pair[1];
          if pair[2] <> 1 then
            factor:= factor * pair[1]^((pair[2]-1)/2);
          fi;
        fi;
      od;
      if   n mod 4 = 1 then
        return factor * ( 2 * EB(n) + 1 );
      elif n mod 4 = 2 then
        return factor * ( E(8) - E(8)^3 ) * ER( n / 2 );
      elif n mod 4 = 3 then
        return factor * (-E(4)) * ( 2 * EB(n) + 1 );
      fi;

    elif IsRat( n ) then
      factor:= DenominatorRat( n );
      return ER( NumeratorRat( n ) * factor ) / factor;
    else
      Error( "argument must be rational" );
    fi;
end );


#############################################################################
##
#F  EI( <n> ) . . . . ATLAS irrationality $i_{<n>}$ (the square root of -<n>)
##
InstallGlobalFunction( EI, n -> E(4) * ER(n) );


#############################################################################
##
#M  Sqrt( <rat> ) . . . . . . . . . . . . . . . . .  square root of rationals
##
InstallMethod( Sqrt,
    "for a rational",
    [ IsRat ],
    ER );


#############################################################################
##
#F  StarCyc( <cyc> )  . . . . the unique nontrivial Galois conjugate of <cyc>
##
InstallGlobalFunction( StarCyc, function( cyc )
    local n, gens, cand, exp, img;

    n:= Conductor( cyc );
    if n = 1 then
      return fail;
    fi;
    gens:= Flat( GeneratorsPrimeResidues( n ).generators );
    cand:= fail;
    for exp in gens do
      img:= GaloisCyc( cyc, exp );
      if img <> cyc then
        if cand = fail then
          cand:= img;
        elif cand <> img then
          return fail;
        fi;
      fi;
    od;
    for exp in gens do
      img:= GaloisCyc( cand, exp );
      if img <> cyc and img <> cand then
        return fail;
      fi;
    od;
    return cand;
    end );


#############################################################################
##
#F  AtlasIrrationality( <irratname> )
##
InstallGlobalFunction( AtlasIrrationality, function( irratname )
    local len, pos, sign, pos2, coeff, letts, funcs, recurse, lpos, N,
          dashes, irrat, qN, gal, oldpos;

    # Check the argument.
    if not IsString( irratname ) or IsEmpty( irratname ) then
      return fail;
    fi;
    len:= Length( irratname );

    # Get the first sign.
    pos:= 1;
    sign:= 1;
    if irratname[1] = '-' then
      sign:= -1;
      pos:= 2;
    elif irratname[1] = '+' then
      pos:= 2;
    fi;

    # Get the first coefficient.
    if pos <= len and IsDigitChar( irratname[ pos ] ) then
      pos2:= pos;
      while pos2 <= len and IsDigitChar( irratname[ pos2 ] ) do
        pos2:= pos2+1;
      od;
      coeff:= sign * Int( irratname{ [ pos .. pos2-1 ] } );
      pos:= pos2;
    else
      coeff:= sign;
    fi;
    if len < pos then
      return coeff;
    fi;

    # Get the first atomic irrationality (with dashes).
    letts:= "bcdefghijklmrstuvwxyz";
    funcs:= List( "BCDEFGHIJKLMRSTUVWXY", x -> ValueGlobal( [ 'E', x ] ) );
    Add( funcs, E );

    # Is the coefficient an integer summand?
    if irratname[ pos ] in "+-" then
      recurse:= AtlasIrrationality( irratname{ [ pos ..
                    Length( irratname ) ] } );
      if recurse = fail then
        return fail;
      else
        return coeff + recurse;
      fi;
    fi;

    lpos:= Position( letts, irratname[ pos ] );
    if lpos = fail then
      return fail;
    fi;
    pos:= pos + 1;
    dashes:= 0;
    while pos <= len and irratname[ pos ] in "\'\"" do
      dashes:= dashes + 1;
      if irratname[ pos ] = '\"' then
        dashes:= dashes + 1;
      fi;
      pos:= pos + 1;
    od;
    pos2:= pos;
    while pos2 <= len and IsDigitChar( irratname[ pos2 ] ) do
      pos2:= pos2+1;
    od;
    if pos2 = pos and lpos = 8 then
      N:= 1;
    else
      N:= Int( irratname{ [ pos .. pos2 - 1 ] } );
    fi;
    if dashes = 0 then
      qN:= funcs[ lpos ]( N );
    else
      qN:= funcs[ lpos ]( N, dashes );
    fi;
    pos:= pos2;
    irrat:= coeff * qN;
    if len < pos then
      return irrat;
    fi;

    # Get the Galois automorphism.
    if irratname[ pos ] = '*' then
      pos:= pos + 1;
      if pos <= len and irratname[ pos ] = '*' then
        pos:= pos + 1;
        pos2:= pos;
        while pos2 <= len and IsDigitChar( irratname[ pos2 ] ) do
          pos2:= pos2 + 1;
        od;
        gal:= Int( irratname{ [ pos .. pos2-1 ] } );
        if gal = 0 then
          irrat:= ComplexConjugate( irrat );
        else
          irrat:= GaloisCyc( irrat, -gal );
        fi;
        pos:= pos2;
      elif len < pos or irratname[ pos ] in "+-&" then
        irrat:= StarCyc( irrat );
      else
        pos2:= pos;
        while pos2 <= len and IsDigitChar( irratname[ pos2 ] ) do
          pos2:= pos2 + 1;
        od;
        gal:= Int( irratname{ [ pos .. pos2-1 ] } );
        irrat:= GaloisCyc( irrat, gal );
        pos:= pos2;
      fi;
    fi;

    while pos <= len do

      # Get ampersand summands.
      if irratname[ pos ] = '&' then
        pos2:= pos + 1;
        while pos2 <= len and IsDigitChar( irratname[ pos2 ] ) do
          pos2:= pos2 + 1;
        od;
        gal:= Int( irratname{ [ pos+1 .. pos2-1 ] } );
        irrat:= irrat + coeff * GaloisCyc( qN, gal );
        pos:= pos2;
      elif irratname[ pos ] in "+-" then
        if irratname[ pos ] = '+' then
          sign:= 1;
          oldpos:= pos+1;
        else
          sign:= -1;
          oldpos:= pos;
        fi;
        pos2:= pos + 1;
        while pos2 <= len and IsDigitChar( irratname[ pos2 ] ) do
          pos2:= pos2 + 1;
        od;
        if pos2 = pos + 1 then
          coeff:= sign;
        else
          coeff:= sign * Int( irratname{ [ pos+1 .. pos2-1 ] } );
        fi;
        pos:= pos2;
        if pos <= len then
          if irratname[ pos ] = '&' then
            pos2:= pos + 1;
            while pos2 <= len and IsDigitChar( irratname[ pos2 ] ) do
              pos2:= pos2 + 1;
            od;
            gal:= Int( irratname{ [ pos+1 .. pos2-1 ] } );
            irrat:= irrat + coeff * GaloisCyc( qN, gal );
            pos:= pos2;
          else
            recurse:= AtlasIrrationality( irratname{ [ oldpos .. len ] } );
            if recurse = fail then
              return fail;
            fi;
            irrat:= irrat + recurse;
            pos:= len+1;
          fi;
        else
          irrat:= irrat + coeff;
        fi;
      else
        return fail;
      fi;

    od;

    # Return the result.
    return irrat;
end );


#############################################################################
##
#F  Quadratic( <cyc>[, <rat>] ) . information about quadratic irrationalities
##
InstallGlobalFunction( Quadratic, function( arg )
    local cyc,
          rat,
          denom,
          coeffs,     # Zumbroich basis coefficients of `cyc'
          facts,      # factors of conductor of `cyc'
          factsset,   # set of `facts'
          two_part,   # 2-part of the conductor of `cyc'
          root,       # `root' component of the result
          a,          # `a'    component of the result
          b,          # `b'    component of the result
          d,          # `d'    component of the result
          ATLAS,      # string that expresses `cyc' in {\sf ATLAS} format
          ATLAS2,     # another string, maybe shorter ...
          display;    # string that shows a way to input `cyc'

    cyc:= arg[1];
    rat:= Length( arg ) = 2 and arg[2] = true;

    if not IsCyc( cyc ) then
      return fail;
    fi;
    denom:= DenominatorCyc( cyc );
    if not ( rat or denom = 1 ) then
      return fail;
    fi;
    if IsRat( cyc ) then
      return rec(
                  a       := NumeratorRat( cyc ),
                  b       := 0,
                  root    := 1,
                  d       := denom,
                  ATLAS   := String( cyc ),
                  display := String( cyc )
                 );
    fi;
    if denom <> 1 then
      cyc:= cyc * denom;
    fi;

    coeffs:= ExtRepOfObj( cyc );
    facts:= Factors(Integers, Length( coeffs ) );
    factsset:= Set( facts );
    two_part:= Number( facts, x -> x = 2 );

    # Compute candidates for `a', `b', `root', `d'.
    if two_part = 0 and Length( facts ) = Length( factsset ) then

      root:= Length( coeffs );
      if root mod 4 = 3 then
        root:= -root;
      fi;
      a:= StarCyc( cyc );
      if a = fail then
        return fail;
      fi;

      # Set `a' to the trace of `cyc' over the rationals.
      a:= cyc + a;

      if Length( factsset ) mod 2 = 0 then
        b:= 2 * coeffs[2] - a;
      else
        b:= 2 * coeffs[2] + a;
      fi;
      if a mod 2 = 0 and b mod 2 = 0 then
        a:= a / 2;
        b:= b / 2;
        d:= 1;
      else
        d:= 2;
      fi;

    elif two_part = 2 and Length( facts ) = Length( factsset ) + 1 then

      root:= Length( coeffs ) / 4;
      if root = 1 then
        a:= coeffs[1];
        b:= - coeffs[2];
      else
        a:= coeffs[5];
        if Length( factsset ) mod 2 = 0 then a:= -a; fi;
        b:= - coeffs[ root + 5 ];
      fi;
      if root mod 4 = 1 then
        root:= -root;
        b:= -b;
      fi;
      d:= 1;

    elif two_part = 3 then

      root:= Length( coeffs ) / 4;
      if root = 2 then
        a:= coeffs[1];
        b:= coeffs[2];
        if b = coeffs[4] then
          root:= -2;
        fi;
      else
        a:= coeffs[9];
        if Length( factsset ) mod 2 = 0 then a:= -a; fi;
        b:= coeffs[ root / 2 + 9 ];
        if b <> - coeffs[ 3 * root / 2 - 7 ] then
          root:= -root;
        elif ( root / 2 ) mod 4 = 3 then
          b:= -b;
        fi;
      fi;
      d:= 1;

    else
      return fail;
    fi;

    # Check whether the candidates `a', `b', `d', `root' are correct.
    if d * cyc <> a + b * ER( root ) then
      return fail;
    fi;

    # Compute a string for the irrationality in {\ATLAS} format.
    if d = 2 then

      # Necessarily `root' is congruent 1 mod 4, only $b_{'root'}$ possible.
      # We have $'cyc' = ('a' + `b') / 2 + `b' b_{'root'}$.
      if a + b = 0 then
        if b = 1 then
          ATLAS:= "";
        elif b = -1 then
          ATLAS:= "-";
        else
          ATLAS:= ShallowCopy( String( b ) );
        fi;
      elif b = 1 then
        ATLAS:= Concatenation( String( ( a + b ) / 2 ), "+" );
      elif b = -1 then
        ATLAS:= Concatenation( String( ( a + b ) / 2 ), "-" );
      elif 0 < b then
        ATLAS:= Concatenation( String( ( a + b ) / 2 ), "+", String( b ) );
      else
        ATLAS:= Concatenation( String( ( a + b ) / 2 ), String( b ) );
      fi;

      Append( ATLAS, "b" );
      if 0 < root then
        Append( ATLAS, String( root ) );
      else
        Append( ATLAS, String( -root ) );
      fi;

    else

      # `d' = 1, so we may use $i_{'root'}$ and $r_{'root'}$.
      if a = 0 then
        ATLAS:= "";
      else
        ATLAS:= ShallowCopy( String( a ) );
      fi;
      if a <> 0 and b > 0 then Append( ATLAS, "+" ); fi;
      if b = -1 then
        Append( ATLAS, "-" );
      elif b <> 1 then
        Append( ATLAS, String( b ) );
      fi;
      if root > 0 then
        ATLAS:= Concatenation( ATLAS, "r", String( root ) );
      elif root = -1 then
        Append( ATLAS, "i" );
      else
        ATLAS:= Concatenation( ATLAS, "i", String( -root ) );
      fi;

      if ( root - 1 ) mod 4 = 0 then

        # In this case, also $b_{|'root'|}$ is possible.
        # Note that here the coefficients are never equal to $\pm 1$.
        if a = -b then
          ATLAS2:= String( 2 * b );
        else
          ATLAS2:= Concatenation( String( a+b ), "+", String( 2*b ) );
        fi;
        if root > 0 then
          ATLAS2:= Concatenation( ATLAS2, "b", String( root ) );
        else
          ATLAS2:= Concatenation( ATLAS2, "b", String( -root ) );
        fi;

        if Length( ATLAS2 ) < Length( ATLAS ) then
          ATLAS:= ATLAS2;
        fi;

      fi;

    fi;
    if denom <> 1 then
      ATLAS:= Concatenation( "(", ATLAS, ")/", String( denom ) );
    fi;
    ConvertToStringRep( ATLAS );

    # Compute a string used by the `Display' function for character tables.
    if a = 0 then
      if b = 1 then
        display:= "";
      elif b = -1 then
        display:= "-";
      else
        display:= Concatenation( String( b ), "*" );
      fi;
    elif b = 1 then
      display:= Concatenation( String( a ), "+" );
    elif b = -1 then
      display:= Concatenation( String( a ), "-" );
    elif 0 < b then
      display:= Concatenation( String( a ), "+", String( b ), "*" );
    else
      display:= Concatenation( String( a ), String( b ), "*" );
    fi;
    Append( display, Concatenation( "Sqrt(", String( root ), ")" ) );
    d:= d * denom;
    if d <> 1 then
      if a <> 0 then
        display:= Concatenation( "(", display, ")" );
      fi;
      display:= Concatenation( display, "/", String( d ) );
    fi;
    ConvertToStringRep( display );

    # Return the result.
    return rec(
                a       := a,
                b       := b,
                root    := root,
                d       := d,
                ATLAS   := ATLAS,
                display := display
               );
end );


#############################################################################
##
#M  GaloisMat( <mat> )  . . . . . . . . . . . . . for a matrix of cyclotomics
##
##  Note that we must not mix up plain lists and class functions, since
##  these objects lie in different families.
##
InstallMethod( GaloisMat,
    "for a matrix of cyclotomics",
    [ IsMatrix and IsCyclotomicCollColl ],
    function( mat )
    local warned,      # at most one warning will be printed if `mat' grows
          ncha,        # number of rows in `mat'
          nccl,        # number of columns in `mat'
          galoisfams,  # list with information about conjugate characters:
                       #       1 means rational character,
                       #      -1 means character with undefs,
                       #       0 means dependent irrational character,
                       #  [ .. ] means leading irrational character.
          n,           # conductor of irrationalities in `mat'
          genexp,      # generators of prime residues mod `n'
          generators,  # permutation of `mat' induced by elements in `genexp'
          X,           # one row of `mat'
          i, j,        # loop over rows of `mat'
          irrats,      # set of irrationalities in `X'
          fusion,      # positions of `irrats' in `X'
          k, l, m,     # loop variables
          generator,
          irratsimages,
          automs,
          family,
          orders,
          exp,
          image,
          oldorder,
          cosets,
          auto,
          conj,
          blocklength,
          innerlength;

    warned := false;
    ncha   := Length( mat );
    mat    := ShallowCopy( mat );

    # Step 1:
    # Find the minimal cyclotomic field $Q_n$ containing all irrational
    # entries of <mat>.

    galoisfams:= [];
    n:= 1;
    for i in [ 1 .. ncha ] do
      if ForAny( mat[i], IsUnknown ) then
        galoisfams[i]:= -1;
      elif ForAll( mat[i], IsRat ) then
        galoisfams[i]:= 1;
      else
        n:= LcmInt( n, Conductor( mat[i] ) );
      fi;
    od;

    # Step 2:
    # Compute generators for Aut( Q(n):Q ), that is,
    # compute generators for (Z/nZ)* and convert them to exponents.

    if 1 < n then

      # Each Galois automorphism induces a permutation of rows.
      # Compute the permutations for each generating automorphism.
      # (Initialize with the identity permutation.)
      genexp:= Flat( GeneratorsPrimeResidues( n ).generators );
      generators:= List( genexp, x -> [ 1 .. ncha ] );

    else

      # The matrix is rational.
      generators:= [];

    fi;

    # Step 3:
    # For each character X, find and complete the family of conjugates.

    if 0 < ncha then
      nccl:= Length( mat[1] );
    fi;

    for i in [ 1 .. ncha ] do
      if not IsBound( galoisfams[i] ) then

        # We have found an independent character that is not integral
        # and contains no unknowns.

        X:= mat[i];
        for j in [ i+1 .. ncha ] do
          if mat[j] = X then
            galoisfams[j]:= Unknown();
            Info( InfoWarning, 1,
                  "GaloisMat: row ", i, " is equal to row ", j );
          fi;
        od;

        # Initialize the list of distinct irrationalities of `X'
        # (not ordered).
        # Each Galois automorphism induces a permutation of that list
        # rather than of the entries of `X' themselves.
        irrats:= [];

        # Store how to distribute the entries of irrats to `X'.
        fusion:= [];

        for j in [ 1 .. nccl ] do
          if IsCyc( X[j] ) and not IsRat( X[j] ) then
            k:= 1;
            while k <= Length( irrats ) and X[j] <> irrats[k] do
              k:= k+1;
            od;
            if k > Length( irrats ) then
              # This is the first appearance of `X[j]' in `X'.
              irrats[k]:= X[j];
            fi;

            # Store the position in `irrats'.
            fusion[j]:= k;
          else
            fusion[j]:= 0;
          fi;
        od;

        irratsimages:= [ irrats ];
        automs:= [ 1 ];
        family:= [ i ]; # indices of family members (same ordering as automs)
        orders:= [];    # orders[k] will be the order of the k-th generator
        for j in [ 1 .. Length( genexp ) ] do
          exp:= genexp[j];
          image:= List( irrats, x -> GaloisCyc( x, exp ) );
          oldorder:= Length( automs );  # group order up to now
          cosets:= [];
          orders[j]:= 1;
          while not image in irratsimages do
            orders[j]:= orders[j] + 1;
            for k in [ 1 .. oldorder ] do
              auto:= ( automs[k] * exp ) mod n;
              image:= List( irrats, x -> GaloisCyc( x, auto ) );
              conj:= [];    # the conjugate character
              for l in [ 1 .. nccl ] do
                if fusion[l] = 0 then
                  conj[l]:= X[l];
                else
                  conj[l]:= image[ fusion[l] ];
                fi;
              od;
              l:= Position( mat, conj, i );
              if l <> fail then

                galoisfams[l]:= 0;
                Add( family, l );
                for m in [ l+1 .. ncha ] do
                  if mat[m] = conj then galoisfams[m]:= 0; fi;
                od;

              else

                if not warned and 500 < Length( mat ) then
                  Info( InfoWarning, 1,
                        "GaloisMat: completion of <mat> will have",
                        " more than 500 rows" );
                  warned:= true;
                fi;

                Add( mat, conj );
                galoisfams[ Length( mat ) ]:= 0;
                Add( family, Length( mat ) );

              fi;
              Add( automs, auto );
              Add( cosets, image );
            od;
            exp:= exp * genexp[j];
            image:= List( irrats, x -> GaloisCyc( x, exp ) );
          od;
          irratsimages:= Concatenation( irratsimages, cosets );
        od;

        # Store the conjugates and automorphisms.
        galoisfams[i]:= [ family, automs ];

        # Now the length of `family' is the size of the Galois family of the
        # row `X'.
        # Compute the permutation operation of the generators on the set of
        # rows in `family'.

        blocklength:= 1;
        for j in [ 1 .. Length( genexp ) ] do

          innerlength:= blocklength;
          blocklength:= blocklength * orders[j];
          generator:= [ 1 .. blocklength ];

          # `genexp[j]' maps the conjugates under the action of
          # $\langle `genexp[1]', \ldots, `genexp[j-1]' \rangle$
          # (a set of length `innerlength') as one block to their images,
          # preserving the order of succession.

          for l in [ 1 .. blocklength - innerlength ] do
            generator[l]:= l + innerlength;
          od;

          # Compute how a power of `genexp[j]' maps back to the block.

          exp:= ( genexp[j] ^ orders[j] ) mod n;
          for l in [ 1 .. innerlength ] do
            generator[ l + blocklength - innerlength ]:=
                 Position( irratsimages, List( irrats,
                             x -> GaloisCyc( x, exp*automs[l] ) ) );
          od;

          # Transfer this operation to the cosets under the operation of
          # $\langle `genexp[j+1]',\ldots,`genexp[Length(genexp)]' \rangle$,
          # and transfer this to <mat> via `family'.

          for k in [ 0 .. Length( family ) / blocklength - 1 ] do
            for l in [ 1 .. blocklength ] do
              generators[j][ family[ l + k*blocklength ] ]:=
                           family[ generator[ l ] + k*blocklength ];
            od;
          od;

        od;

      fi;
    od;

    # Convert the `generators' component to a set of generating permutations.
    generators:= Set( generators, PermList );
    RemoveSet( generators, () );  # `generators' arose from `PermList'
    if IsEmpty( generators ) then
      generators:= [ () ];  # `generators' arose from `PermList'
    fi;

    # Return the result.
    return rec(
                mat        := mat,
                galoisfams := galoisfams,
                generators := generators
               );
    end );


#############################################################################
##
#M  RationalizedMat( <mat> )  . . . . . .  list of rationalized rows of <mat>
##
InstallMethod( RationalizedMat,
    "for a matrix",
    [ IsMatrix ],
    function( mat )
    local i, rationalizedmat, rationalfams;

    rationalfams:= GaloisMat( mat );
    mat:= rationalfams.mat;
    rationalfams:= rationalfams.galoisfams;
    rationalizedmat:= [];
    for i in [ 1 .. Length( mat ) ] do
      if rationalfams[i] = 1 or rationalfams[i] = -1 then
        # The row is rational or contains unknowns.
        Add( rationalizedmat, mat[i] );
      elif IsList( rationalfams[i] ) then
        # The row is a leading character of a family.
        Add( rationalizedmat, Sum( mat{ rationalfams[i][1] } ) );
      fi;
    od;
    return rationalizedmat;
    end );


InstallOtherMethod( RationalizedMat,
    "for an empty list",
    [ IsList and IsEmpty ],
    empty -> [] );


#############################################################################
##
#M  IsGeneratorsOfMagmaWithInverses( <cycs> ) . .  for a coll. of cyclotomics
##
##  Disallow to create groups of cyclotomics because the `\^' operator has
##  a meaning for cyclotomics that makes it not compatible with that for
##  groups:
##  `E(4)^-1' is the *inverse* of `E(4)',
##  but in the group generated by `E(4)',
##  {\GAP} would use this term for the *conjugate* of `E(4)' by `-1'.
##
InstallMethod( IsGeneratorsOfMagmaWithInverses,
    "for a collection of cyclotomics (return false)",
    [ IsCyclotomicCollection ],
    SUM_FLAGS, # override everything else
    function( gens )
    Info( InfoWarning, 1,
          "no groups of cyclotomics allowed because of incompatible ^" );
    return false;
    end );


#############################################################################
##
##  Functions for factorizing polynomials over fields of cyclotomics
##  (For arbitrary number fields, the explicit embedding of the field of
##  rationals may cause that the functions below do not work.)
##


#############################################################################
##
#M  FactorsSquarefree( <R>, <cycpol>, <opt> )
##
##  The function uses Algorithm~3.6.4 in~\cite{Coh93}.
##  (The record <opt> is ignored.)
##
InstallMethod( FactorsSquarefree,
    "for a polynomial over a field of cyclotomics",
    IsCollsElmsX,
    [ IsAbelianNumberFieldPolynomialRing, IsUnivariatePolynomial, IsRecord ],
    function( R, U, opt )
    local coeffring, theta, xind, yind, x, y, T, coeffs, G, powers, pow, i,
          B, c, val, j, k, N, factors;

    if IsRationalsPolynomialRing( R ) then
      TryNextMethod();
    fi;

    # Let $K = \Q(\theta)$ be a number field,
    # $T \in \Q[X]$ the minimal monic polynomial of $\theta$.
    # Let $U(X) be a monic squarefree polynomial in $K[x]$.
    coeffring:= CoefficientsRing( R );
    theta:= PrimitiveElement( coeffring );

    xind:= IndeterminateNumberOfUnivariateRationalFunction( U );
    if xind = 1 then
      yind:= 2;
    else
      yind:= 1;
    fi;
    x:= Indeterminate( Rationals, xind );
    y:= Indeterminate( Rationals, yind );

    # Let $U(X) = \sum_{i=0}^m u_i X^i$ and write $u_i = g_i(\theta)$
    # for some polynomial $g_i \in \Q[X]$.
    # Set $G(X,Y) = \sum_{i=0}^m g_i(Y) X^i \in \Q[X,Y]$.
    coeffs:= CoefficientsOfUnivariatePolynomial( U );
    if ForAll( coeffs, IsRat ) then
      G:= U;
    else
      powers:= [ 1 ];
      pow:= 1;
      for i in [ 2 .. DegreeOverPrimeField( coeffring ) ] do
        pow:= pow * theta;
        powers[i]:= pow;
      od;
      B:= Basis( coeffring, powers );
      G:= Zero( U );
      for i in [ 1 .. Length( coeffs ) ] do
        if IsRat( coeffs[i] ) then
          G:= G + coeffs[i] * x^i;
        else
          c:= Coefficients( B, coeffs[i] );
          val:= c[1];
          for j in [ 2 .. Length( c ) ] do
            val:= val + c[j] * y^(j-1);
          od;
          G:= G + val * x^i;
        fi;
      od;
    fi;

    # Set $k = 0$.
    k:= 0;

    # Compute $N(X) = R_Y( T(Y), G(X - kY,Y) )$
    # where $R_Y$ denotes the resultant with respect to the variable $Y$.
    # If $N(X)$ is not squarefree, increase $k$.
    T:= MinimalPolynomial( Rationals, theta, yind );
    repeat
      k:= k+1;
      N:= Resultant( T, Value( G, [ x, y ], [ x-k*y, y ] ), y );
    until DegreeOfUnivariateLaurentPolynomial( Gcd( N, Derivative(N) ) ) = 0;

    # Let $N = \prod_{i=1}^g N_i$ be a factorization of $N$.
    # For $1 \leq i \leq g$, set $A_i(X) = \gcd( U(X), N_i(X + k \theta) )$.
    # The desired factorization of $U(X)$ is $\prod_{i=1}^g A_i$.
    factors:= List( Factors( PolynomialRing( Rationals, [ xind ] ), N ),
                    f -> Gcd( R, U, Value( f, x + k*theta ) ) );
    return Filtered( factors,
                     x -> DegreeOfUnivariateLaurentPolynomial( x ) <> 0 );
    end );


#############################################################################
##
#M  Factors( <R>, <cycpol> )  .  for a polynomial over a field of cyclotomics
##
InstallMethod( Factors,
    "for a polynomial over a field of cyclotomics",
    IsCollsElms,
    [ IsAbelianNumberFieldPolynomialRing, IsUnivariatePolynomial ],
    function( R, pol )
    local irrfacs, coeffring, i, factors, ind, coeffs, val,
          lc, der, g, factor, q;

    if IsRationalsPolynomialRing( R ) then
      TryNextMethod();
    fi;

    # Check whether the desired factorization is already stored.
    irrfacs:= IrrFacsPol( pol );
    coeffring:= CoefficientsRing( R );
    i:= PositionProperty( irrfacs, pair -> pair[1] = coeffring );
    if i <> fail then
      return ShallowCopy(irrfacs[i][2]);
    fi;

    # Handle (at most) linear polynomials.
    if DegreeOfLaurentPolynomial( pol ) < 2  then
      factors:= [ pol ];
      StoreFactorsPol( coeffring, pol, factors );
      return factors;
    fi;

    # Compute the valuation, split off the indeterminate as a zero.
    ind:= IndeterminateNumberOfLaurentPolynomial( pol );
    coeffs:= CoefficientsOfLaurentPolynomial( pol );
    val:= coeffs[2];
    coeffs:= coeffs[1];
    factors:= ListWithIdenticalEntries( val,
                  IndeterminateOfUnivariateRationalFunction( pol ) );

    if Length( coeffs ) = 1 then

      # The polynomial is a power of the indeterminate.
      factors[1]:= coeffs[1] * factors[1];
      StoreFactorsPol( coeffring, pol, factors );
      return factors;

    elif Length( coeffs ) = 2 then

      # The polynomial is a linear polynomial times a power of the indet.
      factors[1]:= coeffs[2] * factors[1];
      factors[ val+1 ]:= LaurentPolynomialByExtRepNC( FamilyObj( pol ),
                             [ coeffs[1] / coeffs[2], 1 ], 0, ind );
      StoreFactorsPol( coeffring, pol, factors );
      return factors;

    fi;

    # We really have to compute the factorization.
    # First split the polynomial into leading coefficient and monic part.
    lc:= Last(coeffs);
    if not IsOne( lc ) then
      coeffs:= coeffs / lc;
    fi;
    if val = 0 then
      pol:= pol / lc;
    else
      pol:= LaurentPolynomialByExtRepNC( FamilyObj( pol ), coeffs, 0, ind );
    fi;

    # Now compute the quotient of `pol' by the g.c.d. with its derivative,
    # and factorize the squarefree part.
    der:= Derivative( pol );
    g:= Gcd( R, pol, der );
    if DegreeOfLaurentPolynomial( g ) = 0 then
      Append( factors, FactorsSquarefree( R, pol, rec() ) );
    else
      for factor in FactorsSquarefree( R, Quotient( R, pol, g ), rec() ) do
        Add( factors, factor );
        q:= Quotient( R, g, factor );
        while q <> fail do
          Add( factors, factor );
          g:= q;
          q:= Quotient( R, g, factor );
        od;
      od;
    fi;

    # Sort the factorization.
    Sort( factors );

    # Adjust the first factor by the constant term.
    Assert( 2, DegreeOfLaurentPolynomial(g) = 0 );
    if not IsOne( g ) then
      lc:= g * lc;
    fi;
    if not IsOne( lc ) then
      factors[1]:= lc * factors[1];
    fi;

    # Store the factorization.
    Assert( 2, Product( factors ) = lc * pol * IndeterminateOfUnivariateRationalFunction( pol )^val );
    StoreFactorsPol( coeffring, pol, factors );

    # Return the factorization.
    return factors;
    end );


#############################################################################
##
#F  DenominatorCyc( <cyc> )
##
InstallGlobalFunction( DenominatorCyc, function( cyc )
    if IsRat( cyc ) then
      return DenominatorRat( cyc );
    else
      return Lcm( List( COEFFS_CYC( cyc ), DenominatorRat ) );
    fi;
    end );


#############################################################################


#
# The following code is meant to allow comparisons between some select
# cyclotomics domains. So you can do things like
#   Integers = GaussianRationals;
#   IsSubset(Rationals, PositiveIntegers);
# without GAP running into an error. However, this code currently only
# works for a small fixed set of domains. It will not, for example, work
# with cyclotomic field extensions, or manually defined rings over the
# integers such as ClosureRing(1, E(3)).
# It would be nice if we eventually, were able to compare and intersect
# such objects, too.
#

# The following are three lists of equal length. At position i of the first
# list is a certain known cyclotomic (semi)ring. At position i of the second
# list is the corresponding filter. At position i of the third list is a
# finite list which can act as a "proxy" for the corresponding semiring when
# it comes to comparing it for inclusion resp. computing intersections with
# any of the other semirings.
# To simplify the code using these lists, the final entry of each list is
# fail, resp. the trivial filter IsObject.
BindGlobal("CompareCyclotomicCollectionHelper_Semirings", MakeImmutable( [
        PositiveIntegers, NonnegativeIntegers,
        Integers, GaussianIntegers,
        Rationals, GaussianRationals,
        Cyclotomics, fail
] ) );

BindGlobal("CompareCyclotomicCollectionHelper_Filters", MakeImmutable( [
        IsPositiveIntegers, IsNonnegativeIntegers,
        IsIntegers, IsGaussianIntegers,
        IsRationals, IsGaussianRationals,
        HasIsWholeFamily and IsWholeFamily, IsObject
] ) );

BindGlobal("CompareCyclotomicCollectionHelper_Proxies", MakeImmutable( [
        [ 1 ], [ 0, 1 ],
        [ -1, 0, 1 ], [ -1, 0, 1, E(4) ],
        [ -1, 0, 1/2, 1 ], [ -1, 0, 1/2, 1, E(4), 1/2+E(4) ],
        [ -1, 0, 1/2, 1, E(4), 1/2+E(4), E(9) ], fail
] ) );

ForAll(CompareCyclotomicCollectionHelper_Proxies, IsSet);
if IsHPCGAP then
    MakeReadOnlySingleObj(CompareCyclotomicCollectionHelper_Semirings);
    MakeReadOnlyObj(CompareCyclotomicCollectionHelper_Filters);
fi;


BindGlobal("CompareCyclotomicCollectionHelper", function (A, B)
  local a, b;
  a := PositionProperty( CompareCyclotomicCollectionHelper_Filters, p -> p(A) );
  b := PositionProperty( CompareCyclotomicCollectionHelper_Filters, p -> p(B) );
  return CompareCyclotomicCollectionHelper_Proxies{[a,b]};
end );


InstallMethod( \=, "for certain cyclotomic semirings",
             [IsCyclotomicCollection and IsSemiringWithOne,
              IsCyclotomicCollection and IsSemiringWithOne],
function (A,B)
  local ab;
  ab := CompareCyclotomicCollectionHelper(A, B);
  # It suffices if we "recognize" at least on of A and B; but if we
  # recognize neither, we give up.
  if ab = [fail,fail] then TryNextMethod(); fi;
  return ab[1] = ab[2];
end );


InstallMethod( IsSubset, "for certain cyclotomic semirings",
             [IsCyclotomicCollection and IsSemiringWithOne,
              IsCyclotomicCollection and IsSemiringWithOne],
function (A,B)
  local ab;
  ab := CompareCyclotomicCollectionHelper(A, B);
  # Verify that we recognized both A and B, otherwise give up.
  if fail in ab then TryNextMethod(); fi;
  return IsSubset(ab[1], ab[2]);
end );


InstallMethod( Intersection2, "for certain cyclotomic semirings",
             [IsCyclotomicCollection and IsSemiringWithOne,
              IsCyclotomicCollection and IsSemiringWithOne],
function (A,B)
  local ab, i;
  ab := CompareCyclotomicCollectionHelper(A, B);
  # Verify that we recognized both A and B, otherwise give up.
  if fail in ab then TryNextMethod(); fi;
  i := Position( CompareCyclotomicCollectionHelper_Proxies, Intersection2( ab[1], ab[2] ) );
  return CompareCyclotomicCollectionHelper_Semirings[i];
end );


InstallMethod( Union2, "for certain cyclotomic semirings",
             [IsCyclotomicCollection and IsSemiringWithOne,
              IsCyclotomicCollection and IsSemiringWithOne],
function (A,B)
  local ab, i;
  ab := CompareCyclotomicCollectionHelper(A, B);
  # Verify that we recognized both A and B, otherwise give up.
  if fail in ab then TryNextMethod(); fi;
  i := Position( CompareCyclotomicCollectionHelper_Proxies, Union2( ab[1], ab[2] ) );
  if i = fail then TryNextMethod(); fi;
  return CompareCyclotomicCollectionHelper_Semirings[i];
end );