// Functions used in paper 1:
// _Some computational experiments in number theory_, by Wieb Bosma

// Section 2: Covering systems
forward addm, addr;

try2cover := function(S, D)
  Z := Integers();
  for tries := 1 to 50 do
    T := [ ];
    Q := [ 1 : i in [1..D] ];
    for i in [1..#S] do
      addm(~Q, ~T, S[i]);
      if &+[ Z | D div s : s in S[i+1..#S] ] lt &+Q then
        if &+Q/D gt 0.1 then
          break tries;
        end if;
        break;
      elif &+Q eq 0 then
        return true, T;
      end if;
    end for;
  end for;
  return false, _;
end function;

addm := procedure(~Q, ~T, m)
  if &+Q eq #Q then
    addr(~Q, ~T, 1, m);
  else
    mx := [ &+Q[[i..#Q by m]] : i in [1..m] ];
    mm := Max(mx);
    im := Random([ i : i in [1..#mx] | mx[i] eq mm ]);
    addr(~Q, ~T, im-1, m);
  end if;
end procedure;

// This next function (addr) was not given in the text
addr := procedure(~Q, ~T, r, m)
  N := #Q;
  Append(~T, [r, m]);
  if N eq 0 then
    Q := [ 1 : i in [0..m-1] ];
    Q[r+1] := 0;
  else
    if N mod m ne 0 then
      N := Lcm(N, m);
      Q := &cat[ Q : i in [1..(N div #Q)] ];
    end if;
    for i := 0 to (N div m)-1 do
      Q[ r + i*m + 1 ] := 0;
    end for;
  end if;
end procedure;


// Section 3: Covering systems and primality proving
forward cut, get, cubicsymbol;

minfind := function(h, bound, PX)
  i := 0;
  K := [1];
  repeat
    i +:= 1;
    N := h*3^i - 1;
    _, I := get(N, PX);
    if IsEmpty(I) then return 0; end if;
    J := cut(I);
    K := Sort([ Lcm(k,j) : k in K, j in J | Lcm(k,j) lt bound ]);
    if IsEmpty(K) then
      return 0;
    else
      K := cut(K);
    end if;
  until i in K;
  return i;
end function;

// This next function (plusfind) was not given in the text, although
// it was stated to be a trivial modification of minfind.
plusfind := function(h, bound, PX)
  i := 0;
  K := [1];
  repeat
    i +:= 1;
    N := h*3^i + 1;
    _, I := get(N, PX);
    if IsEmpty(I) then return 0; end if;
    J := cut(I);
    K := Sort([ Lcm(k,j) : k in K, j in J | Lcm(k,j) lt bound ]);
    if IsEmpty(K) then
      return 0;
    else
      K := cut(K);
    end if;
  until i in K;
  return i;
end function;

Qz<zeta> := CyclotomicField(3);

get := function(N, PX)
  S := [ Parent(zeta) | ];
  OS := [ ];
  for i in [1..#PX] do
    if IsEmpty(PX[i]) then
      continue;
    end if;
    for x in PX[i] do
      if cubicsymbol(x[1], x[2], N, 0) ne 1 then
        if x notin S then
          Append(~S, x);
          Append(~OS, i);
        end if;
      end if;
    end for;
  end for;
  return S, OS;
end function;

// This next function (cut) was not given in the text
cut := function(I)
  J := [ I[1] ];
  for i := 2 to #I do
    for j in J do
      if I[i] mod j eq 0 then
        continue i;
      end if;
    end for;
    Append(~J, I[i]);
  end for;
  return J;
end function;

// This next function (cubicsymbol) was not given in the text
// It uses the functions gcd, primary, quotrem, div3 as part of its
// operation over Z[zeta] and an explicit representation of elements
// as a + b*zeta.
forward gcd, primary, quotrem, div3;

cubicsymbol := function(a, b, x, y, zeta)
/* determine the cubic residue symbol as the integer 0, or 1, or
   the first or second power of zeta_3
   assume that x + y*zeta_3 is not 0 and not divisible by 1-zeta_3
*/
  res := 1;
  nd := x^2 - x*y + y^2;        // norm denominator
  if a eq 0 and b eq 0 then
    return 0;
  elif nd eq 1 then
    return 1;
  end if;
  g1, g2 := gcd(a,b, x,y);
  if not (g2 eq 0 and g1 eq -1) then
    return 0;
  end if;
  u := a;  v := b;
  s,t := primary(x,y);
  while nd ge 4 do
    temp1 := s; temp2 := t;
    _, _, s, t := quotrem(u,v, s,t);
    u := temp1; v := temp2;
    k, s,t := div3(s,t);
    if k gt 0 then
      res *:= zeta^(k*2*((u+1) div 3));
    end if;
    s,t, sg, ex := primary(s,t);
    if ex gt 0 then
      res *:= zeta^(ex*(((u+1) div 3) + (v div 3)));
    end if;
    nd := s^2 - s*t + t^2;
  end while;
  return res;
end function;

gcd := function(a,b,r,s)
 nm := r^2 - r*s + s^2;
 while nm gt 0 do
  u := r; v := s;
  _,_, r,s := quotrem(a,b, r,s);
  a := u; b := v;
  nm := r^2 - r*s + s^2;
 end while;
 a,b := primary(a,b);
 return a,b;
end function;

primary := function(x, y)
/* given element x + y*zeta_3, with x,y in Z
   return integers a,b and s,e, such that
   (x + y*zeta_3) = s*zeta_3^e*(a + b*zeta_3) and
   s in {1,-1}, e in {0,1,2} and a  + b*zeta_3 primary,
   that is, it is 2 mod 3
*/
  if (x mod 3) eq 0 then
    if (y mod 3) eq 1 then return x-y, x, -1, 1;
    elif (y mod 3) eq 2 then return -x+y, -x, 1, 1;
    end if;
  elif (y mod 3) eq 0 then
    if (x mod 3) eq 1 then return -x, -y, -1, 0;
    elif (x mod 3) eq 2 then return x, y, 1, 0;
    end if;
  else
    if (y mod 3) eq 1 then return -y, x-y, 1, 2;
    elif (y mod 3) eq 2 then return y, y-x, -1, 2;
    end if;
  end if;
end function;

quotrem := function(x,y, z,w)
  n := z^2 + w^2 - z*w;
  p1 := x*(z-w) + y*w;
  p2 := -x*w + y*(z-w) + y*w;
  q1 := Round(p1/n);
  q2 := Round(p2/n);
  r1 := x - (q1*z - q2*w);
  r2 := y - (q2*z + q1*w - q2*w);
  return q1, q2, r1, r2;
end function;

div3 := function(a, b)
/* for a + b*zeta_3 this returns non-negative k and integers x, y
   such that (a + b*zeta_3) = (1 - zeta_3)^k*(x + y*zeta_3) and
   x + y*zeta_3) not divisible by (1 - zeta_3) [ the prime over 3 ]
*/
  x := a; y := b;
  k := 0;
  while x mod 3 eq 0 and y mod 3 eq 0 do
    x := (x - y); y := (x + y);
    x := x div 3; y := y div 3;
    k +:= 2;
  end while;
  while (x^2 - x*y + y^2) mod 3 eq 0 do
    y := (x + y); x := (-y + 3*x);
    x := x div 3; y := y div 3;
    k +:= 1;
  end while;
  return k, x,y;
end function;

// This next function (check) was not given in the text
check := function(O)
  n := Lcm([ o[2] : o in O ]);
  N := [ 1 : i in [1..n] ];
  for o in O do
    for r := o[1] to n by o[2] do N[r+1] := 0; end for;
  end for;
  return IsEmpty([ n-1 : n in [1..#N] | N[n] eq 1 ]);
end function;


// Section 4: The totient function
inv := function(m)
  mfact := Factorization(m);
  if IsEven(m) then
    twopows := {@ 2^i : i in [0..mfact[1][2]] @};
  else
    if m gt 1 then return [ ]; end if;
    twopows := {@ 1 @};
  end if;
  D := Divisors(mfact);
  P := [ ];
  for d in D do
   if d eq 1 then continue; end if;
   if IsPrime(d+1) then
     Append(~P, d+1);
   end if;
  end for;
  S := [ <SeqFact([ ]), m> ];
  for p in Reverse(P) do
    for s in S do
      if s[2] eq 1 then continue; end if;
      k := 1;
      d, mmod := Quotrem(s[2], p-1);
      while mmod eq 0 do
        if IsEven(d) or d eq 1 then
          Append(~S, <SeqFact([<p, k>])*s[1], d>);
        end if;
        k +:= 1;
        d, mmod := Quotrem(d, p);
      end while;
    end for;
  end for;
  R := { };
  for s in S do
    j := Index(twopows, s[2]);
    if j gt 0 then
      Include(~R, SeqFact([<2, j>] cat s[1]));
      if j eq 1 then
        Include(~R, s[1]);
      end if;
    end if;
  end for;
  return Sort([ Facint(nf) : nf in R ]);
end function;


// Section 5: Class number relations
permcharmat := function(G)
  nc := #ConjugacyClasses(G);
  subs := Subgroups(G);
  M := Hom(RSpace(Integers(), #subs), RSpace(Integers(), nc));
  return M !
    &cat[ Eltseq(PermutationCharacter(G, s`subgroup)) : s in subs ];
end function;