15 Number Theory
  
  GAP  provides  a  couple  of  elementary number theoretic functions. Most of
  these deal with the group of integers coprime to m, called the prime residue
  group.  The  order  of  this  group  is  ϕ(m)  (see Phi  (15.2-2)), and λ(m)
  (see Lambda (15.2-3)) is its exponent. This group is cyclic if and only if m
  is  2, 4, an odd prime power p^n, or twice an odd prime power 2 p^n. In this
  case  the  generators of the group, i.e., elements of order ϕ(m), are called
  primitive roots (see PrimitiveRootMod (15.3-3)).
  
  Note  that  neither  the  arguments  nor  the return values of the functions
  listed below are groups or group elements in the sense of GAP. The arguments
  are simply integers.
  
  
  15.1 InfoNumtheor (Info Class)
  
  15.1-1 InfoNumtheor
  
  InfoNumtheor info class
  
  InfoNumtheor  is  the  info  class (see 7.4) for the functions in the number
  theory chapter.
  
  
  15.2 Prime Residues
  
  15.2-1 PrimeResidues
  
  PrimeResidues( m )  function
  
  PrimeResidues returns the set of integers from the range [ 0 .. Abs( m )-1 ]
  that are coprime to the integer m.
  
  Abs(m) must be less than 2^28, otherwise the set would probably be too large
  anyhow.
  
    Example  
    gap> PrimeResidues( 0 );  PrimeResidues( 1 );  PrimeResidues( 20 );
    [  ]
    [ 0 ]
    [ 1, 3, 7, 9, 11, 13, 17, 19 ]
  
  
  15.2-2 Phi
  
  Phi( m )  operation
  
  Phi  returns  the  number  ϕ(m)  of positive integers less than the positive
  integer m that are coprime to m.
  
  Suppose  that  m = p_1^{e_1} p_2^{e_2} ⋯ p_k^{e_k}. Then ϕ(m) is p_1^{e_1-1}
  (p_1-1) p_2^{e_2-1} (p_2-1) ⋯ p_k^{e_k-1} (p_k-1).
  
    Example  
    gap> Phi( 12 );
    4
    gap> Phi( 2^13-1 );  # this proves that 2^(13)-1 is a prime
    8190
    gap> Phi( 2^15-1 );
    27000
  
  
  15.2-3 Lambda
  
  Lambda( m )  operation
  
  Lambda  returns  the exponent λ(m) of the group of prime residues modulo the
  integer m.
  
  λ(m)  is  the  smallest  positive integer l such that for every a relatively
  prime  to m we have a^l ≡ 1 mod m. Fermat's theorem asserts a^{ϕ(m)} ≡ 1 mod
  m; thus λ(m) divides ϕ(m) (see Phi (15.2-2)).
  
  Carmichael's  theorem  states  that  λ can be computed as follows: λ(2) = 1,
  λ(4)  =  2 and λ(2^e) = 2^{e-2} if 3 ≤ e, λ(p^e) = (p-1) p^{e-1} (i.e. ϕ(m))
  if p is an odd prime and λ(m*n) =Lcm( λ(m), λ(n) ) if m, n are coprime.
  
  Composites  for  which  λ(m)  divides m - 1 are called Carmichaels. If 6k+1,
  12k+1  and  18k+1  are primes their product is such a number. There are only
  1547 Carmichaels below 10^10 but 455052511 primes.
  
    Example  
    gap> Lambda( 10 );
    4
    gap> Lambda( 30 );
    4
    gap> Lambda( 561 );  # 561 is the smallest Carmichael number
    80
  
  
  15.2-4 GeneratorsPrimeResidues
  
  GeneratorsPrimeResidues( n )  function
  
  Let  n  be a positive integer. GeneratorsPrimeResidues returns a description
  of generators of the group of prime residues modulo n. The return value is a
  record with components
  
  primes: 
        a list of the prime factors of n,
  
  exponents: 
        a list of the exponents of these primes in the factorization of n, and
  
  generators: 
        a  list  describing generators of the group of prime residues; for the
        prime factor 2, either a primitive root or a list of two generators is
        stored, for each other prime factor of n, a primitive root is stored.
  
    Example  
    gap> GeneratorsPrimeResidues( 1 );
    rec( exponents := [  ], generators := [  ], primes := [  ] )
    gap> GeneratorsPrimeResidues( 4*3 );
    rec( exponents := [ 2, 1 ], generators := [ 7, 5 ], 
      primes := [ 2, 3 ] )
    gap> GeneratorsPrimeResidues( 8*9*5 );
    rec( exponents := [ 3, 2, 1 ], 
      generators := [ [ 271, 181 ], 281, 217 ], primes := [ 2, 3, 5 ] )
  
  
  
  15.3 Primitive Roots and Discrete Logarithms
  
  15.3-1 OrderMod
  
  OrderMod( n, m )  function
  
  OrderMod  returns  the  multiplicative  order  of  the  integer n modulo the
  positive integer m. If n and m are not coprime the order of n is not defined
  and OrderMod will return 0.
  
  If  n  and  m are relatively prime the multiplicative order of n modulo m is
  the  smallest  positive  integer  i such that n^i ≡ 1 mod m. If the group of
  prime  residues  modulo  m  is  cyclic then each element of maximal order is
  called a primitive root modulo m (see IsPrimitiveRootMod (15.3-4)).
  
  OrderMod   usually   spends   most   of   its  time  factoring  m  and  ϕ(m)
  (see FactorsInt (14.4-7)).
  
    Example  
    gap> OrderMod( 2, 7 );
    3
    gap> OrderMod( 3, 7 );  # 3 is a primitive root modulo 7
    6
  
  
  15.3-2 LogMod
  
  LogMod( n, r, m )  function
  LogModShanks( n, r, m )  function
  
  computes  the discrete r-logarithm of the integer n modulo the integer m. It
  returns  a  number  l  such  that  r^l  ≡  n  mod m if such a number exists.
  Otherwise fail is returned.
  
  LogModShanks  uses  the  Baby  Step  -  Giant Step Method of Shanks (see for
  example  [Coh93,  section 5.4.1]) and in general requires more memory than a
  call to LogMod.
  
    Example  
    gap> l:= LogMod( 2, 5, 7 );  5^l mod 7 = 2;
    4
    true
    gap> LogMod( 1, 3, 3 );  LogMod( 2, 3, 3 );
    0
    fail
  
  
  15.3-3 PrimitiveRootMod
  
  PrimitiveRootMod( m[, start] )  function
  
  PrimitiveRootMod  returns  the  smallest  primitive root modulo the positive
  integer  m and fail if no such primitive root exists. If the optional second
  integer  argument  start  is  given  PrimitiveRootMod  returns  the smallest
  primitive root that is strictly larger than start.
  
    Example  
    gap> # largest primitive root for a prime less than 2000:
    gap> PrimitiveRootMod( 409 ); 
    21
    gap> PrimitiveRootMod( 541, 2 );
    10
    gap> # 327 is the largest primitive root mod 337:
    gap> PrimitiveRootMod( 337, 327 );
    fail
    gap> # there exists no primitive root modulo 30:
    gap> PrimitiveRootMod( 30 );
    fail
  
  
  15.3-4 IsPrimitiveRootMod
  
  IsPrimitiveRootMod( r, m )  function
  
  IsPrimitiveRootMod  returns true if the integer r is a primitive root modulo
  the  positive  integer m, and false otherwise. If r is less than 0 or larger
  than m it is replaced by its remainder.
  
    Example  
    gap> IsPrimitiveRootMod( 2, 541 );
    true
    gap> IsPrimitiveRootMod( -539, 541 );  # same computation as above;
    true
    gap> IsPrimitiveRootMod( 4, 541 );
    false
    gap> ForAny( [1..29], r -> IsPrimitiveRootMod( r, 30 ) );
    false
    gap> # there is no a primitive root modulo 30
  
  
  
  15.4 Roots Modulo Integers
  
  15.4-1 Jacobi
  
  Jacobi( n, m )  function
  
  Jacobi  returns  the  value  of  the  Kronecker-Jacobi  symbol J(n,m) of the
  integer n modulo the integer m. It is defined as follows:
  
  If  n  and  m  are  not coprime then J(n,m) = 0. Furthermore, J(n,1) = 1 and
  J(n,-1)  = -1 if m < 0 and +1 otherwise. And for odd n it is J(n,2) = (-1)^k
  with  k  =  (n^2-1)/8.  For  odd  primes  m  which  are  coprime  to  n  the
  Kronecker-Jacobi   symbol   has  the  same  value  as  the  Legendre  symbol
  (see Legendre (15.4-2)).
  
  For the general case suppose that m = p_1 ⋅ p_2 ⋯ p_k is a product of -1 and
  of primes, not necessarily distinct, and that n is coprime to m. Then J(n,m)
  = J(n,p_1) ⋅ J(n,p_2) ⋯ J(n,p_k).
  
  Note  that the Kronecker-Jacobi symbol coincides with the Jacobi symbol that
  is  defined  for  odd  m in many number theory books. For odd primes m and n
  coprime to m it coincides with the Legendre symbol.
  
  Jacobi  is  very efficient, even for large values of n and m, it is about as
  fast as the Euclidean algorithm (see Gcd (56.7-1)).
  
    Example  
    gap> Jacobi( 11, 35 );  # 9^2 = 11 mod 35
    1
    gap> # this is -1, thus there is no r such that r^2 = 6 mod 35
    gap> Jacobi( 6, 35 );
    -1
    gap> # this is 1 even though there is no r with r^2 = 3 mod 35
    gap> Jacobi( 3, 35 );
    1
  
  
  15.4-2 Legendre
  
  Legendre( n, m )  function
  
  Legendre  returns  the  value of the Legendre symbol of the integer n modulo
  the positive integer m.
  
  The  value  of  the  Legendre symbol L(n/m) is 1 if n is a quadratic residue
  modulo  m, i.e., if there exists an integer r such that r^2 ≡ n mod m and -1
  otherwise.
  
  If a root of n exists it can be found by RootMod (15.4-3).
  
  While  the  value  of  the  Legendre  symbol usually is only defined for m a
  prime,  we have extended the definition to include composite moduli too. The
  Jacobi  symbol  (see  Jacobi  (15.4-1))  is  another  generalization  of the
  Legendre  symbol  for  composite  moduli  that  is  much cheaper to compute,
  because it does not need the factorization of m (see FactorsInt (14.4-7)).
  
  A  description of the Jacobi symbol, the Legendre symbol, and related topics
  can be found in [Bak84].
  
    Example  
    gap> Legendre( 5, 11 );  # 4^2 = 5 mod 11
    1
    gap> # this is -1, thus there is no r such that r^2 = 6 mod 11
    gap> Legendre( 6, 11 );
    -1
    gap> # this is -1, thus there is no r such that r^2 = 3 mod 35
    gap> Legendre( 3, 35 );
    -1
  
  
  15.4-3 RootMod
  
  RootMod( n[, k], m )  function
  
  RootMod  computes a kth root of the integer n modulo the positive integer m,
  i.e.,  a  r  such that r^k ≡ n mod m. If no such root exists RootMod returns
  fail. If only the arguments n and m are given, the default value for k is 2.
  
  A square root of n exists only if Legendre(n,m) = 1 (see Legendre (15.4-2)).
  If  m  has r different prime factors then there are 2^r different roots of n
  mod  m.  It  is unspecified which one RootMod returns. You can, however, use
  RootsMod (15.4-4) to compute the full set of roots.
  
  RootMod  is  efficient  even for large values of m, in fact the most time is
  usually spent factoring m (see FactorsInt (14.4-7)).
  
    Example  
    gap> # note 'RootMod' does not return 8 in this case but -8:
    gap> RootMod( 64, 1009 );
    1001
    gap> RootMod( 64, 3, 1009 );
    518
    gap> RootMod( 64, 5, 1009 );
    656
    gap> List( RootMod( 64, 1009 ) * RootsUnityMod( 1009 ),
    >       x -> x mod 1009 );  # set of all square roots of 64 mod 1009
    [ 1001, 8 ]
  
  
  15.4-4 RootsMod
  
  RootsMod( n[, k], m )  function
  
  RootsMod  computes the set of kth roots of the integer n modulo the positive
  integer  m,  i.e.,  the  list  of all r such that r^k ≡ n mod m. If only the
  arguments n and m are given, the default value for k is 2.
  
    Example  
    gap> RootsMod( 1, 7*31 );  # the same as `RootsUnityMod( 7*31 )'
    [ 1, 92, 125, 216 ]
    gap> RootsMod( 7, 7*31 );
    [ 21, 196 ]
    gap> RootsMod( 5, 7*31 );
    [  ]
    gap> RootsMod( 1, 5, 7*31 );
    [ 1, 8, 64, 78, 190 ]
  
  
  15.4-5 RootsUnityMod
  
  RootsUnityMod( [k, ]m )  function
  
  RootsUnityMod  returns  the  set  of k-th roots of unity modulo the positive
  integer  m,  i.e., the list of all solutions r of r^k ≡ n mod m. If only the
  argument m is given, the default value for k is 2.
  
  In  general  there are k^n such roots if the modulus m has n different prime
  factors  p  such  that  p ≡ 1 mod k. If k^2 divides m then there are k^{n+1}
  such  roots;  and especially if k = 2 and 8 divides m there are 2^{n+2} such
  roots.
  
  In the current implementation k must be a prime.
  
    Example  
    gap> RootsUnityMod( 7*31 );  RootsUnityMod( 3, 7*31 );
    [ 1, 92, 125, 216 ]
    [ 1, 25, 32, 36, 67, 149, 156, 191, 211 ]
    gap> RootsUnityMod( 5, 7*31 );
    [ 1, 8, 64, 78, 190 ]
    gap> List( RootMod( 64, 1009 ) * RootsUnityMod( 1009 ),
    >          x -> x mod 1009 );  # set of all square roots of 64 mod 1009
    [ 1001, 8 ]
  
  
  
  15.5 Multiplicative Arithmetic Functions
  
  15.5-1 Sigma
  
  Sigma( n )  operation
  
  Sigma returns the sum of the positive divisors of the nonzero integer n.
  
  Sigma  is  a  multiplicative  arithmetic  function,  i.e.,  if  n  and m are
  relatively prime we have that σ(n ⋅ m) = σ(n) σ(m).
  
  Together  with  the  formula  σ(p^k) = (p^{k+1}-1) / (p-1) this allows us to
  compute σ(n).
  
  Integers  n  for  which σ(n) = 2 n are called perfect. Even perfect integers
  are  exactly  of the form 2^{n-1}(2^n-1) where 2^n-1 is prime. Primes of the
  form  2^n-1  are  called  Mersenne  primes,  and 42 among the known Mersenne
  primes  are  obtained  for n = 2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127,
  521,  607,  1279,  2203,  2281,  3217, 4253, 4423, 9689, 9941, 11213, 19937,
  21701, 23209, 44497, 86243, 110503, 132049, 216091, 756839, 859433, 1257787,
  1398269,   2976221,  3021377,  6972593,  13466917,  20996011,  24036583  and
  25964951.  Please  find more up to date information about Mersenne primes at
  http://www.mersenne.org. It is not known whether odd perfect integers exist,
  however [BC89]  show  that  any  such integer must have at least 300 decimal
  digits.
  
  Sigma usually spends most of its time factoring n (see FactorsInt (14.4-7)).
  
    Example  
    gap> Sigma( 1 );
    1
    gap> Sigma( 1009 );  # 1009 is a prime
    1010
    gap> Sigma( 8128 ) = 2*8128;  # 8128 is a perfect number
    true
  
  
  15.5-2 Tau
  
  Tau( n )  operation
  
  Tau returns the number of the positive divisors of the nonzero integer n.
  
  Tau  is  a multiplicative arithmetic function, i.e., if n and m are relative
  prime  we  have τ(n ⋅ m) = τ(n) τ(m). Together with the formula τ(p^k) = k+1
  this allows us to compute τ(n).
  
  Tau usually spends most of its time factoring n (see FactorsInt (14.4-7)).
  
    Example  
    gap> Tau( 1 );
    1
    gap> Tau( 1013 );  # thus 1013 is a prime
    2
    gap> Tau( 8128 );
    14
    gap> # result is odd if and only if argument is a perfect square:
    gap> Tau( 36 );
    9
  
  
  15.5-3 MoebiusMu
  
  MoebiusMu( n )  function
  
  MoebiusMu  computes  the value of Moebius inversion function for the nonzero
  integer  n. This is 0 for integers which are not squarefree, i.e., which are
  divided  by a square r^2. Otherwise it is 1 if n has a even number and -1 if
  n has an odd number of prime factors.
  
  The importance of μ stems from the so called inversion formula. Suppose f is
  a  multiplicative  arithmetic  function defined on the positive integers and
  let  g(n)  = ∑_{d ∣ n} f(d). Then f(n) = ∑_{d ∣ n} μ(d) g(n/d). As a special
  case  we  have  ϕ(n)  = ∑_{d ∣ n} μ(d) n/d since n = ∑_{d ∣ n} ϕ(d) (see Phi
  (15.2-2)).
  
  MoebiusMu  usually  spends  all  of  its  time  factoring  n (see FactorsInt
  (14.4-7)).
  
    Example  
    gap> MoebiusMu( 60 );  MoebiusMu( 61 );  MoebiusMu( 62 );
    0
    -1
    1
  
  
  
  15.6 Continued Fractions
  
  15.6-1 ContinuedFractionExpansionOfRoot
  
  ContinuedFractionExpansionOfRoot( f, n )  function
  
  The  first  n terms of the continued fraction expansion of the only positive
  real  root  of  the  polynomial  f  with  integer  coefficients. The leading
  coefficient  of f must be positive and the value of f at 0 must be negative.
  If  the  degree of f is 2 and n = 0, the function computes one period of the
  continued fraction expansion of the root in question. Anything may happen if
  f has three or more positive real roots.
  
    Example  
    gap> x := Indeterminate(Integers);;
    gap> ContinuedFractionExpansionOfRoot(x^2-7,20);
    [ 2, 1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1 ]
    gap> ContinuedFractionExpansionOfRoot(x^2-7,0);
    [ 2, 1, 1, 1, 4 ]
    gap> ContinuedFractionExpansionOfRoot(x^3-2,20);
    [ 1, 3, 1, 5, 1, 1, 4, 1, 1, 8, 1, 14, 1, 10, 2, 1, 4, 12, 2, 3 ]
    gap> ContinuedFractionExpansionOfRoot(x^5-x-1,50);
    [ 1, 5, 1, 42, 1, 3, 24, 2, 2, 1, 16, 1, 11, 1, 1, 2, 31, 1, 12, 5, 
      1, 7, 11, 1, 4, 1, 4, 2, 2, 3, 4, 2, 1, 1, 11, 1, 41, 12, 1, 8, 1, 
      1, 1, 1, 1, 9, 2, 1, 5, 4 ]
  
  
  15.6-2 ContinuedFractionApproximationOfRoot
  
  ContinuedFractionApproximationOfRoot( f, n )  function
  
  The  nth  continued fraction approximation of the only positive real root of
  the  polynomial  f  with  integer coefficients. The leading coefficient of f
  must  be  positive  and  the  value of f at 0 must be negative. Anything may
  happen if f has three or more positive real roots.
  
    Example  
    gap> ContinuedFractionApproximationOfRoot(x^2-2,10);
    3363/2378
    gap> 3363^2-2*2378^2;
    1
    gap> z := ContinuedFractionApproximationOfRoot(x^5-x-1,20);
    499898783527/428250732317
    gap> z^5-z-1;
    486192462527432755459620441970617283/
    14404247382319842421697357558805709031116987826242631261357
  
  
  
  15.7 Miscellaneous
  
  15.7-1 TwoSquares
  
  TwoSquares( n )  function
  
  TwoSquares  returns  a  list  of two integers x ≤ y such that the sum of the
  squares  of  x  and y is equal to the nonnegative integer n, i.e., n = x^2 +
  y^2.   If  no  such  representation  exists  TwoSquares  will  return  fail.
  TwoSquares  will  return a representation for which the gcd of x and y is as
  small  as  possible.  It  is  not  specified which representation TwoSquares
  returns if there is more than one.
  
  Let  a  be  the  product  of  all  maximal powers of primes of the form 4k+3
  dividing n. A representation of n as a sum of two squares exists if and only
  if  a is a perfect square. Let b be the maximal power of 2 dividing n or its
  half,  whichever is a perfect square. Then the minimal possible gcd of x and
  y  is  the  square  root  c  of  a  ⋅  b.  The  number  of different minimal
  representation  with  x  ≤  y is 2^{l-1}, where l is the number of different
  prime factors of the form 4k+1 of n.
  
  The  algorithm  first  finds a square root r of -1 modulo n / (a ⋅ b), which
  must  exist,  and  applies  the  Euclidean  algorithm  to r and n. The first
  residues in the sequence that are smaller than sqrt{n/(a ⋅ b)} times c are a
  possible pair x and y.
  
  Better  descriptions  of  the  algorithm  and related topics can be found in
  [Wag90] and [Zag90].
  
    Example  
    gap> TwoSquares( 5 );
    [ 1, 2 ]
    gap> TwoSquares( 11 );  # there is no representation
    fail
    gap> TwoSquares( 16 );
    [ 0, 4 ]
    gap> # 3 is the minimal possible gcd because 9 divides 45:
    gap> TwoSquares( 45 );
    [ 3, 6 ]
    gap> # it is not [5,10] because their gcd is not minimal:
    gap> TwoSquares( 125 );
    [ 2, 11 ]
    gap> # [10,11] would be the other possible representation:
    gap> TwoSquares( 13*17 );
    [ 5, 14 ]
    gap> TwoSquares( 848654483879497562821 );  # argument is prime
    [ 6305894639, 28440994650 ]