home *** CD-ROM | disk | FTP | other *** search
/ Math Solutions 1995 October / Math_Solutions_CD-ROM_Walnut_Creek_October_1995.iso / pc / mac / discrete / lib / numfield.g < prev    next >
Encoding:
Text File  |  1993-05-05  |  52.5 KB  |  1,674 lines
  1. #############################################################################
  2. ##
  3. #A  numfield.g                  GAP library                     Thomas Breuer
  4. ##
  5. #A  @(#)$Id: numfield.g,v 3.9 1993/02/09 14:25:55 martin Rel $
  6. ##
  7. #Y  Copyright 1990-1992,  Lehrstuhl D fuer Mathematik,  RWTH Aachen,  Germany
  8. ##
  9. ##  This file contains those functions   that deal with  number fields,  that
  10. ##  is   subfields  of cyclotomic  fields,   and  their  elements  (which are
  11. ##  cyclotomics, see 'cyclotom.g')
  12. ##
  13. #H  $Log: numfield.g,v $
  14. #H  Revision 3.9  1993/02/09  14:25:55  martin
  15. #H  made undefined globals local
  16. #H
  17. #H  Revision 3.8  1993/02/04  11:10:55  sam
  18. #H  corrections mainly concerning CF(4) and GaussianRationals
  19. #H
  20. #H  Revision 3.7  1992/12/16  19:47:27  martin
  21. #H  replaced quoted record names with escaped ones
  22. #H
  23. #H  Revision 3.6  1992/11/02  12:58:14  sam
  24. #H  fixed bug in 'OrderCyc'
  25. #H
  26. #H  Revision 3.5  1992/08/11  17:21:45  sam
  27. #H  changed 'Print' functions of number fields, fixed minor bug
  28. #H
  29. #H  Revision 3.4  1992/03/27  11:14:51  martin
  30. #H  changed mapping to general mapping and function to mapping
  31. #H
  32. #H  Revision 3.3  1992/02/13  14:22:16  sam
  33. #H  added 'GaloisGroup'
  34. #H
  35. #H  Revision 3.2  1991/12/04  11:14:09  sam
  36. #H  corrected errors in 'NumberRingOps'
  37. #H
  38. #H  Revision 3.1  1991/12/04  10:22:20  sam
  39. #H  changed succession of functions (necess. to overlay correctly)
  40. #H
  41. #H  Revision 3.0  1991/12/04  09:14:35  martin
  42. #H  initial revision under RCS
  43. #H
  44. ##
  45.  
  46.  
  47. #############################################################################
  48. ##
  49. #F  IsNumberField( <obj> )  . . .  test if <obj> is a cyclotomic field record
  50. ##
  51. IsNumberField := function ( obj )
  52.     return IsRec( obj )
  53.            and IsBound( obj.isField ) and obj.isField
  54.            and obj.char = 0
  55.            and ForAll( obj.generators, IsCyc );
  56.     end;
  57.  
  58.  
  59. #############################################################################
  60. ##
  61. #F  IsCyclotomicField( <obj> )  .  test if <obj> is a cyclotomic field record
  62. ##
  63. IsCyclotomicField := function ( obj )
  64.     return IsNumberField( obj )
  65.            and IsBound( obj.isCyclotomicField ) and obj.isCyclotomicField;
  66.     end;
  67.  
  68.  
  69. #############################################################################
  70. ##
  71. #V  NumberFieldOps  . . . . . . . . . . . operations record for number fields
  72. ##
  73. NumberFieldOps := Copy( FieldOps );
  74.  
  75.  
  76. #############################################################################
  77. ##
  78. #F  NumberFieldOps.\in( <z>, <F> ) . . . . . . . .   test if <z> lies in <F>
  79. ##
  80. NumberFieldOps.\in := function ( z, F )
  81.     return IsCyc( z )
  82.            and NofCyc( F.generators ) mod NofCyc( z ) = 0  # same envelopping
  83.                                                            # cyclotomic field
  84.            and ForAll( F.stabilizer, x -> GaloisCyc( z, x ) = z );
  85.     end;
  86.  
  87.  
  88. #############################################################################
  89. ##
  90. #F  NumberFieldOps.Intersection( <F>, <G> ) . . intersection of number fields
  91. ##
  92. NumberFieldOps.Intersection := function ( F, G )
  93.     local i, j, N, stab, stabF, stabG;
  94.     
  95.     if IsNumberField( F ) and IsNumberField( G ) then
  96.     
  97.       if IsCyclotomicField( F ) then
  98.     
  99.         if IsCyclotomicField( G ) then
  100.     
  101.           return CF( GcdInt( NofCyc(F.generators), NofCyc(G.generators) ) );
  102.     
  103.         else
  104.     
  105.           # intersection of cyclotomic field 'F = FC(N)' with number field
  106.           # 'G'; simply take the elements of 'G.stabilizer' modulo 'N'.
  107.           # (If a reduction is necessary, 'NF' will do.)
  108.     
  109.           N:= NofCyc( F.generators );
  110.           return NF( N, Set( List( G.stabilizer, x -> x mod N ) ) );
  111.     
  112.         fi;
  113.     
  114.       elif IsCyclotomicField( G ) then
  115.     
  116.         # intersection of cyclotomic field 'G = FC(N)' with number field 'F';
  117.         # simply take the elements of 'F.stabilizer' modulo 'N'.
  118.         # (If a reduction is necessary, 'NF' will do.)
  119.     
  120.         N:= NofCyc( G.generators );
  121.         return NF( N, Set( List( F.stabilizer, x -> x mod N ) ) );
  122.     
  123.       else
  124.     
  125.         # first compute 'N' where 'CF(N)' contains the intersection;
  126.         # reduce the elements of the stabilizers modulo 'N', i.e. intersect
  127.         # 'F' and 'G' with 'CF(N)';
  128.         # then compute the corresponding stabilizer, i.e. the product of
  129.         # stabilizers.
  130.         N:= GcdInt( NofCyc( F.generators), NofCyc( G.generators ) );
  131.         stabF:= Set( List( F.stabilizer, x -> x mod N ) );
  132.         stabG:= Set( List( G.stabilizer, x -> x mod N ) );
  133.         stab:= [];
  134.         for i in stabF do
  135.           for j in stabG do
  136.             AddSet( stab, ( i * j ) mod N );
  137.           od;
  138.         od;
  139.     
  140.         # (If a reduction is necessary, 'NF' will do.)
  141.     
  142.         return NF( N, stab );
  143.     
  144.       fi;
  145.     else
  146.       return FieldOps.Intersection( F, G );
  147.     fi;
  148.     end;
  149.  
  150.  
  151. #############################################################################
  152. ##
  153. #F  NumberFieldOps.Random( <F> )  . . . . .  random element of a number field
  154. ##
  155. NumberFieldOps.Random := function ( F )
  156.  
  157.     # get random elements of the subfield and multiply with the base
  158.     if IsBound( F.base ) then
  159.       if F.field = 0 then
  160.         return List( F.base, x->Random(Rationals) ) * F.base;
  161.       else
  162.         return List( F.base, x->Random( F.field ) ) * F.base;
  163.       fi;
  164.     else
  165.       return FieldOps.Random( F );
  166.     fi;
  167.     end;
  168.  
  169.  
  170. #############################################################################
  171. ##
  172. #F  NumberFieldOps.Print( <F> ) . . . . . . . . . . . .  print a number field
  173. ##
  174. NumberFieldOps.Print := function ( F )
  175.     if IsBound( F.name ) then
  176.       Print( F.name );
  177.     elif F.field = 0 then
  178.       if IsBound( F.basechangemat ) then
  179.         Print( "NF( 0,", F.base, ")" );
  180.       else
  181.         Print( "NF(", NofCyc( F.generators ), ",", F.stabilizer, ")" );
  182.       fi;
  183.     else
  184.       if IsBound( F.basechangemat ) then
  185.         Print( "NF(", F.field, ",", F.base, ")" );
  186.       else
  187.         Print( "NF(", NofCyc( F.generators ), ",", F.stabilizer, ")/",
  188.                F.field );
  189.       fi;
  190.     fi;
  191.     end;
  192.  
  193.  
  194. #############################################################################
  195. ##
  196. #F  NumberFieldOps.\/( <F>, <G> )  . . . . . . . . quotient of number fields
  197. ##
  198. NumberFieldOps.\/ := function( F, G )
  199.  
  200.     if not ForAll( G.generators, x -> x in F ) then
  201.       Error( "<G> is not a subfield of <F>" );
  202.     fi;
  203.  
  204.     if IsCyclotomicField( F ) then
  205.       return CF( G, NofCyc( F.generators ) );
  206.     else
  207.       if F.field = 0 or F.field = Rationals then
  208.         return NF( G, F.base );
  209.       else
  210.         return NF( G, Concatenation( List( F.base, x -> x*F.field.base ) ) );
  211.       fi;
  212.     fi;
  213.     end;
  214.  
  215.  
  216. #############################################################################
  217. ##
  218. #F  NumberFieldOps.GaloisGroup(<F>) . . . . .  Galois group of a number field
  219. ##
  220. NumberFieldOps.GaloisGroup := function ( F )
  221.     local group;
  222.     group:= Group( List( Flat(
  223.                 GeneratorsPrimeResidues( NofCyc(F.generators) ).generators ),
  224.                          x -> NFAutomorphism( F, x ) ),
  225.                    IdentityMapping( F ) );
  226.  
  227.     if IsInt( F.field )  then
  228.  
  229.         # Galois group of a number field over the rationals,
  230.         # simply a factor group of the Galois group of the conductor
  231.         return group;
  232.     else
  233.  
  234.         # take the subgroup of the prime residue group that fixes 'F.field'
  235.         # pointwise; the required group is a factor group of this group
  236.  
  237.         return Group( Stabilizer( group, F.field.generators, OnTuples ) );
  238.     fi;
  239. end;
  240.  
  241.  
  242. #############################################################################
  243. ##
  244. #F  IsNFAutomorphism(<obj>)  . . . . test if an object is a number field aut.
  245. ##
  246. IsNFAutomorphism := function ( obj )
  247.     return IsRec( obj )
  248.        and IsBound( obj.isNFAutomorphism )
  249.        and obj.isNFAutomorphism;
  250. end;
  251.  
  252.  
  253. #############################################################################
  254. ##
  255. #F  NFAutomorphism( <F>, <k> ) . . . . . . . . automorphism of a number field
  256. #V  NFAutomorphismOps  . . . . . . . operations record for number field auts.
  257. ##
  258. NFAutomorphism := function ( F, k )
  259.  
  260.     # check the arguments
  261.     if not IsNumberField( F ) then Error("<F> must be a number field"); fi;
  262.  
  263.     # Let $F / K$ be a field extension where $Q_n$ is the conductor of $F$;
  264.     # let $S(F)$ be the group of those prime residues mod $n$ that fix $F$
  265.     # pointwise.  The Galois group of $F / K$ is in natural correspondence
  266.     # to $S(K) / S(F)$.  Thus each automorphism of $F / K$ corresponds to
  267.     # a coset $c$ of $S(K)$, and it acts on $F$ like each element of $c$.
  268.     # The automorphism 'NFAutomorphism( F/K, k )' maps $x\in F / K$ to
  269.     # 'GaloisCyc( <x>, k )'.
  270.  
  271.  
  272.     # make the mapping record
  273.     return rec(
  274.         isGeneralMapping        := true,
  275.         domain                  := Mappings,
  276.  
  277.         isNFAutomorphism        := true,
  278.         galois                  := Set(List(F.stabilizer,
  279.                                        x->x*k mod NofCyc(F.generators)))[1],
  280.         source                  := F,
  281.         range                   := F,
  282.  
  283.         isMapping               := true,
  284.         isHomomorphism          := true,
  285.         isSurjective            := true,
  286.         isInjective             := true,
  287.         isBijection             := true,
  288.  
  289.         operations              := NFAutomorphismOps );
  290.  
  291. end;
  292.  
  293. NFAutomorphismOps := Copy( FieldHomomorphismOps );
  294.  
  295. NFAutomorphismOps.\= := function ( aut1, aut2 )
  296.     if IsNFAutomorphism( aut1 )  then
  297.         if IsNFAutomorphism( aut2 )  then
  298.             return aut1.source  = aut2.source
  299.                and aut1.galois  = aut2.galois;
  300.         else
  301.             return false;
  302.         fi;
  303.     else
  304.         if IsNFAutomorphism( aut2 )  then
  305.             return false;
  306.         else
  307.             Error("panic: neither <aut1> nor <aut2> is a NF aut.");
  308.         fi;
  309.     fi;
  310. end;
  311.  
  312. NFAutomorphismOps.ImageElm := function ( aut, elm )
  313.     return GaloisCyc( elm, aut.galois );
  314. end;
  315.  
  316. NFAutomorphismOps.ImagesElm := function ( aut, elm )
  317.     return [ GaloisCyc( elm, aut.galois ) ];
  318. end;
  319.  
  320. NFAutomorphismOps.ImagesSet := function ( aut, elms )
  321.     if IsField( elms )  and  IsSubset( aut.source, elms )  then
  322.       if IsInt( elms.field ) then
  323.         return NF( List( elms.generators,
  324.                          x->GaloisCyc(x,aut.galois) ) );
  325.       else
  326.         return NF( List( elms.generators,
  327.                          x->GaloisCyc(x,aut.galois) ) ) / 
  328.                NF( List( elms.field.generators,
  329.                          x->GaloisCyc(x,aut.galois) ) );
  330.       fi;
  331.     else
  332.         return FieldHomomorphismOps.ImagesSet( aut, elms );
  333.     fi;
  334. end;
  335.  
  336. NFAutomorphismOps.ImagesRepresentative := function ( aut, elm )
  337.     return GaloisCyc( elm, aut.galois );
  338. end;
  339.  
  340. NFAutomorphismOps.CompositionMapping := function ( aut1, aut2 )
  341.     if IsNFAutomorphism( aut1 )  then
  342.         if IsNFAutomorphism( aut2 ) and aut1.source = aut2.source then
  343.             return NFAutomorphism( aut1.source,
  344.                                            aut1.galois * aut2.galois );
  345.         else
  346.             return FieldHomomorphismOps.CompositionMapping( aut1, aut2 );
  347.         fi;
  348.     else
  349.         return FieldHomomorphismOps.CompositionMapping( aut1, aut2 );
  350.     fi;
  351. end;
  352.  
  353. NFAutomorphismOps.InverseMapping := function ( aut )
  354.     return NFAutomorphism( aut.source, 1 / aut.galois );
  355. end;
  356.  
  357. NFAutomorphismOps.PowerMapping := function ( aut, i )
  358.     return NFAutomorphism( aut.source, aut.galois^i );
  359. end;
  360.  
  361. NFAutomorphismOps.\< := function ( aut1, aut2 )
  362.     if IsNFAutomorphism( aut1 )  then
  363.         if IsNFAutomorphism( aut2 )  then
  364.             if aut1.source <> aut2.source  then
  365.                 return aut1.source < aut2.source;
  366.             else
  367.                 return aut1.galois < aut2.galois;
  368.             fi;
  369.         else
  370.             return FieldHomomorphismOps.\<( aut1, aut2 );
  371.         fi;
  372.     else
  373.         return FieldHomomorphismOps.\<( aut1, aut2 );
  374.     fi;
  375. end;
  376.  
  377. NFAutomorphismOps.Print := function ( aut )
  378.     if aut.galois = 1  then
  379.         Print( "IdentityMapping( ", aut.source, " )" );
  380.     else
  381.         Print( "NFAutomorphism( ", aut.source, " , ", aut.galois, " )" );
  382.     fi;
  383. end;
  384.  
  385. NumberFieldOps.IdentityMapping := F -> NFAutomorphism( F, 1 );
  386.  
  387.  
  388. #############################################################################
  389. ##
  390. #F  NumberFieldOps.Conjugates( <F>, <z> ) .  all conjugates of an alg. number
  391. ##
  392. NumberFieldOps.Conjugates := function ( F, z )
  393.     local N, gal, result, pnt;
  394.     
  395.     N:= NofCyc( F.generators );
  396.     
  397.     # automorphisms of the envelopping cyclotomic field
  398.     gal:= PrimeResidues( N );
  399.     
  400.     if IsNumberField( F.field ) and F.field <> Rationals then
  401.     
  402.       # take only the subgroup of 'gal' that fixes 'F.field'
  403.       gal:= Filtered( gal,
  404.                       x -> ForAll(F.field.generators,y->GaloisCyc(y,x)=y) );
  405.     fi;
  406.     
  407.     # get representatives of cosets of 'F.stabilizer'
  408.     result:= [];
  409.     while gal <> [] do
  410.       pnt:= gal[1];
  411.       Add( result, GaloisCyc( z, pnt ) );
  412.       SubtractSet( gal, List( F.stabilizer, x -> ( x * pnt ) mod N ) );
  413.     od;
  414.     
  415.     return result;
  416.     end;
  417.  
  418.  
  419. #############################################################################
  420. ##
  421. #F  NumberFieldOps.Norm( <F>, <z> ) . . . . . . . . .  norm of an alg. number
  422. ##
  423. NumberFieldOps.Norm := function ( F, z )
  424.     local N, gal, result, pnt;
  425.     
  426.     N:= NofCyc( F.generators );
  427.     
  428.     # automorphisms of the envelopping cyclotomic field
  429.     gal:= PrimeResidues( N );
  430.     
  431.     if IsNumberField( F.field ) and F.field <> Rationals then
  432.     
  433.       # take only the subgroup of 'gal' that fixes 'F.field'
  434.       gal:= Filtered( gal,
  435.                       x -> ForAll(F.field.generators,y->GaloisCyc(y,x)=y) );
  436.     fi;
  437.     
  438.     # get representatives of cosets of 'F.stabilizer'
  439.     result:= 1;
  440.     while gal <> [] do
  441.       pnt:= gal[1];
  442.       result:= result * GaloisCyc( z, pnt );
  443.       SubtractSet( gal, List( F.stabilizer, x -> ( x * pnt ) mod N ) );
  444.     od;
  445.     
  446.     return result;
  447.     end;
  448.  
  449.  
  450. #############################################################################
  451. ##
  452. #F  NumberFieldOps.Trace( <F>, <z> )  . . . . . . . . trace of an alg. number
  453. ##
  454. NumberFieldOps.Trace := function ( F, z )
  455.     local N, gal, result, pnt;
  456.     
  457.     N:= NofCyc( F.generators );
  458.     
  459.     # automorphisms of the envelopping cyclotomic field
  460.     gal:= PrimeResidues( N );
  461.     
  462.     if IsNumberField( F.field ) and F.field <> Rationals then
  463.     
  464.       # take only the subgroup of 'gal' that fixes 'F.field'
  465.       gal:= Filtered( gal,
  466.                       x -> ForAll(F.field.generators,y->GaloisCyc(y,x)=y) );
  467.     fi;
  468.     
  469.     # get representatives of cosets of 'F.stabilizer'
  470.     result:= 0;
  471.     while gal <> [] do
  472.       pnt:= gal[1];
  473.       result:= result + GaloisCyc( z, pnt );
  474.       SubtractSet( gal, List( F.stabilizer, x -> ( x * pnt ) mod N ) );
  475.     od;
  476.     
  477.     return result;
  478.     end;
  479.  
  480.  
  481. #############################################################################
  482. ##
  483. #F  NumberFieldOps.Order( <F>, <z> )  . . . . . . . . order of an alg. number
  484. ##
  485. OrderCyc := function ( cyc )
  486.     local ord, val;
  487.     if not IsCyc( cyc ) or cyc = 0 then
  488.       Error( "argument must be a nonzero cyclotomic" );
  489.     elif cyc ^ ( 2 * NofCyc( cyc ) ) <> 1 then   # not a root of unity
  490.       return "infinity";
  491.     else
  492.       ord:= 1;
  493.       val:= cyc;
  494.       while val <> 1 do
  495.         val:= val * cyc;
  496.         ord:= ord + 1;
  497.       od;
  498.       return ord;
  499.     fi;
  500.     end;
  501.  
  502. NumberFieldOps.Order := function ( F, z ) return OrderCyc( z ); end;
  503.  
  504.  
  505. #############################################################################
  506. ##
  507. #F  NormalBaseNumberField(<F>[,<x>])  . . . . . normal base of a number field
  508. ##
  509. ##  returns a normal base of the number field 'F' with respect to 'F.subfield'
  510. ##  (using the algorithm given in E. Artin, Galoissche Theorie, p. 65 f.).
  511. ##
  512. ##  The second form tries to evaluate the polynomial at <x>.
  513. ##
  514. NormalBaseNumberField := function( arg )
  515.     local F, alpha, poly, En, normal, i, val, stabilizer;
  516.     
  517.     if Length(arg) < 1 or Length(arg) > 2 or not IsNumberField(arg[1]) then
  518.       Error( "usage: NormalBaseNumberField(<F>,[<x>]) for number field <F>" );
  519.     fi;
  520.     F:= arg[1];
  521.     
  522.     if IsBound( F.stabilizer ) then
  523.       stabilizer:= F.stabilizer;
  524.     else
  525.       stabilizer:= Filtered( PrimeResidues( NofCyc( F.generators ) ),
  526.                              x -> ForAll(F.generators,y->GaloisCyc(y,x)=y) );
  527.     fi;
  528.     
  529.     # get a primitive element 'alpha'
  530.     if Length( F.generators ) = 1 then
  531.       alpha:= F.generators[1];
  532.     else
  533.       En:= E( NofCyc( F.generators ) );
  534.       alpha:= Sum( stabilizer, x -> GaloisCyc( En, x ) );
  535.     fi;
  536.     
  537.     # the polynomial
  538.     # $\prod_{\sigma\in 'Gal( alpha )'\setminus \{1\} } (x-\sigma('alpha') )
  539.     # for the primitive element 'alpha'
  540.     poly:= [ 1 ];
  541.     for i in Difference( F.operations.Conjugates( F, alpha ), [ alpha ] ) do
  542.       poly:= ProductPol( poly, [ -i, 1 ] );
  543.     od;
  544.     
  545.     # the denominator\: eval 'poly' at 'a'
  546.     val:= 1 / ValuePol( poly, alpha );
  547.     
  548.     # there are only finitely many values 'x' in 'F.field' for which
  549.     # 'poly(x) \* val' is not an element of a normal base.
  550.     if Length( arg ) = 2 then
  551.       i:= arg[2];
  552.       if not i in F.field then
  553.         Error( "<x> must be an element of <F>.field" );
  554.       fi;
  555.     else
  556.       i:= 1;
  557.     fi;
  558.  
  559.     repeat
  560.       normal:= F.operations.Conjugates( F, ValuePol( poly, i ) * val );
  561.       i:= i + 1;
  562.     until RankMat( List( normal, COEFFSCYC ) ) = F.dimension;
  563.     
  564.     return normal;
  565.     end;
  566.  
  567.  
  568. #############################################################################
  569. ##
  570. #F  ZumbroichBase( <n>, <m> )
  571. ##
  572. ##  returns the set of exponents 'e' for which 'E(n)^e' belongs to the
  573. ##  (generalized) Zumbroich base of the cyclotomic field $Q_n$,
  574. ##  viewed as vector space over $Q_m$.
  575. ##
  576. ##  *Note* that for $n \equiv 2 \bmod 4$ we have
  577. ##  'ZumbroichBase( <n>, 1 ) = 2 * ZumbroichBase( <n>/2, 1 )' but
  578. ##  'List( ZumbroichBase(  <n>, 1  ), x -> E(  <n>  )^x ) =
  579. ##   List( ZumbroichBase( <n>/2, 1 ), x -> E( <n>/2 )^x )'.
  580. ##
  581. ZumbroichBase := function( n, m )
  582.     local nn, base, basefactor, factsn, exponsn, factsm, exponsm, primes,
  583.           p, pos, i, k;
  584.     
  585.     if not n mod m = 0 then
  586.       Error( "<m> must be a divisor of <n>" );
  587.     fi;
  588.     
  589.     factsn:= FactorsInt( n );
  590.     primes:= Set( factsn );
  591.     exponsn:= List( primes, x -> 0 );   # Product(List( [1..Length(primes)],
  592.                                         #         x->primes[i]^exponsn[i]))=n
  593.     p:= factsn[1];
  594.     pos:= 1;
  595.     for i in factsn do
  596.       if i > p then
  597.         p:= i;
  598.         pos:= pos + 1;
  599.       fi;
  600.       exponsn[ pos ]:= exponsn[ pos ] + 1;
  601.     od;
  602.     
  603.     factsm:= FactorsInt( m );
  604.     exponsm:= List( primes, x -> 0 );    # Product(List( [1..Length(primes)],
  605.                                          #         x->primes[i]^exponsm[i]))=m
  606.     if m <> 1 then
  607.       p:= factsm[1];
  608.       pos:= Position( primes, p );
  609.       for i in factsm do
  610.         if i > p then
  611.           p:= i;
  612.           pos:= Position( primes, p );
  613.         fi;
  614.         exponsm[ pos ]:= exponsm[ pos ] + 1;
  615.       od; 
  616.     fi;
  617.     
  618.     base:= [ 0 ];
  619.     if n = 1 then return base; fi;
  620.       
  621.     if primes[1] = 2 then             # special case: $J_{k,2} = \{ 0, 1 \}$
  622.     
  623.       if exponsm[1] = 0 then exponsm[1]:= 1; fi;    # $J_{0,2} = \{ 0 \}$
  624.     
  625.       nn:= n / 2^( exponsm[1] + 1 );
  626.     
  627.       for k in [ exponsm[1] .. exponsn[1] - 1 ] do
  628.         Append( base, base + nn );
  629.         nn:= nn / 2;
  630.       od;
  631.       pos:= 2;
  632.     else
  633.       pos:= 1;
  634.     fi;
  635.     
  636.     for i in [ pos .. Length( primes ) ] do
  637.     
  638.       if m mod primes[i] <> 0 then
  639.         basefactor:= [ 1 .. primes[i] - 1 ] * ( n / primes[i] );
  640.         base:= Concatenation( List( base, x -> x + basefactor ) );
  641.         exponsm[i]:= 1;
  642.       fi;
  643.     
  644.       basefactor:= [ - ( primes[i] - 1 ) / 2 .. ( primes[i] - 1 ) / 2 ]
  645.                      * n / primes[i]^exponsm[i];
  646.     
  647.       for k in [ exponsm[i] .. exponsn[i] - 1 ] do
  648.         basefactor:= basefactor / primes[i];
  649.         base:= Concatenation( List( base, x -> x + basefactor ) );
  650.       od;
  651.     od;
  652.     return Set( List( base, x -> x mod n ) );
  653.     end;
  654.  
  655.  
  656. #############################################################################
  657. ##
  658. #F  LenstraBase(<N>,<stabilizer>,<super>) .   integral base of a number field
  659. ##
  660. ##  returns a list of lists of integers; each list indexing the exponents of
  661. ##  an orbit of a subgroup of <stabilizer> on <N>-th roots of unity.
  662. ##
  663. ##  *Note* that the elements are in general not sets, since the first element
  664. ##  is always an element of 'ZumbroichBase(<N>,1)'; this is used by 'NF' and
  665. ##  'Coefficients'.
  666. ##
  667. ##  <super> is a list representing a supergroup of <stabilizer> which
  668. ##  shall act consistently with the action of <stabilizer>, i.e. each orbit
  669. ##  of <supergroup> is a union of orbits of <stabilizer>.
  670. ##
  671. ##  ( Shall there be a test if this is possible ? )
  672. ##
  673. ##  *Note* that <stabilizer> must not contain the stabilizer of a proper
  674. ##  cyclotomic subfield of the <N>-th cyclotomic field.
  675. ##
  676. LenstraBase := function( N, stabilizer, supergroup )
  677.     local i, k, primes, NN, zumb, orbits, pnt, orb, d, N2, No, ppnt, ord, a,
  678.           neworbits, rep, super, H1;
  679.     
  680.     primes:= Set( FactorsInt( N ) );
  681.     NN:= Product( primes );                # squarefree part of 'N'
  682.     zumb:= ZumbroichBase( N, 1 );          # exps of roots in base of 'CF(N)'
  683.     stabilizer:= Set( stabilizer );
  684.     
  685.     if N = NN then
  686.     
  687.       # 'N' is squarefree, we have the normal base, 'stabilizer' acts on
  688.       # 'zumb'; do not consider equivalence classes since they are all
  689.       # trivial.  'supergroup' is obsolete since 'zumb' is a normal base.
  690.     
  691.       # *Note* that for even 'N' 'zumb' does not consist of at least 'NN'-th
  692.       # roots!
  693.     
  694.       orbits:= [];
  695.       while zumb <> [] do
  696.         pnt:= zumb[1];
  697.         orb:= List( stabilizer, x -> pnt * x mod N );
  698.         Add( orbits, orb );
  699.         SubtractSet( zumb, orb );
  700.       od;
  701.     
  702.     else
  703.     
  704.       # Let $d(i)$ be the largest squarefree number whose square divides the
  705.       # order of $e_N^i$, that is 'N / gcd(N,i)'.
  706.       # Define an equivalence relation on the set $S$ of at least 'NN'-th
  707.       # roots of unity\:
  708.       # $i$ and $j$ are equivalent iff 'N' divides $( i - j ) d(i)$.  The
  709.       # equivalence class $(i)$ of $i$ is $\{ i+kN/d(i) ; 0\leq k\< d(i)\}$.
  710.     
  711.       # For the case that 'NN' is even, replace those roots in $S$ with order
  712.       # not divisible by 4 by their negatives. (Equivalently\: Replace *all*
  713.       # elements in $S$ by their negatives.)
  714.     
  715.       # If 8 does not divide 'N' and $'N' \not= 4$, 'zumb' is a subset of $S$,
  716.       # the intersection of $(i)$ with 'zumb' is of order $\varphi( d(i) )$,
  717.       # it is a basis for the $Z$--submodule spanned by $(i)$.
  718.       # Furthermore, the minimality of 'N' yields that 'stabilizer' acts fixed
  719.       # point freely on the set of equivalence classes.
  720.     
  721.       # More exactly, fixed points occur exactly if there is an element 's' in
  722.       # 'stabilizer' which is congruent $-1$ modulo 'N2' and congruent $+1$
  723.       # modulo 'No'.
  724.     
  725.       # The base is constructed as follows\:
  726.       #
  727.       # Until all classes are touched:
  728.       # 1. Take a point 'pnt' (in 'zumb').
  729.       # 2. Choose a maximal linear independent set 'pnts' in the equivalence
  730.       #    class of 'pnt' (the intersection of the class with 'zumb').
  731.       # 3. Take the 'stabilizer'--orbits of 'pnts' as base elements;
  732.       #    remove the touched equivalence classes.
  733.       # 4. For the representatives 'rep' in 'supergroup'\:
  734.       #    If 'rep' maps 'pnt' to an equivalence class that was not yet
  735.       #    touched, take the 'stabilizer'--orbits of the images of 'pnts'
  736.       #    under 'rep' as base elements;
  737.       #    remove the touched equivalence classes.
  738.     
  739.       # Compute nontriv. representatives of 'supergroup' over 'stabilizer'
  740.       super:= Difference( supergroup, stabilizer );
  741.       supergroup:= [];
  742.       while super <> [] do
  743.         pnt:= super[1];
  744.         Add( supergroup, pnt );
  745.         SubtractSet( super, List( stabilizer, x -> ( x * pnt ) mod N ) );
  746.       od;
  747.     
  748.       N2:= 1; No:= N;
  749.       while No mod 2 = 0 do
  750.         N2:= 2 * N2; No:= No / 2;
  751.       od;
  752.     
  753.       H1:= [];    # will be the subgroup of 'stabilizer' that acts fixed point
  754.                   # freely on the set of equivalence classes
  755.     
  756.       a:= 0;
  757.       for k in stabilizer do
  758.         if k mod 4 = 1 then
  759.           Add( H1, k );
  760.         elif ( k -1 ) mod No = 0
  761.              and ( ( k + 1 ) mod N2 = 0 or ( k + 1 - N2/2 ) mod N2 = 0 ) then
  762.           a:= k;
  763.         fi;
  764.       od;
  765.       if a = 0 then H1:= stabilizer; fi;
  766.     
  767.       orbits:= [];
  768.     
  769.       while zumb <> [] do
  770.         neworbits:= [];
  771.         pnt:= zumb[1];
  772.         d:= 1;
  773.         ord:= N / GcdInt( N, pnt );
  774.         for i in primes do
  775.           if ord mod i^2 = 0 then d:= d * i; fi;
  776.         od;
  777.     
  778.         if ( a = 0 ) or ( ord mod 8 = 0 ) then
  779.     
  780.           # the orbit of 'H1' cannot be a fixed point of 'a'
  781.           for k in [ 0 .. d-1 ] do
  782.             ppnt:= pnt + k * N/d;
  783.             if ppnt in zumb then
  784.               orb:= List( stabilizer, x -> ( ppnt * x ) mod N );
  785.               Add( neworbits, orb );
  786.             fi;
  787.           od;
  788.     
  789.         elif ord mod 4 = 0 then
  790.     
  791.           # 'a' maps each point in the orbit of 'H1' to its inverse
  792.           # (ignore all points)
  793.           orb:= List( stabilizer, x -> ( pnt * x ) mod N );
  794.     
  795.         else
  796.     
  797.           # the orbit of 'H1' is pointwise fixed by 'a'
  798.           for k in [ 0 .. d-1 ] do
  799.             ppnt:= pnt + k * N/d;
  800.             if ppnt in zumb then
  801.               orb:= List( H1, x -> ( ppnt * x ) mod N );
  802.               Add( neworbits, orb );
  803.             fi;
  804.           od;
  805.     
  806.         fi;
  807.     
  808.         for pnt in orb do        # take any of the latest orbits
  809.     
  810.           # remove the equivalence class of 'pnt'
  811.           SubtractSet( zumb, List( [ 0 .. d-1 ], k -> ( pnt+k*N/d ) mod N ) );
  812.         od;
  813.     
  814.         Append( orbits, neworbits );
  815.     
  816.         # use 'supergroup'\:
  817.         # Is there a point in 'zumb' that is not equivalent to
  818.         # '( pnt * rep ) mod N' ?
  819.         # (Note that the factor group 'supergroup / stabilizer' acts on the
  820.         # set of unions of orbits with equivalent elements.)
  821.     
  822.         for rep in supergroup do
  823.     
  824.           # is there an 'x' in 'zumb' that is equivalent to 'pnt * rep' ?
  825.           if ForAny( zumb, x -> ( ( x - pnt * rep ) * d ) mod N = 0 ) then
  826.             Append( orbits, List( neworbits, x->List(x,y->(y*rep) mod N) ) );
  827.             for ppnt in orbits[ Length( orbits ) ] do
  828.               SubtractSet( zumb, List( [ 0..d-1 ], k->(ppnt+k*N/d) mod N ) );
  829.             od;
  830.           fi;
  831.         od;
  832.     
  833.       od;
  834.     fi;
  835.     
  836.     return orbits;
  837.     end;
  838.  
  839.  
  840. #############################################################################
  841. ##
  842. #F  NumberFieldOps.Coefficents( <F>, <z> ) . . coefficients of an alg. number
  843. ##
  844. ##  returns the coefficient vector of the cyclotomic <cyc> relative to
  845. ##  '<F>.base' for a number field <F>.
  846. ##
  847. ##  We have 'Coefficients( <F>, <z> ) * <F>.base = <z>'.
  848. ##
  849. NumberFieldOps.Coefficients := function ( F, z )
  850.     local coeffs;
  851.     
  852.     if IsBound( F.base ) then
  853.    
  854.       if IsBound( F.coeffslist ) then 
  855.  
  856.         # The information about the Lenstra base coefficients is stored
  857.         # in 'F.coeffslist'.
  858.         coeffs:= Sublist( CoeffsCyc( z, NofCyc( F.generators ) ),
  859.                           F.coeffslist );
  860.  
  861.         if IsBound( F.coeffsmat ) then
  862.  
  863.           # basechange from Lenstra base to another one
  864.           coeffs:= coeffs * F.coeffsmat;
  865.         fi;
  866.  
  867.       else
  868.  
  869.         Error( "no <F>.coeffslist" );
  870.  
  871.       fi;
  872.     
  873.       # change to the base given by 'F.base', if necessary\:
  874.       if IsBound( F.basechangemat ) then
  875.         coeffs:= coeffs * F.basechangemat;
  876.       fi;
  877.     
  878.     else
  879.     
  880.       Error( "no base stored, no 'Coefficients' function for that field" );
  881.     
  882.     fi;
  883.     
  884.     return coeffs;
  885.     end;
  886.  
  887.  
  888. #############################################################################
  889. ##
  890. #F  NumberField( <gens> ) . . . . . . . . . . .  create a number field record
  891. #F  NumberField( <n>, <stab> )
  892. #F  NumberField( <subfield>, <poly> )
  893. #F  NumberField( <subfield>, <base> )
  894. ##
  895. ##  number field generated by the elements of the list <gens>,
  896. ##
  897. ##  fixed field of the group generated by <stab> (prime residues modulo <n>)
  898. ##  in the cyclotomic field 'CF( <n> )',
  899. ##
  900. NF := function ( arg )
  901.     local i, j, k, F, subfield, xtension, d, gens, N, stabilizer, gen, base,
  902.           lenst, lenstset, zumb, newlenst, image, a, H1, No, N2, NN, aut,
  903.           pos, p, coeffslist, root, m, C, val, F_base, l;
  904.     
  905.     if Length( arg ) = 1 then                  # NumberField( gens )
  906.       gens:= arg[1];
  907.       N:= NofCyc( gens );                      # envelopping cyclotomic field
  908.       if N = 1 then
  909.         return Rationals;
  910.       elif N = 4 then
  911.         return ShallowCopy( GaussianRationals );
  912.       fi;
  913.     
  914.       stabilizer:= PrimeResidues( N );
  915.       for gen in gens do
  916.         stabilizer:= Filtered( stabilizer, x -> GaloisCyc( gen, x ) = gen );
  917.       od;
  918.     
  919.       F:= NF( N, stabilizer );
  920.       F.generators:= gens;
  921.     
  922.       return F;
  923.     
  924.     elif Length( arg ) = 2 then
  925.     
  926.       # NF( N, stabilizer )   or
  927.       # NF( subfield, base )  or
  928.       # NF( subfield, poly )
  929.     
  930.       if IsInt( arg[1] ) and arg[1] <> 0 and IsList( arg[2] ) then
  931.     
  932.         # NF( N, stabilizer )
  933.     
  934.         N:= arg[1];
  935.         if N mod 2 = 0 and ( N / 2 ) mod 2 <> 0 then
  936.           N:= N / 2;
  937.         fi;
  938.         stabilizer:= arg[2];
  939.     
  940.         if N <= 2 then return Rationals; fi;
  941.     
  942.         d:= Phi(N) / Length( stabilizer );
  943.     
  944.         # reduce the pair '( N, stabilizer )' such that afterwards 'N'
  945.         # describes the envelopping cyclotomic field of the required field;
  946.     
  947.         gens:= GeneratorsPrimeResidues( N );
  948.         NN:= 1;
  949.         if gens.primes[1] = 2 then
  950.     
  951.           if gens.exponents[1] < 3 then
  952.             if not gens.generators[1] in stabilizer then
  953.               NN:= NN * 4;
  954.             fi;
  955.     
  956.           else
  957.     
  958.             # the only case where 'gens.generators[i]' is a list;
  959.             # it contains the generators corresponding to '**' and '*5';
  960.             # the first one is irrelevant for the envelopping cyclotomic
  961.             # field, except if also the other generator is contained.
  962.             if gens.generators[1][2] in stabilizer then
  963.               if not gens.generators[1][1] in stabilizer then
  964.                 NN:= NN * 4;
  965.               fi;
  966.             else
  967.               NN:= NN * 4;
  968.               aut:= gens.generators[1][2];
  969.               while not aut in stabilizer do
  970.                 aut:= ( aut ^ 2 ) mod N;
  971.                 NN:= NN * 2;
  972.               od;
  973.     
  974.             fi;
  975.           fi;
  976.           pos:= 2;
  977.         else
  978.           pos:= 1;
  979.         fi;
  980.     
  981.         for i in [ pos .. Length( gens.primes ) ] do
  982.           p:= gens.primes[i];
  983.           if not gens.generators[i] in stabilizer then
  984.             NN:= NN * p;
  985.             aut:= ( gens.generators[i] ^ ( p - 1 ) ) mod N;
  986.             while not aut in stabilizer do
  987.               aut:= ( aut ^ p ) mod N;
  988.               NN:= NN * p;
  989.             od;
  990.           fi;
  991.         od;
  992.     
  993.         N:= NN;
  994.     
  995.         stabilizer:= Set( List( stabilizer, x -> x mod N ) );
  996.     
  997.         if stabilizer = [ 1 ] then return CF( N ); fi;
  998.     
  999.         # compute the standard Lenstra base and 'F.coeffslist'\:
  1000.         # If 'stabilizer' acts fixed point freely on the equivalence classes
  1001.         # we must change from the Zumbroich base to a 'stabilizer'--normal
  1002.         # base and afterwards choose coefficients with respect to that base.
  1003.         # In the case of fixed points, only the subgroup 'H1' of index 2 in
  1004.         # stabilizer acts fixed point freely; we change to a 'H1'--normal
  1005.         # base and afterwards choose coefficients.
  1006.     
  1007.         N2:= 1; No:= N;
  1008.         while No mod 2 = 0 do
  1009.           N2:= 2 * N2; No:= No / 2;
  1010.         od;
  1011.     
  1012.         H1:= [];    # will be the subgroup of 'stabilizer' that acts fixed
  1013.                     # point freely on the set of equivalence classes
  1014.     
  1015.         a:= 0;
  1016.         for k in stabilizer do
  1017.           if k mod 4 = 1 then
  1018.             Add( H1, k );
  1019.           elif ( k -1 ) mod No = 0
  1020.                and ( ( k+1 ) mod N2 = 0 or ( k+1-N2/2 ) mod N2 = 0 ) then
  1021.             a:= k;
  1022.           fi;
  1023.         od;
  1024.         if a = 0 then H1:= stabilizer; fi;
  1025.     
  1026.         zumb := ZumbroichBase( N, 1 );
  1027.         lenst:= LenstraBase( N, H1, stabilizer );
  1028.         
  1029.         # We want 'Sublist(CoeffsCyc(z,N),F.coeffslist) = Coefficients(F,z)'
  1030.         # ( and   'Coefficients( F, z ) * F.base = z' )
  1031.         # with respect to the standard Lenstra base.
  1032.     
  1033.         if H1 <> stabilizer then
  1034.     
  1035.           # let 'a' act on 'lenst' to get the right base
  1036.           newlenst:= [];
  1037.           lenstset:= List( lenst, Set );
  1038.           for i in [ 1 .. Length( lenst ) ] do
  1039.             if ( lenst[i][1] * ( a - 1 ) ) mod N = 0 then
  1040.     
  1041.               # pointwise fixed
  1042.               Add( newlenst, lenst[i] );
  1043.     
  1044.             elif ( lenst[i][1] * ( a - 1 ) - N/2 ) mod N <> 0 then
  1045.     
  1046.               # 'a' joins two 'H1'--orbits
  1047.               image:= Set( List( lenst[i], x -> ( x * a ) mod N ) );
  1048.     
  1049.               # *Note* that the elements of 'image' need not be in an element
  1050.               # of 'lenst', only a member of the equivalence class must be
  1051.               # contained;
  1052.               if Position( lenstset, image ) >= i then
  1053.                 Add( newlenst, Concatenation( lenst[i], image ) );
  1054.               fi;
  1055.             fi;
  1056.           od;
  1057.     
  1058.           lenst:= newlenst;
  1059.     
  1060.         fi;
  1061.     
  1062.         coeffslist:= List( lenst, x -> x[1] + 1 );
  1063.         base:=       List( lenst, x -> Sum( List( x, y -> E(N)^y ) ) );
  1064.     
  1065.         return rec( isDomain                := true,
  1066.                     isField                 := true,
  1067.     
  1068.                     char                    := 0,
  1069.                     degree                  := d,
  1070.                     generators              := base,
  1071.     
  1072.                     zero                    := 0,
  1073.                     one                     := 1,
  1074.     
  1075.                     size                    := "infinity",
  1076.                     isFinite                := false,
  1077.     
  1078.                     stabilizer              := stabilizer,
  1079.                     field                   := 0,
  1080.                     dimension               := d,
  1081.                     base                    := base,
  1082.                     isIntegralBase          := true,
  1083.                     coeffslist              := coeffslist,
  1084.     
  1085.                     operations              := NumberFieldOps );
  1086.     
  1087.       elif IsVector( arg[2] ) then
  1088.     
  1089.         subfield:= arg[1];
  1090.         xtension:= arg[2];
  1091.     
  1092.         if ( subfield = 0 and ForAll( xtension, IsRat ) ) or
  1093.            ( IsNumberField( subfield ) and ForAll( xtension,
  1094.                                                    x -> x in subfield ) ) then
  1095.     
  1096.           # NF( subfield, poly )
  1097.     
  1098.           if Length( xtension ) > 3 then
  1099.             Error("NF(<subfield>,<poly>) for polynomial of degree at most 2");
  1100.           fi;
  1101.           if Length( xtension ) = 2 then   # linear polynomial
  1102.  
  1103.             if subfield = 0 then
  1104.               return Rationals;
  1105.             else
  1106.               return NF( subfield.generators );
  1107.             fi;
  1108.  
  1109.           else
  1110.             
  1111.             # The roots of 'a*x^2 + b*x + c' are
  1112.             # $\frac{ -b \pm \sqrt{ b^2 - 4ac } }{2a}$.
  1113.     
  1114.             root:= ( ER( xtension[2]^2 - 4*xtension[1]*xtension[3] )
  1115.                            -xtension[2] ) / 2*xtension[3];
  1116.     
  1117.             if subfield = 0 then
  1118.               return NF( [ root ] );
  1119.             elif root in subfield then
  1120.               return NF( subfield, subfield.base );
  1121.             else
  1122.               return NF( subfield, [ 1, root ] );
  1123.             fi;
  1124.  
  1125.           fi;
  1126.  
  1127.         else
  1128.     
  1129.           # 'NF( 0, base )' or 'NF( subfield, base )'
  1130.           # where 'base' at least contains a vector space base
  1131.           if subfield = 0 or subfield = Rationals then
  1132.             F:= NF( xtension );
  1133.           else
  1134.             F:= NF( Union( subfield.generators, xtension ) );
  1135.             if IsCyclotomicField( F ) then
  1136.               return CF( subfield, xtension );
  1137.             fi;
  1138.    
  1139.             # general case\:\ extension of number field; do not
  1140.             #                 ask if <subfield> is a cyclotomic field 
  1141.  
  1142.             # Let $(v_1, \ldots, v_m)$ denote 'subfield.base' and
  1143.             #     $(w_1, \ldots, w_k)$ denote 'F.base';
  1144.             # Define $u_{i+m(j-1)} = v_i w_j$.  Then $(u_l; 1\leq l\leq mk)$
  1145.             # is a $Q$--base of 'F'.  First change from the Lenstra base to
  1146.             # this base; the matrix is 'C'\:
  1147.  
  1148.             F.dimension    := F.degree / subfield.degree;
  1149.             F.field        := subfield;
  1150.             F_base         := NormalBaseNumberField( F );
  1151.  
  1152.             m:= Length( subfield.base );
  1153.             k:= Length( F_base );
  1154.             N:= NofCyc( F_base );
  1155.             C:= [];
  1156.             for j in F_base do
  1157.               for i in subfield.base do
  1158.                 Add( C, F.operations.Coefficients( F, i*j ) );
  1159.                 # (These are the Lenstra base coefficients!)
  1160.               od;
  1161.             od;
  1162.             C:= C^(-1);
  1163.  
  1164.             # Let $(c_1, \ldots, c_{mk})$ denote the coefficients with respect
  1165.             # to the new base.  To achieve '<coeffs> \* F_base = <z>' we have
  1166.             # to take $\sum_{i=1}^m c_{i+m(j-1)} v_i$ as $j$--th coefficient\:
  1167.  
  1168.             F.coeffsmat:= [];
  1169.             for i in [ 1 .. Length( C ) ] do     # for all rows
  1170.               F.coeffsmat[i]:= [];
  1171.               for j in [ 1 .. k ] do
  1172.                 val:= 0;
  1173.                 for l in [ 1 .. m ] do
  1174.                   val:= val + C[i][ m*(j-1)+l ] * subfield.base[l];
  1175.                 od;
  1176.                 F.coeffsmat[i][j]:= val;
  1177.               od;
  1178.             od;
  1179.  
  1180.             # Multiplication of a Lenstra base coefficient vector with
  1181.             # 'F.coeffsmat' means first changing to the base of products
  1182.             # $v_i w_j$ and then summation over the parts of the $v_i$.
  1183.  
  1184.             F.base         := F_base;
  1185.             F.isNormalBase := true;
  1186.             Unbind( F.isIntegralBase );
  1187.             Unbind( F.name );
  1188.  
  1189.           fi;
  1190.     
  1191.           # If 'xtension' just contains a base but is no base, do nothing;
  1192.           # 'xtension' is a base if and only if the basechange is regular:
  1193.           if Length( xtension ) = Length( F.base ) and xtension <> F.base then
  1194.  
  1195.             F.basechangemat:= List( xtension,
  1196.                                   x -> F.operations.Coefficients(F,x) )^(-1);
  1197.             F.base         := Copy( xtension );
  1198.             Unbind( F.isNormalBase );
  1199.             Unbind( F.name );
  1200.  
  1201.           fi;
  1202.     
  1203.         fi;
  1204.     
  1205.       else
  1206.         Error("usage: NF(<gens>),NF(<N>,<stab>),NF(<S>,<base>),NF(<S>,<poly>)");
  1207.       fi;
  1208.     
  1209.     else
  1210.       Error("<subfield> must be zero or a number field");
  1211.     fi;
  1212.     
  1213.     # return the number field record
  1214.     return F;
  1215.     end;
  1216.  
  1217. NumberField := NF;
  1218.  
  1219.  
  1220. #############################################################################
  1221. ##
  1222. #V  CyclotomicFieldOps  . . . . . . . operations record for cyclotomic fields
  1223. ##
  1224. CyclotomicFieldOps := Copy( NumberFieldOps );
  1225.  
  1226.  
  1227. #############################################################################
  1228. ##
  1229. #F  CyclotomicFieldOps.\in( <z>, <F> ) . . . . . .   test if <z> lies in <F>
  1230. ##
  1231. CyclotomicFieldOps.\in := function ( z, F )
  1232.     return IsCyc( z ) and NofCyc( F.generators ) mod NofCyc( z ) = 0; end;
  1233.  
  1234.  
  1235. #############################################################################
  1236. ##
  1237. #F  CyclotomicFieldOps.Print( <F> ) . . . . . . . .  print a cyclotomic field
  1238. ##
  1239. CyclotomicFieldOps.Print := function( F )
  1240.     if IsBound( F.name ) then
  1241.       Print( F.name );
  1242.     elif IsBound( F.basechangemat ) then
  1243.       Print( "CF( ", F.field, ",", F.base, ")" );
  1244.     else
  1245.       Print( "CF(", NofCyc( F.generators ), ")" );
  1246.       if F.field <> 0 and F.field <> Rationals then
  1247.         Print( "/", F.field );
  1248.       fi;
  1249.     fi;
  1250.     end;
  1251.  
  1252.  
  1253. #############################################################################
  1254. ##
  1255. #F  CyclotomicFieldOps.Conjugates( <F>, <z> ) . . .  conjugates of <z> in <F>
  1256. ##
  1257. CyclotomicFieldOps.Conjugates := function ( F, z )
  1258.     local i, conj;
  1259.     if not z in F then Error( "<z> must lie in <F>" ); fi;
  1260.     if F.field = 0 or F.field = Rationals then
  1261.       conj:= List( PrimeResidues( NofCyc( F.generators ) ),
  1262.                    i -> GaloisCyc( z, i ) );
  1263.     else
  1264.       conj:= [];
  1265.       for i in PrimeResidues( NofCyc( F.generators ) ) do
  1266.         if ForAll( F.field.generators, x -> GaloisCyc( x, i ) = x ) then
  1267.           Add( conj, GaloisCyc( z, i ) );
  1268.         fi;
  1269.       od;
  1270.     fi;
  1271.     return conj;
  1272.     end;
  1273.  
  1274.  
  1275. #############################################################################
  1276. ##
  1277. #F  CyclotomicFieldOps.Norm( <F>, <z> ) . . . . . . . .  norm of a cyclotomic
  1278. ##
  1279. CyclotomicFieldOps.Norm := function ( F, z )
  1280.     local i, result;
  1281.     result:= 1;
  1282.     if F.field = 0 or F.field = Rationals then
  1283.       for i in PrimeResidues( NofCyc( F.generators ) ) do
  1284.         result:= result * GaloisCyc( z, i );
  1285.       od;
  1286.     else
  1287.       for i in PrimeResidues( NofCyc( F.generators ) ) do
  1288.         if ForAll( F.field.generators, x -> GaloisCyc( x, i ) = x ) then
  1289.           result:= result * GaloisCyc( z, i );
  1290.         fi;
  1291.       od;
  1292.     fi;
  1293.     return result;
  1294.     end;
  1295.  
  1296.  
  1297. #############################################################################
  1298. ##
  1299. #F  CyclotomicFieldOps.Trace( <F>, <z> )  . . . . . . . trace of a cyclotomic
  1300. ##
  1301. CyclotomicFieldOps.Trace := function ( F, z )
  1302.     local i, result;
  1303.     result:= 0;
  1304.     if F.field = 0 or F.field = Rationals then
  1305.       for i in PrimeResidues( NofCyc( F.generators ) ) do
  1306.         result:= result + GaloisCyc( z, i );
  1307.       od;
  1308.     else
  1309.       for i in PrimeResidues( NofCyc( F.generators ) ) do
  1310.         if ForAll( F.field.generators, x -> GaloisCyc( x, i ) = x ) then
  1311.           result:= result + GaloisCyc( z, i );
  1312.         fi;
  1313.       od;
  1314.     fi;
  1315.     return result;
  1316.     end;
  1317.  
  1318.  
  1319. #############################################################################
  1320. ##
  1321. #F  CyclotomicFieldOps.Coefficents( <F>, <z> ) . coefficients of a cyclotomic
  1322. ##
  1323. ##  returns the coefficient vector of the cyclotomic <z> relative to
  1324. ##  '<F>.base' for a proper cyclotomic field <F>.
  1325. ##
  1326. ##  We have 'Coefficients( <F>, <z> ) * <F>.base = <z>'.
  1327. ##
  1328. CyclotomicFieldOps.Coefficients := function ( F, z )
  1329.     local N, m, j, k, p, NN, coeffs, value, zumb, Em;
  1330.     
  1331.     N:= NofCyc( F.generators );
  1332.     
  1333.     if not z in F then
  1334.       Error( "no representation of <z> in ",Ordinal(N)," roots of unity" );
  1335.     fi;
  1336.     
  1337.     # get the Zumbroich base representation of <z> in 'N'-th roots
  1338.     coeffs:= CoeffsCyc( z, N );
  1339.     
  1340.     if F.field = 0 or F.field = Rationals then
  1341.     
  1342.       # the Zumbroich base coefficients
  1343.       coeffs:= Sublist( coeffs, F.zumbroichbase + 1 );
  1344.     
  1345.     elif IsCyclotomicField( F.field ) then
  1346.     
  1347.       # direct computation of Zumbroich base coefficients
  1348.       # (without usage of 'F.coeffslist' or 'F.basechangemat')
  1349.       m:= NofCyc( F.field.generators );
  1350.       zumb:= F.field.zumbroichbase;
  1351.       NN:= N/m;
  1352.       Em:= E(m);
  1353.       coeffs:= List( F.zumbroichbase,
  1354.                      j->Sum( zumb, k->coeffs[ ((k*NN+j) mod N )+1 ]*Em^k ) );
  1355.     else
  1356.     
  1357.       # 'F.field' is only a number field.
  1358.       # The necessary information is stored in 'F.coeffsmat'.
  1359.       coeffs:= Sublist( coeffs, F.zumbroichbase + 1 );
  1360.     
  1361.     fi;
  1362.     
  1363.     # basechange from Lenstra base to another one
  1364.     if IsBound( F.coeffsmat ) then
  1365.       coeffs:= coeffs * F.coeffsmat;
  1366.     fi;
  1367.  
  1368.     # change the base according to the value of <F>.base
  1369.     if IsBound( F.basechangemat ) then
  1370.       coeffs:= coeffs * F.basechangemat;
  1371.     fi;
  1372.     
  1373.     return coeffs;
  1374.     end;
  1375.  
  1376.  
  1377. #############################################################################
  1378. ##
  1379. #F  CyclotomicField( <n> )  . . . . . . . .  create a cyclotomic field record
  1380. #F  CyclotomicField( <gens> )
  1381. #F  CyclotomicField( <subfield>, <n> )
  1382. #F  CyclotomicField( <subfield>, <base> )
  1383. ##
  1384. ##  Note: unorthodox generators and base if n is even and n/2 odd
  1385. ##
  1386. CF := function ( arg )
  1387.     local i, j, k, l, m, C, F, subfield, xtension, d, zumb, N, val;
  1388.     
  1389.     # if necessary split the arguments
  1390.     if Length( arg ) = 1  then
  1391.     
  1392.       # CF( N ) or CF( <gens> )
  1393.       subfield:= 0;
  1394.       xtension:= arg[1];
  1395.       if IsVector( xtension ) then xtension:= NofCyc( xtension ); fi;
  1396.  
  1397.       if xtension = 1 then
  1398.         return Rationals;
  1399.       elif xtension = 4 then
  1400.         return ShallowCopy( GaussianRationals );
  1401.       fi;
  1402.     
  1403.     elif Length( arg ) = 2  then
  1404.       subfield:= arg[1];
  1405.       xtension:= arg[2];
  1406.     else
  1407.       Error("usage: CF( <n> ) or CF( <subfield>, <extension> )");
  1408.     fi;
  1409.     
  1410.     # the 'NofCyc' of the extension is given
  1411.     if IsInt( xtension )  then
  1412.     
  1413.       # the subfield is given by '0' or 'Rationals' denoting the prime field
  1414.       if subfield = Rationals or subfield = 0 then
  1415.     
  1416.         # CF( N ) or CF( 0, N ) or CF( Rationals, N )
  1417.         if xtension > 1 then
  1418.           d:= Phi( xtension );
  1419.         else
  1420.           d:= 1;
  1421.         fi;
  1422.  
  1423.         zumb:= ZumbroichBase( xtension, 1 );
  1424.     
  1425.         return rec( isDomain                := true,
  1426.                     isField                 := true,
  1427.                     isCyclotomicField       := true,
  1428.     
  1429.                     char                    := 0,
  1430.                     degree                  := d,
  1431.                     generators              := [ E( xtension ) ],
  1432.     
  1433.                     zero                    := 0,
  1434.                     one                     := 1,
  1435.     
  1436.                     size                    := "infinity",
  1437.                     isFinite                := false,
  1438.     
  1439.                     field                   := 0,
  1440.                     dimension               := d,
  1441.                     base                    := List( zumb,
  1442.                                                      x -> E( xtension )^x ),
  1443.                     isIntegralBase          := true,
  1444.                     zumbroichbase           := zumb,
  1445.                     stabilizer              := [ 1 ],
  1446.     
  1447.                     operations              := CyclotomicFieldOps );
  1448.     
  1449.       elif IsNumberField( subfield ) then
  1450.     
  1451.         # CF( subfield, N )
  1452.         if not xtension mod NofCyc( subfield.generators ) = 0 then
  1453.           Error( "<subfield> is not contained in CF( <xtension> )" );
  1454.         fi;
  1455.     
  1456.         F:= CF( xtension );
  1457.         F.field:= subfield;
  1458.         F.dimension:= F.dimension / subfield.dimension;
  1459.     
  1460.         if IsCyclotomicField( subfield ) then
  1461.           zumb:= ZumbroichBase( xtension, NofCyc( subfield.generators ) );
  1462.           F.base:= List( zumb, x -> E( xtension )^x );
  1463.           F.zumbroichbase:= zumb;
  1464.         else
  1465.           F.base:= NormalBaseNumberField( F );
  1466.           Unbind( F.isIntegralBase );
  1467.           F.isNormalBase:= true;
  1468.           F.zumbroichbase:= ZumbroichBase( xtension, 1 );
  1469.           zumb:= F.zumbroichbase + 1;
  1470.  
  1471.           # Let $(v_1, \ldots, v_m)$ denote 'subfield.base' and
  1472.           #     $(w_1, \ldots, w_k)$ denote 'F.base';
  1473.           # Define $u_{i+m(j-1)} = v_i w_j$.  Then $(u_l; 1\leq l\leq mk)$ is
  1474.           # a $Q$--base of 'F'.  First change from the Zumbroich base to
  1475.           # this base; the matrix is 'C'\:
  1476.  
  1477.           m:= Length( subfield.base );
  1478.           k:= Length( F.base );
  1479.           C:= [];
  1480.           for j in F.base do
  1481.             for i in subfield.base do
  1482.               Add( C, Sublist( CoeffsCyc( i*j, xtension ), zumb ) );
  1483.             od;
  1484.           od;
  1485.           C:= C^(-1);
  1486.  
  1487.           # Let $(c_1, \ldots, c_{mk})$ denote the coefficients with respect
  1488.           # to the new base.  To achieve '<coeffs> \* F.base = <z>' we have
  1489.           # to take $\sum_{i=1}^m c_{i+m(j-1)} v_i$ as $j$--th coefficient\:
  1490.  
  1491.           F.coeffsmat:= [];
  1492.           for i in [ 1 .. Length( C ) ] do     # for all rows
  1493.             F.coeffsmat[i]:= [];
  1494.             for j in [ 1 .. k ] do
  1495.               val:= 0;
  1496.               for l in [ 1 .. m ] do
  1497.                 val:= val + C[i][ m*(j-1)+l ] * subfield.base[l];
  1498.               od;
  1499.               F.coeffsmat[i][j]:= val;
  1500.             od;
  1501.           od;
  1502.  
  1503.         fi;
  1504.         return F;
  1505.     
  1506.       else
  1507.         Error("<subfield> must be zero, Rationals or a number field");
  1508.       fi;
  1509.     
  1510.     elif IsVector( xtension ) then
  1511.     
  1512.       # first create the right field extension, afterwards change the base
  1513.     
  1514.       if subfield = 0 or subfield = Rationals then
  1515.     
  1516.         # CF( 0, base )
  1517.         F:= CF( NofCyc( xtension ) );
  1518.     
  1519.       elif IsNumberField( subfield ) then
  1520.     
  1521.         # CF( F, base )
  1522.         F:= CF( subfield, NofCyc( Union( subfield.generators, xtension ) ) );
  1523.     
  1524.       else
  1525.         Error("<subfield> must be zero, Rationals or a number field");
  1526.       fi;
  1527.     
  1528.       # If 'xtension' just contains a base but is no base, do nothing;
  1529.       # 'xtension' is a base if and only if the basechange is regular:
  1530.       if Length( xtension ) = Length( F.base ) and xtension <> F.base then
  1531.         F.basechangemat:= List( xtension,
  1532.                               x -> F.operations.Coefficients( F, x ) )^(-1);
  1533.         F.base         := Copy( xtension );
  1534.         Unbind( F.isNormalBase );
  1535.         Unbind( F.name );
  1536.       fi;
  1537.       if subfield = 0 or subfield = Rationals then
  1538.         F.generators   := Copy( xtension );
  1539.       else
  1540.         F.generators   := Union( subfield.generators, xtension );
  1541.       fi;
  1542.       return F;
  1543.     
  1544.     # otherwise we don't know how to handle the parameters
  1545.     else
  1546.       Error("usage: CF(<n>) or CF(<subfield>,<n>) or CF(<subfield>,<base>)");
  1547.     fi;
  1548.     end;
  1549.  
  1550. CyclotomicField := CF;
  1551.  
  1552.  
  1553. #############################################################################
  1554. ##
  1555. #F  NumberRing( gens )  . . . . . . . . . . . integral ring in a number field
  1556. #F  CyclotomicRing( gens )  . . . . . . . integral ring in a cyclotomic field
  1557. #V  NumberRingOps . . . . . . . . . . . . operations records for number rings
  1558. ##
  1559. NumberRingOps := Copy( RingOps );
  1560.  
  1561. NumberRingOps.\in := function( z, R )
  1562.     return IsCycInt( z )
  1563.            and NofCyc( R.generators ) mod NofCyc( z ) = 0
  1564.            and ForAll( R.stabilizer, x -> GaloisCyc( z, x ) = z );
  1565.     end;
  1566.  
  1567. NumberRingOps.Quotient := function( R, x, y )
  1568.     local q;
  1569.     q:= x/y;
  1570.     if q in R then return x / y; else return false; fi;
  1571.     end;
  1572.  
  1573. NumberRingOps.IsUnit := function( R, z ) return 1 / z in R; end;
  1574. NumberRingOps.Random := R -> List( R.base, x -> Random(Integers) ) * R.base;
  1575.  
  1576. NumberRing := function( gens )
  1577.     local F;
  1578.     if ForAll( gens, IsCycInt ) then
  1579.  
  1580.       if NofCyc( gens ) = 4 then
  1581.         return ShallowCopy( GaussianIntegers );
  1582.       fi;
  1583.  
  1584.       F:= NF( gens );
  1585.       return rec( isDomain       := true,
  1586.                   isRing         := true,
  1587.  
  1588.                   generators     := gens,
  1589.                   zero           := 0,
  1590.                   one            := 1,
  1591.  
  1592.                   size           := "infinity",
  1593.                   isFinite       := false,
  1594.                   isCommutative  := true,
  1595.                   isIntegralRing := true,
  1596.  
  1597.                   base           := F.base,
  1598.                   stabilizer     := F.stabilizer,
  1599.                   rank           := F.dimension,
  1600.  
  1601.                   operations     := NumberRingOps  );
  1602.     else
  1603.       Error( "<gens> must be cyclotomic integers" );
  1604.     fi;
  1605.     end;
  1606.  
  1607. CyclotomicRing := function( gens )
  1608.     local F;
  1609.     if ForAll( gens, IsCycInt ) then
  1610.  
  1611.       if NofCyc( gens ) = 4 then
  1612.         return ShallowCopy( GaussianIntegers );
  1613.       fi;
  1614.  
  1615.       F:= CF( gens );
  1616.       return rec( isDomain       := true,
  1617.                   isRing         := true,
  1618.  
  1619.                   generators     := F.generators,
  1620.                   zero           := 0,
  1621.                   one            := 1,
  1622.  
  1623.                   size           := "infinity",
  1624.                   isFinite       := false,
  1625.                   isCommutative  := true,
  1626.                   isIntegralRing := true,
  1627.  
  1628.                   base           := F.base,
  1629.                   stabilizer     := F.stabilizer,
  1630.                   rank           := F.dimension,
  1631.  
  1632.                   operations     := NumberRingOps  );
  1633.       else
  1634.         Error( "<gens> must be cyclotomic integers" );
  1635.       fi;
  1636.     end;
  1637.  
  1638.  
  1639. #############################################################################
  1640. ##
  1641. #V  Cyclotomics . . . . . . . . . . . . . . . . . . domain of all cyclotomics
  1642. #V  CyclotomicsOps  . . . . . . . . . . . . operations record for this domain
  1643. ##
  1644. CyclotomicsOps              := Copy( FieldElementsOps );
  1645.     CyclotomicsOps.\in         := function( z, F ) return IsCyc( z ); end;
  1646.     CyclotomicsOps.Order        := NumberFieldOps.Order;
  1647.     CyclotomicsOps.Field        := ( gens -> NumberField( gens ) );
  1648.     CyclotomicsOps.DefaultField := ( gens -> CF( NofCyc( gens ) ) );
  1649.     CyclotomicsOps.Ring         := ( gens -> NumberRing( gens ) );
  1650.     CyclotomicsOps.DefaultRing  := ( gens -> CyclotomicRing( gens ) );
  1651.  
  1652. Cyclotomics           := rec( isDomain  := true,
  1653.                               name      := "Cyclotomics",
  1654.                               isFinite  := false,
  1655.                               size      := "infinity",
  1656.                               operations:= CyclotomicsOps );
  1657.  
  1658.  
  1659. #############################################################################
  1660. ##
  1661. #E  Emacs . . . . . . . . . . . . . . . . . . . . . . . local emacs variables
  1662. ##
  1663. ##  Local Variables:
  1664. ##  mode:               outline
  1665. ##  outline-regexp:     "#F\\|#V\\|#E"
  1666. ##  fill-column:        73
  1667. ##  fill-prefix:        "##  "
  1668. ##  eval:               (hide-body)
  1669. ##  End:
  1670. ##
  1671.  
  1672.  
  1673.  
  1674.