home *** CD-ROM | disk | FTP | other *** search
/ Math Solutions 1995 October / Math_Solutions_CD-ROM_Walnut_Creek_October_1995.iso / pc / mac / discrete / lib / cyclotom.g < prev    next >
Encoding:
Text File  |  1993-05-05  |  29.3 KB  |  896 lines
  1. #############################################################################
  2. ##
  3. #A  cyclotom.g                  GAP library                     Thomas Breuer
  4. ##
  5. #A  @(#)$Id: cyclotom.g,v 3.12 1992/02/24 18:00:26 sam Rel $
  6. ##
  7. #Y  Copyright 1990-1992,  Lehrstuhl D fuer Mathematik,  RWTH Aachen,  Germany
  8. ##
  9. ##  This file contains  those  functions  that  mainly  deal  with cyclotomics,
  10. ##  that is sums of rational multiples of roots of unity.
  11. ##
  12. #H  $Log: cyclotom.g,v $
  13. #H  Revision 3.12  1992/02/24  18:00:26  sam
  14. #H  fixed bug in 'NK'
  15. #H
  16. #H  Revision 3.11  1992/02/13  09:41:38  sam
  17. #H  fixed bug in 'NK'
  18. #H
  19. #H  Revision 3.10  1992/01/03  10:53:48  sam
  20. #H  improved 'Quadratic'
  21. #H
  22. #H  Revision 3.9  1991/12/12  11:36:56  sam
  23. #H  fixed bug in 'Quadratic'
  24. #H
  25. #H  Revision 3.8  1991/12/04  09:29:30  sam
  26. #H  inserted a <space> after each function name
  27. #H
  28. #H  Revision 3.7  1991/12/03  17:09:15  sam
  29. #H  added some functions for rings, corrected some bugs
  30. #H
  31. #H  Revision 3.6  1991/11/29  14:48:23  sam
  32. #H  completely revised, domain dependent functions added
  33. #H
  34. #H  Revision 3.5  1991/10/15  07:43:27  sam
  35. #H  little changes in 'Quadratic'
  36. #H
  37. #H  Revision 3.4  1991/10/01  14:29:00  casper
  38. #H  changed 'Gcd' to 'GcdInt'
  39. #H
  40. #H  Revision 3.3  1991/09/05  15:54:23  sam
  41. #H  changed 'LcmList' to 'Lcm'
  42. #H
  43. #H  Revision 3.2  1991/09/05  15:36:45  sam
  44. #H  changed 'PrimitiveRoot' to 'PrimitiveRootMod'
  45. #H
  46. #H  Revision 3.1  1991/09/05  14:53:21  sam
  47. #H  changed 'Quo' to 'QuoInt'
  48. #H
  49. #H  Revision 3.0  1991/09/03  13:58:53  goetz
  50. #H  Initial Revision.
  51. #H
  52. ##
  53.  
  54.  
  55. #############################################################################
  56. ##
  57. #F  IntCyc( <cyc> ) . . . . . . . . . . . .  cyclotomic integer near to <cyc>
  58. ##
  59. IntCyc := function ( x )
  60.     local i, int, n, cfs;
  61.     n:= NofCyc( x );
  62.     cfs:= COEFFSCYC( x );
  63.     int:= 0;
  64.     for i in [ 1 .. n ] do int:= int + Int( cfs[i] ) * E(n)^(i-1); od;
  65.     return int;
  66.     end;
  67.  
  68.  
  69. #############################################################################
  70. ##
  71. #F  RoundCyc( <cyc> ) . . . . . . . . . . .  cyclotomic integer near to <cyc>
  72. ##
  73. RoundCyc := function ( x )
  74.     local i, int, n, cfs;
  75.     n:= NofCyc( x );
  76.     cfs:= COEFFSCYC( x );
  77.     int:= 0;
  78.     for i in [ 1 .. n ]  do
  79.       if cfs[i] < 0 then
  80.         int:= int + Int( cfs[i]-1/2 ) * E(n)^(i-1);
  81.       else
  82.         int:= int + Int( cfs[i]+1/2 ) * E(n)^(i-1);
  83.       fi;
  84.     od;
  85.     return int;
  86.     end;
  87.  
  88.  
  89. #############################################################################
  90. ##
  91. #F  'CoeffsCyc( <z>, <N> )'
  92. ##
  93. ##  If <z> is a cyclotomic that lies in the field of <N>-th roots of unity,
  94. ##  'CoeffsCyc' returns a list of length <N> which is the Zumbroich base 
  95. ##  representation of <z> in that field, i.e. at position 'i' the
  96. ##  coefficient of 'E(N)^(i-1)' is stored.
  97. ##
  98. CoeffsCyc := function( z, N )
  99.     local j, k, p, n, quo, coeffs, newcoeffs;
  100.     coeffs:= COEFFSCYC( z );     # the internal function:
  101.                                  # returns 'CoeffsCyc( z, NofCyc( z ) )'
  102.     quo:= N / NofCyc( z );
  103.     if not IsInt( quo ) then
  104.       Error( "no representation of <z> in ",Ordinal(N)," roots of unity" );
  105.     fi;
  106.     
  107.     # blow up step by step
  108.     
  109.     for p in FactorsInt( quo ) do
  110.       n:= Length( coeffs );
  111.       newcoeffs:= List( [ 1 .. p*n ], x -> 0 );
  112.       if p = 2 or n mod p = 0 then
  113.     
  114.         # simply blow up the base elements:
  115.         # for a base element 'E(n)^k' we have 'E(n*p)^(k*p)' in the base, too
  116.     
  117.         for k in [ 1 .. Length( coeffs ) ] do
  118.           if coeffs[k] <> 0 then newcoeffs[ (k-1)*p + 1 ]:= coeffs[k]; fi;
  119.         od;
  120.       else
  121.     
  122.         # replace 'E(n)^k' by $ - \sum_{j=1}^p 'E(n*p)^(p*k+j*n)'$
  123.     
  124.         for k in [ 1 .. Length( coeffs ) ] do
  125.           if coeffs[k] <> 0 then
  126.             for j in [ 1 .. p-1 ] do
  127.               newcoeffs[ ( ( (k-1)*p + j*n  ) mod (n*p) ) + 1 ]:= -coeffs[k];
  128.             od;
  129.           fi;
  130.         od;
  131.       fi;
  132.       coeffs:= newcoeffs;
  133.     od;
  134.     return coeffs;
  135.     end;
  136.  
  137.  
  138. #############################################################################
  139. ##
  140. #F  CycList( <coeffs> ) . . . . .  cyclotomic of Zumbroich base coeff. vector
  141. ##
  142. ##  (mainly used to read tables produced by 'ctoc')
  143. ##
  144. ##  *Note*\: 'CycList( COEFFSCYC( <cyc> ) )' = <cyc>, but
  145. ##           'COEFFSCYC( CycList( <coeffs> ))' need not be equal to <coeffs>.
  146. ##
  147. CycList := function( coeffs )
  148.     local i, n;
  149.     n:= Length( coeffs );
  150.     return Sum( [ 1 .. n ], i -> coeffs[i] * E(n)^(i-1) );
  151.     end;
  152.  
  153.  
  154. #############################################################################
  155. ##
  156. #F  Atlas1( <n>, <i> )  . . . . . . . . . . . . . . . utility for EB, ..., EH
  157. ##
  158. ##  calculates the value $\frac{1}{i}\sum{j=1}^{n-1}z_n^{j^i}$ for
  159. ##  $2 \leq i \leq 8$ and $<n> \equiv 1 \pmod{i}$; if $i > 2$, <n> should
  160. ##  be a prime to get sure that the result is well defined; 'ATLAS' returns
  161. ##  the value given above if it is a cyclotomic integer.
  162. ##  (see: Conway et al, ATLAS of finite groups, Oxford University Press 1985;
  163. ##  here: chapter 7, section 10)
  164. ##
  165. Atlas1 := function( n, i )
  166.     local k, kpow, resultlist, atlas;
  167.     if not IsInt( n ) or n < 1 then
  168.       Error("usage: EB(<n>),EC(<n>),..,EH(<n>) with appropriate integer <n>");
  169.     elif n mod i <> 1 then
  170.       Error("<n> not congruent 1 mod ", i );
  171.     fi;
  172.     if n = 1 then return 0; fi;
  173.     atlas:= 0;
  174.     if i mod 2 = 0 then
  175.       for k in [ 1 .. QuoInt( n-1, 2 ) ] do
  176.         atlas:= atlas + 2 * E( n ) ^ ( (k^i) mod n );
  177.       od;
  178.     else
  179.       for k in [ 1 .. QuoInt( n-1, 2 ) ] do
  180.         atlas:= atlas + E(n)^( ( k^i ) mod n ) + E(n)^( n - ( k^i ) mod n );
  181.       od;
  182.       if n mod 2 = 0 then
  183.         atlas:= atlas + E( n ) ^ ( ( ( n / 2 ) ^ i )  mod  n );
  184.       fi;
  185.     fi;
  186.     atlas:= atlas / i;
  187.     if not IsCycInt( atlas ) then
  188.       Error( "result divided by ", i, " is not a cyclotomic integer" );
  189.     fi;
  190.     return atlas;
  191.     end;
  192.  
  193.  
  194. #############################################################################
  195. ##
  196. #F  EB( <n> ), EC( <n> ), \ldots, EH( <n> ) . . .  some ATLAS irrationalities
  197. ##
  198. EB := n -> Atlas1( n, 2 );
  199. EC := n -> Atlas1( n, 3 );
  200. ED := n -> Atlas1( n, 4 );
  201. EE := n -> Atlas1( n, 5 );
  202. EF := n -> Atlas1( n, 6 );
  203. EG := n -> Atlas1( n, 7 );
  204. EH := n -> Atlas1( n, 8 );
  205.  
  206.  
  207. #############################################################################
  208. ##
  209. #F  NK( <n>, <k>, <deriv> ) . . . . . . . . utility for ATLAS irrationalities
  210. ##
  211. ##  'NK( <n>, <k>, <deriv> )' is the $(<deriv>+1)$-th number of multiplicative
  212. ##  order exactly <k> modulo <N>, chosen in the order of preference
  213. ##  \[ 1, -1, 2, -2, 3, -3, 4, -4, \ldots .\]
  214. ##  (see: Conway et al, ATLAS of finite groups, Oxford University Press 1985;
  215. ##  here: chapter 7, section 10)
  216. ##
  217. NK := function( n, k, deriv )
  218.     local i, nk;
  219.     if n <= 2 or ( k in [ 2, 3, 5, 6, 7 ] and Phi( n ) mod k <> 0 )
  220.               or k < 2 or k > 8 then
  221.       Error( "value does not exist" );
  222.     fi;
  223.     if k mod 4 = 0 then
  224.  
  225.       # an automorphism of order 4 exists if 4 divides $p-1$ for an odd
  226.       # prime divisor $p$ of 'n', or if 16 divides 'n';
  227.       # an automorphism of order 8 exists if 8 divides $p-1$ for an odd
  228.       # prime divisor $p$ of 'n', or if 32 divides 'n';
  229.       if ForAll(Set(FactorsInt(n)),x->(x-1) mod k<>0) and n mod (4*k)<>0 then
  230.         Error( "value does not exist" );
  231.       fi;
  232.     fi;
  233.     nk:= 1;
  234.     if k in [ 2, 3, 5, 7 ] then        # for primes
  235.       while true do
  236.         if ( nk ^ k ) mod n = 1 and nk mod n <> 1 then
  237.           if deriv = 0 then return nk; fi;
  238.           deriv:= deriv - 1;
  239.         fi;
  240.         if ( ( (-nk) ^ k ) - 1 ) mod n = 0 and ( -nk -1 ) mod n <> 0 then
  241.           if deriv = 0 then return -nk; fi;
  242.           deriv:= deriv - 1;
  243.         fi;
  244.         nk:= nk + 1;
  245.       od;
  246.     elif k = 4 then
  247.       while true do
  248.         if ( nk ^ 4 ) mod n = 1 and ( nk ^ 2 ) mod n <> 1 then
  249.           if deriv = 0 then return nk; fi;
  250.           deriv:= deriv - 1;
  251.           if deriv = 0 then return -nk; fi;
  252.           deriv:= deriv - 1;
  253.         fi;
  254.         nk:= nk + 1;
  255.       od;
  256.     elif k = 6 then
  257.       while true do
  258.         if (nk^6) mod n = 1 and (nk^2) mod n <> 1 and (nk^3) mod n <> 1 then
  259.           if deriv = 0 then return nk; fi;
  260.           deriv:= deriv - 1;
  261.         fi;
  262.         if (nk^6) mod n=1 and (nk^2) mod n<>1 and (-(nk^3) mod n)+n<>1 then
  263.           if deriv = 0 then return -nk; fi;
  264.           deriv:= deriv - 1;
  265.         fi;
  266.         nk:= nk + 1;
  267.       od;
  268.     elif k = 8 then
  269.       while true do
  270.         if ( nk ^ 8 ) mod n = 1 and ( nk ^ 4 ) mod n <> 1 then
  271.           if deriv = 0 then return nk; fi;
  272.           deriv:= deriv - 1;
  273.           if deriv = 0 then return -nk; fi;
  274.           deriv:= deriv - 1;
  275.         fi;
  276.         nk:= nk + 1;
  277.       od;
  278.     fi; 
  279.     end;
  280.  
  281.  
  282. #############################################################################
  283. ##
  284. #F  Atlas2( <n>, <k>, <deriv> ) . . . . . . utility for ATLAS irrationalities
  285. ##
  286. Atlas2 := function( n, k, deriv )
  287.     local i, e, nk, result;
  288.     
  289.     if not ( IsInt( n ) and IsInt( k ) and IsInt( deriv ) ) then
  290.       Error( "usage: ATLAS irrationalities require integral arguments" );
  291.     fi;
  292.     
  293.     nk:= NK( n, k, deriv );
  294.     e:= E(n);
  295.     result:= 0;
  296.     for i in [ 0 .. k-1 ] do
  297.       result:= result + e^( (nk^i) mod n );
  298.     od;
  299.     return result;
  300.     end;
  301.  
  302.  
  303. #############################################################################
  304. ##
  305. #F  EY(<n>), EY(<n>,<deriv>) . . . . . . .  ATLAS irrationalities $y_n$ resp.
  306. #F                                          $y_n^{<deriv>}$
  307. #F  ... ES(<n>), ES(<n>,<deriv>)              ... $s_n$ resp. $s_n^{<deriv>}$
  308. ##
  309. EY :=function(arg) if   Length(arg)=1 then return Atlas2(arg[1],2,0);
  310.                   elif Length(arg)=2 then return Atlas2(arg[1],2,arg[2]);
  311.                   else Error( "usage: EY(n) resp. EY(n,deriv)" ); fi; end;
  312.  
  313. EX :=function(arg) if   Length(arg)=1 then return Atlas2(arg[1],3,0);
  314.                   elif Length(arg)=2 then return Atlas2(arg[1],3,arg[2]);
  315.                   else Error( "usage: EX(n) resp. EX(n,deriv)" ); fi; end;
  316.  
  317. EW :=function(arg) if   Length(arg)=1 then return Atlas2(arg[1],4,0);
  318.                   elif Length(arg)=2 then return Atlas2(arg[1],4,arg[2]);
  319.                   else Error( "usage: EW(n) resp. EW(n,deriv)" ); fi; end;
  320.  
  321. EV :=function(arg) if   Length(arg)=1 then return Atlas2(arg[1],5,0);
  322.                   elif Length(arg)=2 then return Atlas2(arg[1],5,arg[2]);
  323.                   else Error( "usage: EV(n) resp. EV(n,deriv)" ); fi; end;
  324.  
  325. EU :=function(arg) if   Length(arg)=1 then return Atlas2(arg[1],6,0);
  326.                   elif Length(arg)=2 then return Atlas2(arg[1],6,arg[2]);
  327.                   else Error( "usage: EU(n) resp. EU(n,deriv)" ); fi; end;
  328.  
  329. ET :=function(arg) if   Length(arg)=1 then return Atlas2(arg[1],7,0);
  330.                   elif Length(arg)=2 then return Atlas2(arg[1],7,arg[2]);
  331.                   else Error( "usage: ET(n) resp. ET(n,deriv)" ); fi; end;
  332.  
  333. ES :=function(arg) if   Length(arg)=1 then return Atlas2(arg[1],8,0);
  334.                   elif Length(arg)=2 then return Atlas2(arg[1],8,arg[2]);
  335.                   else Error( "usage: ES(n) resp. ES(n,deriv)" ); fi; end;
  336.  
  337.  
  338. #############################################################################
  339. ##
  340. #F  EM( <n> ), EM( <n>, <deriv> ) . . . . ATLAS irrationality $m_{<n>}$ resp.
  341. ##                                                        $m_{<n>}^{<deriv>}$
  342. EM := function( arg )
  343.     local n;
  344.     n:= arg[1];
  345.     if Length( arg ) = 1 then
  346.       return E(n) - E(n)^(-1);
  347.     elif Length( arg ) = 2 and IsInt( n ) then
  348.       return E(n) - E(n)^( NK( n, 2, arg[2] ) mod n );
  349.     else
  350.       Error( "usage: EM(<n>) resp. EM(<n>,<deriv>)" );
  351.     fi;
  352.     end;
  353.  
  354.  
  355. #############################################################################
  356. ##
  357. #F  EL( <n> ), EL( <n>, <deriv> ) . . . . ATLAS irrationality $l_{<n>}$ resp.
  358. ##                                                        $l_{<n>}^{<deriv>}$
  359. EL := function( arg )
  360.     local n, nk;
  361.     n:= arg[1];
  362.     if Length( arg ) > 2 or not IsInt( n ) then
  363.       Error( "usage: EL( <n> ) resp. EL( <n>, <deriv> )" );
  364.     fi;
  365.     if Length(arg)=1 then
  366.       nk:= NK(n,4,0) mod n;
  367.     else
  368.       nk:= NK(n,4,arg[2]) mod n;
  369.     fi;
  370.     return E(n)-E(n)^nk+E(n)^((nk^2) mod n)-E(n)^((nk^3) mod n);
  371.     end;
  372.  
  373.  
  374. #############################################################################
  375. ##
  376. #F  EK( <n> ), EK( <n>, <deriv> ) . . . . ATLAS irrationality $k_{<n>}$ resp.
  377. ##                                                        $k_{<n>}^{<deriv>}$
  378. EK := function( arg )
  379.     local n, nk;
  380.     n:= arg[1];
  381.     if Length( arg ) > 2 or not IsInt( n ) then
  382.       Error( "usage: EK( <n> ) resp. EK( <n>, <deriv> )" );
  383.     fi;
  384.     if Length(arg)=1 then
  385.       nk:= NK(n,6,0) mod n;
  386.     else
  387.       nk:= NK(n,6,arg[2]) mod n;
  388.     fi;
  389.     return E(n)-E(n)^nk+E(n)^((nk^2) mod n)-E(n)^((nk^3) mod n)+
  390.            E(n)^((nk^4) mod n)-E(n)^((nk^5) mod n);
  391.     end;
  392.  
  393.  
  394. #############################################################################
  395. ##
  396. #F  EJ( <n> ), EJ( <n>, <deriv> ) . . . . ATLAS irrationality $j_{<n>}$ resp.
  397. ##                                                        $j_{<n>}^{<deriv>}$
  398. EJ := function( arg )
  399.     local n, nk;
  400.     n:= arg[1];
  401.     if Length( arg ) > 2 or not IsInt( n ) then
  402.       Error( "usage: EJ( <n> ) resp. EJ( <n>, <deriv> )" );
  403.     fi;
  404.     if Length(arg)=1 then
  405.       nk:= NK(n,8,0) mod n;
  406.     else
  407.       nk:= NK(n,8,arg[2]) mod n;
  408.     fi;
  409.     return E(n)-E(n)^nk+E(n)^((nk^2) mod n)-E(n)^((nk^3) mod n)+
  410.            E(n)^((nk^4) mod n)-E(n)^((nk^5) mod n)+E(n)^((nk^6) mod n)-
  411.            E(n)^((nk^7) mod n);
  412.     end;
  413.  
  414.  
  415. #############################################################################
  416. ##
  417. #F  ER( <n> ) . . . . ATLAS irrationality $r_{<n>}$ (pos. square root of <n>)
  418. ##
  419. ER := function( n )
  420.     local factor;
  421.     if not IsInt( n ) then Error( "argument must be integer valued" ); fi;
  422.     if n = 0 then
  423.       return 0;
  424.     elif n < 0 then
  425.       factor:= E(4);
  426.       n:= -n;
  427.     else
  428.       factor:= 1;
  429.     fi;
  430.     if   n mod 4 = 1 then return         factor * ( 2 * EB(n) + 1 );
  431.     elif n mod 4 = 2 then return ( E(8) - E(8)^3 ) * factor * ER( n / 2 );
  432.     elif n mod 4 = 3 then return -E(4) * factor * ( 2 * EB(n) + 1 );
  433.     else
  434.       return 2 * ER( n / 4 );
  435.     fi;
  436.     end;
  437.  
  438.  
  439. #############################################################################
  440. ##
  441. #F  EI( <n> ) . . . . ATLAS irrationality $i_{<n>}$ (the square root of -<n>)
  442. ##
  443. EI := ( n -> E(4) * ER(n) );
  444.  
  445.  
  446. #############################################################################
  447. ##
  448. #F  StarCyc( <cyc> )  . . . . the unique nontrivial galois conjugate of <cyc>
  449. ##
  450. StarCyc := function( cyc )
  451.     local i, conj;
  452.     conj:= [];
  453.     for i in PrimeResidues( NofCyc( cyc ) ) do
  454.       AddSet( conj, GaloisCyc( cyc, i ) );
  455.     od;
  456.     if Length( conj ) = 2 then
  457.       return Difference( conj, [ cyc ] )[1];
  458.     else
  459.       return false;
  460.     fi;
  461.     end;
  462.  
  463.  
  464. #############################################################################
  465. ##
  466. #F  Quadratic( <cyc> ) . . . . . informations about quadratic irrationalities
  467. ##
  468. ##  If <cyc> is a quadratic irrationality, Quadratic( <cyc> ) calculates the
  469. ##  representation $<cyc> = \frac{ a + b \sqrt{ 'root' } }{'d'}$ and a
  470. ##  (not necessarily shortest) representation by a combination of the \ATLAS\
  471. ##  irrationalities $b_{'root'}, i_{'root'}$ and $r_{'root'}$.
  472. ##  Otherwise 'false' is returned.
  473. ##
  474. ##  1. If the denominator 'd' is 2, necessarily 'root' is congruent 1 mod 4,
  475. ##     and $r_n$, $i_n$ are not possible;
  476. ##     '<cyc> = x + y * EB( root )' with y = b, x = ( a + b ) / 2.
  477. ##  2. If the denominator 'd' is 1, we have the possibilities
  478. ##     $i_n$ for $'root' \< -1$, 'a + b * i' for 'root' = -1, $a + b * r_n$
  479. ##     for $'root' > 0$. Furthermore if 'root' is congruent 1 modulo 4, also
  480. ##     '<cyc> = (a+b) + 2 * b * EB( root )' is possible; the shortest string
  481. ##     of these is taken as value for the field ATLAS.
  482. ##
  483. Quadratic := function( cyc )
  484.     local i, root, a, b, d, list, facts, two_part, ATLAS, ATLAS2;
  485.  
  486.     if not IsCycInt( cyc ) then
  487.       return false;
  488.     elif IsInt( cyc ) then
  489.       return rec( a:= cyc, b:= 0, root:= 1, d:= 1, ATLAS:= String( cyc ) );
  490.     fi;
  491.  
  492.     list:= COEFFSCYC( cyc );
  493.     facts:= FactorsInt( Length( list ) );
  494.     two_part:= 0;
  495.     for i in facts do if i = 2 then two_part:= two_part + 1; fi; od;
  496.  
  497.     if two_part > 3 or
  498.        ( two_part = 2 and Length( facts ) <> Length( Set( facts ) ) + 1 ) or
  499.        ( two_part < 2 and Length( facts ) <> Length( Set( facts ) ) ) then
  500.       return false;
  501.     fi;
  502.  
  503.     facts:= Set( facts );
  504.     if two_part = 0 then
  505.       root:= Length( list );
  506.       if root mod 4 = 3 then root:= -root; fi;
  507.       a:= StarCyc( cyc );
  508.       if a = false then return false; fi;
  509.       a:= cyc + a;                       # trace of 'cyc' over the rationals
  510.       if Length( facts ) mod 2 = 0 then
  511.         b:= 2 * list[2] - a;
  512.       else
  513.         b:= 2 * list[2] + a;
  514.       fi;
  515.       if a mod 2 = 0 and b mod 2 = 0 then
  516.         a:= a / 2;
  517.         b:= b / 2;
  518.         d:= 1;
  519.       else
  520.         d:= 2;
  521.       fi;
  522.     elif two_part = 2 then
  523.       root:= Length( list ) / 4;
  524.       if root = 1 then
  525.         a:= list[1];
  526.         b:= - list[2];
  527.       else
  528.         a:= list[5];
  529.         if Length( facts ) mod 2 = 0 then a:= -a; fi;
  530.         b:= - list[ root + 5 ];
  531.       fi;
  532.       if root mod 4 = 1 then
  533.         root:= -root;
  534.         b:= -b;
  535.       fi;
  536.       d:= 1;
  537.     else        # two_part = 3
  538.       root:= Length( list ) / 4;
  539.       if root = 2 then
  540.         a:= list[1];
  541.         b:= list[2];
  542.         if b = list[4] then
  543.           root:= -2;
  544.         fi;
  545.       else
  546.         a:= list[9];
  547.         if Length( facts ) mod 2 = 0 then a:= -a; fi;
  548.         b:= list[ root / 2 + 9 ];
  549.         if b <> - list[ 3 * root / 2 - 7 ] then
  550.           root:= -root;
  551.         elif ( root / 2 ) mod 4 = 3 then
  552.           b:= -b;
  553.         fi;
  554.       fi;
  555.       d:= 1;
  556.     fi;
  557.     
  558.     if d * cyc <> a + b * ER( root ) then return false; fi;
  559.     if d = 2 then        # necessarily 'root' congruent 1 mod 4,
  560.                 # only $b_{'root'}$ possible
  561.       if a = -b then
  562.         ATLAS:= "";
  563.       else
  564.         ATLAS:= String( ( a + b ) / 2 );
  565.       fi;
  566.       if b = -1 then
  567.         ATLAS:= ConcatenationString( ATLAS, "-" );
  568.       elif b < 0 then
  569.         ATLAS:= ConcatenationString( ATLAS, String( b ) );
  570.       else
  571.         if ATLAS <> "" then
  572.           ATLAS:= ConcatenationString( ATLAS, "+" );
  573.         fi;
  574.         if b <> 1 then
  575.           ATLAS:= ConcatenationString( ATLAS, String( b ) );
  576.         fi;
  577.       fi;
  578.       if root > 0 then
  579.         ATLAS:= ConcatenationString( ATLAS, "b", String( root ) );
  580.       else
  581.         ATLAS:= ConcatenationString( ATLAS, "b", String( -root ) );
  582.       fi;
  583.     else
  584.       if a = 0 then
  585.         ATLAS:= "";
  586.       else
  587.         ATLAS:= String( a );
  588.       fi;
  589.       if a <> 0 and b > 0 then ATLAS:= ConcatenationString( ATLAS, "+" ); fi;
  590.       if b = -1 then
  591.         ATLAS:= ConcatenationString( ATLAS, "-" );
  592.       elif b <> 1 then
  593.         ATLAS:= ConcatenationString( ATLAS, String( b ) );
  594.       fi;
  595.       if root > 0 then
  596.         ATLAS:= ConcatenationString( ATLAS, "r", String( root ) );
  597.       elif root = -1 then
  598.         ATLAS:= ConcatenationString( ATLAS, "i" );
  599.       else
  600.         ATLAS:= ConcatenationString( ATLAS, "i", String( -root ) );
  601.       fi;
  602.       if ( root -1 ) mod 4 = 0 then
  603.         if a = -b then
  604.           ATLAS2:= String( 2 * b );
  605.         else
  606.           ATLAS2:= ConcatenationString( String( a+b ), "+", String( 2*b ) );
  607.         fi;
  608.         if root > 0 then
  609.           ATLAS2:= ConcatenationString( ATLAS2, "b", String( root ) );
  610.         else
  611.           ATLAS2:= ConcatenationString( ATLAS2, "b", String( -root ) );
  612.         fi;
  613.         if LengthString(ATLAS2) < LengthString(ATLAS) then ATLAS:= ATLAS2; fi;
  614.       fi;
  615.     fi;
  616.     return rec( a:= a, b:= b, root:= root, d:= d, ATLAS:= ATLAS );
  617.     end;
  618.  
  619.  
  620. #############################################################################
  621. ##
  622. #F  GeneratorsPrimeResidues( <N> ) . . . . . . generators of the Galois group
  623. ##
  624. ##  returns a record with fields 'primes' (the list of prime divisors of
  625. ##  'N'), 'exponents' (the list of exponents of these primes) and
  626. ##  'generators' (the list of generators of the prime--parts of the
  627. ##  group of prime residues; for p = 2 a list of two generators, otherwise
  628. ##  a primitive root).
  629. ##
  630. GeneratorsPrimeResidues := function( n )
  631.     local factors, exponents, primes, generators, ppart, rest, gcd, pos,
  632.           root, i;
  633.     
  634.     if n = 1 then return rec( primes:=[], exponent:=[], generators:=[] ); fi;
  635.     
  636.     factors:= Collected( FactorsInt( n ) );
  637.     
  638.     exponents := [];
  639.     primes    := [];
  640.     generators:= [];   # For each prime part 'ppart', the generator must be
  641.                        # congruent to a primitive root modulo 'ppart' and
  642.                        # congruent 1 modulo the rest 'N/ppart'.
  643.     
  644.     if factors[1][1] = 2 then
  645.       primes[1]   := 2;
  646.       exponents[1]:= factors[1][2];
  647.       ppart:= 2 ^ factors[1][2];
  648.       rest:= n / ppart;
  649.       gcd:= Gcdex( ppart, rest );
  650.       generators[1]:=   ( -2 * gcd.coeff2 * rest + 1 ) mod n; # generator '**'
  651.       if ppart mod 8 = 0 then 
  652.         generators[2]:= (  4 * gcd.coeff2 * rest + 1 ) mod n; # generator '*5'
  653.         generators:= [ generators ];
  654.       fi;
  655.       pos:= 2;
  656.     else
  657.       pos:= 1;
  658.     fi;
  659.     
  660.     for i in [ pos .. Length( factors ) ] do
  661.       ppart:= factors[i][1] ^ factors[i][2];
  662.       rest:= n / ppart;
  663.       gcd:= Gcdex( ppart, rest );
  664.       root:= PrimitiveRootMod( ppart );
  665.       primes[i]    := factors[i][1];
  666.       exponents[i] := factors[i][2];
  667.       generators[i]:= ( ( root - 1 ) * gcd.coeff2 * rest + 1 ) mod n;
  668.     od;
  669.     
  670.     return rec( primes    := primes,
  671.                 exponents := exponents,
  672.                 generators:= generators );
  673.     end;
  674.  
  675.  
  676. #############################################################################
  677. ##
  678. #F  GaloisMat( <mat> )
  679. ##
  680. ##  calculates the completions of orbits under the operation of the galois
  681. ##  group of the irrationalities of <mat>, and the permutations of rows
  682. ##  corresponding to the generators of the galois group.
  683. ##
  684. ##  If some rows of <mat> are identical, only the first one is considered
  685. ##  for the permutations, and a warning will be printed.
  686. ##
  687. GaloisMat := function( mat )
  688.     local i, j, k, l, m, generator, ncha, nccl, n, galoisfams, warned,
  689.         genexp, x, y, generators, X,
  690.         irrats, fusion, value, irratsimages, automs, family, orders, exp,
  691.         image, oldorder, cosets, auto, conj, blocklength, innerlength;
  692.     #
  693.     warned:= false;     # at most one warning will be printed if 'mat' grows
  694.     ncha:= Length( mat );
  695.     nccl:= Length( mat[1] );
  696.     mat:= ShallowCopy( mat );
  697.     
  698.     # step 1: find minimal cyclotomic field $Q_n$ containing all irrational
  699.     #         entries of <mat>:
  700.     
  701.     galoisfams:= [];     # will contain information about conjugate characters:
  702.                 # 1 means rational character,
  703.                 # -1 means character with undefs
  704.                 # 0 means dependent irrational character
  705.                 # [ .. ] means leading irrational character
  706.     n:= 1;
  707.     for i in [ 1 .. ncha ] do
  708.       if ForAny( mat[i], IsUnknown ) then
  709.         galoisfams[i]:= -1;
  710.       elif ForAll( mat[i], IsRat ) then
  711.         galoisfams[i]:= 1;
  712.       else
  713.         n:= LcmInt( n, NofCyc( mat[i] ) );
  714.       fi;
  715.     od;
  716.     
  717.     # step 2: compute generators for Aut( Q(n):Q )
  718.     #       ( compute generators for (Z/nZ)* and convert them to exponents )
  719.     
  720.     if n > 1 then
  721.     
  722.       genexp:= GeneratorsPrimeResidues( n ).generators;
  723.       if IsList( genexp[1] ) then
  724.         genexp:= Concatenation( genexp[1],
  725.                                 Sublist( genexp, [ 2..Length( genexp ) ] ) );
  726.       fi;
  727.     
  728.       generators:= []; # every galois automorphism induces a permutation of
  729.                        # rows, compute the permutations for each generating
  730.                        # automorphism
  731.       for i in [ 1 .. ncha ] do generators[i]:= i; od;
  732.       generators:= List( genexp, x -> Copy( generators ) );   # idperm
  733.     
  734.     else
  735.       generators:= [];    # rational matrix
  736.     fi;
  737.     
  738.     # step 3: for each character X find and complete the family of conjugates
  739.     
  740.     for i in [ 1 .. ncha ] do
  741.       if not IsBound( galoisfams[i] ) then      # independent character,
  742.                                                 # not integral, no undefs
  743.         X:= mat[i];
  744.         for j in [ i+1 .. ncha ] do
  745.           if mat[j] = X then
  746.             galoisfams[j]:= Unknown();
  747.     
  748.             Print( "#E GaloisMat: row ", i, " is equal to row ", j, "\n" );
  749.     
  750.           fi;
  751.         od;
  752.         irrats:= []; # list of distinct irrationalities of X (not ordered);
  753.              # every galois automorphism induces a permutation of that
  754.              # list rather than of X's entries themselves.
  755.         fusion:= []; # how to distribute the entries of irrats to X
  756.         for j in [ 1 .. nccl ] do
  757.           if IsCyc( X[j] ) and not IsRat( X[j] ) then
  758.             k:= 1;
  759.             while k <= Length( irrats ) and X[j] <> irrats[k] do k:= k+1; od;
  760.             if k > Length( irrats ) then      # first appearance of X[j] in X
  761.               irrats[k]:= X[j];
  762.             fi;
  763.             fusion[j]:= k;    # position in irrats
  764.           else
  765.             fusion[j]:= 0;
  766.           fi;
  767.         od;
  768.         irratsimages:= [ irrats ];
  769.         automs:= [ 1 ];
  770.         family:= [ i ];    # indices of family members (same ordering as automs)
  771.         orders:= [];    # orders[k] will be the order of the k-th generator
  772.         for j in [ 1 .. Length( genexp ) ] do
  773.           exp:= genexp[j];
  774.           image:= List( irrats, x -> GaloisCyc( x, exp mod NofCyc( x ) ) );
  775.           oldorder:= Length( automs );    # group order up to now
  776.           cosets:= [];
  777.           orders[j]:= 1;
  778.           while not image in irratsimages do
  779.             orders[j]:= orders[j] + 1;
  780.             for k in [ 1 .. oldorder ] do
  781.               auto:= ( automs[k] * exp ) mod n;
  782.               image:= List( irrats, x -> GaloisCyc( x, auto mod NofCyc(x) ) );
  783.               conj:= [];            # the conjugate character
  784.               for l in [ 1 .. nccl ] do
  785.             if fusion[l] = 0 then
  786.                   conj[l]:= X[l];
  787.             else
  788.                   conj[l]:= image[ fusion[l] ];
  789.             fi;
  790.               od;
  791.               l:= i+1;
  792.               while l <= ncha and mat[l] <> conj do l:= l+1; od;
  793.               if l <= ncha then
  794.               galoisfams[l]:= 0;
  795.                 Add( family, l );
  796.             for m in [ l+1 .. ncha ] do
  797.               if mat[m] = conj then galoisfams[m]:= 0; fi;
  798.             od;
  799.               else
  800.     
  801.                 if not warned and Length( mat ) > 500 then
  802.                   Print( "#I GaloisMat: completion of <mat> will have",
  803.                          " more than 500 rows\n" );
  804.                   warned:= true;
  805.                 fi;
  806.     
  807.              Add( mat, conj );
  808.             galoisfams[ Length( mat ) ]:= 0;
  809.                 Add( family, Length( mat ) );
  810.               fi;    
  811.               Add( automs, auto );
  812.               Add( cosets, image );
  813.             od;
  814.             exp:= exp * genexp[j];
  815.             image:= List( irrats, x -> GaloisCyc( x, exp mod NofCyc( x ) ) );
  816.           od;
  817.           irratsimages:= Concatenation( irratsimages, cosets );
  818.         od;
  819.         galoisfams[i]:= [ family, automs ];    # conjugates and automorphisms
  820.     
  821.         # Now the length of 'family' is the order of the Galois family of the
  822.         # row 'X'.
  823.     
  824.         # compute the permutation operation of the generators on the set of
  825.         # rows in 'family'\:
  826.     
  827.         blocklength:= 1;
  828.         for j in [ 1 .. Length( genexp ) ] do
  829.           innerlength:= blocklength;
  830.           blocklength:= blocklength * orders[j];
  831.           generator:= [ 1 .. blocklength ];
  832.     
  833.           # 'genexp[j]' maps the conjugates under the action of
  834.           # $\langle 'genexp[1]',\ldots,'genexp[j-1]' \rangle$
  835.           # (a set of length 'innerlength') as
  836.           # one block to their images and preserves the order of
  837.           # succession.
  838.     
  839.           for l in [ 1 .. blocklength - innerlength ] do
  840.             generator[l]:= l + innerlength;
  841.           od;
  842.     
  843.           # How does a power of 'genexp[j]' map back to the block:
  844.     
  845.           exp:= ( genexp[j] ^ orders[j] ) mod n;  # operates on the block
  846.           for l in [ 1 .. innerlength ] do
  847.             generator[ l + blocklength - innerlength ]:=
  848.                  Position( irratsimages, List( irrats,
  849.                            x->GaloisCyc( x, exp*automs[l] mod NofCyc(x) ) ) );
  850.           od;
  851.     
  852.           # Transfer this operation to the cosets under the operation of
  853.           # $\langle 'genexp[j+1],\ldots,'genexp[ Length( genexp ) ]' \rangle$
  854.           # and transfer this to <mat> via 'family'\:
  855.     
  856.           for k in [ 0 .. Length( family ) / blocklength - 1 ] do
  857.             for l in [ 1 .. blocklength ] do
  858.               generators[j][ family[ l + k*blocklength ] ]:=
  859.                            family[ generator[ l ] + k*blocklength ];
  860.             od;
  861.           od;
  862.         od;
  863.     
  864.       fi;
  865.     od;
  866.     generators:= Set( List( generators, x-> PermList(x) ) );
  867.     RemoveSet( generators, () );
  868.     if generators = [] then generators:= [ () ]; fi;
  869.     return rec( mat:= mat,
  870.             galoisfams:= galoisfams,
  871.             generators:= generators );
  872.     end;
  873.       
  874.  
  875. #############################################################################
  876. ##
  877. #F  RationalizedMat( <mat> ) . .  list of rationalized rows of <mat>
  878. ##
  879. RationalizedMat := function( mat )
  880.     local i, x, rationalizedmat, rationalfams, mat;
  881.     rationalizedmat:= [];
  882.     rationalfams:= GaloisMat( mat );
  883.     mat:= rationalfams.mat;
  884.     rationalfams:= rationalfams.galoisfams;
  885.     for i in [ 1 .. Length( mat ) ] do
  886.       if rationalfams[i] in [ 1, -1 ] then     # rational or with undefs
  887.         Add( rationalizedmat, Copy( mat[i] ) );
  888.       elif IsList( rationalfams[i] ) then      # leading character
  889.         Add( rationalizedmat, Sum( Sublist( mat, rationalfams[i][1] ) ) );
  890.       fi;
  891.     od;
  892.     return rationalizedmat;
  893.     end;
  894.  
  895.  
  896.