// Example code from paper 1:
// _Some computational experiments in number theory_ by Wieb Bosma
//
print "Examples from: Some computational experiments in number theory";
load "Paper01_funcs";
ei := GetEchoInput();
SetEchoInput(true);
// Section 2: Covering systems
for c := 2 to 10 do
  D := 0;
  done := false;
  repeat
    D +:= 4;
    S := [ x : x in Divisors(D) | x ge c ];
    if  &+[ Integers() | D div s : s in S ] ge D then
      done, F := try2cover(S, D);
    end if;
  until done;
  print < Min([ f[2] : f in F ]), D, F >;
end for;


m := [ 0, 0, 1, 3, 7, 23 ];
P := [ 3, 7, 5, 17, 13, 241 ];
N := CRT([ 2^j : j in m ], P); N;
N + &*P;
Q := &*P; Q;
H := InverseMod(-N, Q); H;
for k := 0 to 100 do
  print k, [ p : p in P | (H*2^k+1) mod p eq 0 ];
end for;


// Section 3: Covering systems and explicit primality tests
S := [ 2,3,4,9,12,18,36 ];  CS := [ ];
K := CartesianProduct([ Integers(i) : i in S ]);
for x in K do
  C := [ [ Integers() ! x[i], S[i] ]: i in [1..#S] ];
  if check(C) then
    Append(~CS, C);
  end if;
end for;
P := [ 3, 7, 5, 73, 13, 19, 109 ];
H := CRT([ -Modexp(2, -C[i][1], P[i]) : i in [1..#C] ], P);
P := [ 3, 7, 5, 73, 13, 19, 37 ];
H := CRT([ -Modexp(2, -C[i][1], P[i]) : i in [1..#C] ], P);
C := [ [1,3], [2,4], [2,6], [3,8], [0,12], [7,16], [5,16], [18,24],
       [30,48], [6,48] ];
P := [ 13, 5, 7, 41, 73, 17, 193, 6481, 97, 577 ];
H := CRT([ -Modexp(3, -C[i][1], P[i]) : i in [1..#C] ], P);


// Section 4: The totient function
A := 35;
P := [ n : n in [2..A] | IsPrime(n) ];
E := [ Floor(A/p) : p in P ];
C := CartesianProduct([ [e..0 by -1] : e in E ]);
for c in C do
  nfn := SeqFact([ <P[i], c[i]> : i in [1..#P] | c[i] ne 0 ]);
  if EulerPhi(DivisorSigma(1, nfn)) eq Facint(nfn) then
    print nfn;
  end if;
end for;


// Section 5: Class number relations
permcharmat( Alt(4) );
relations := func< G | Kernel( permcharmat(G) ) >;
relations( Alt(4) );


G := Sym(6);
U := sub< G | (1,2)(3,4), (1,3)(2,4), (1,4)(2,3) >;
V := sub< G | (1,2)(3,4), (1,2)(5,6), (3,4)(5,6) >;
PermutationCharacter(G, U);
PermutationCharacter(G, V);
Induction(PrincipalCharacter(U), G) eq PermutationCharacter(G, U);


for n := 1 to 12 do
  for k := 1 to NumberOfTransitiveGroups(n) do
    G := TransitiveGroup(n, k);
    U := Stabilizer(G, 1);
    chi := PermutationCharacter(G, U);
    S := Subgroups(G : OrderEqual := Order(G) div n);
    if exists(i){ i :  i in [1..#S] |
                  PermutationCharacter(G, V) eq chi and IsEmpty(Fix(V))
                  where V := S[i]`subgroup }
    then
      < n, k, Order(G) >;
    end if;
  end for;
end for;
// The next line is a modified print line for the above, as described
// at the bottom of page 24:
// < n, k, G, U, S[i]`subgroup >;


F<s,t> := FunctionField(Rationals(), 2);
Q<x> := PolynomialRing(F);
f := x^7 + (-6*t + 2)*x^6 + (8*t^2 + 4*t - 3)*x^5 +
     (-s - 14*t^2 + 6*t - 2)*x^4 + (s + 6*t^2 - 8*t^3 - 4*t + 2)*x^3 +
     (8*t^3 + 16*t^2)*x^2 + (8*t^3 - 12*t^2)*x - 8*t^3;

h1 := (-1)^Degree(f)*Evaluate(f, -x);
Y<y> := PolynomialRing(Q);
r := Resultant(Evaluate(h1, y-x), Evaluate(f, y));
h2 := 2^7*Evaluate(f, x/2);
q1 := SquareRoot(r div h2);

q := Resultant(Evaluate(h1, y-x), Evaluate(h2, y));
h3 := 3^7*Evaluate(f, x/3);
q2 := q div h3;
R := Resultant(Evaluate(h1, y-x), Evaluate(q1, y));
P := Root(R div q2, 3);

fP := Factorization(P);
g := fP[1][1]; g;
&+[ #Terms(Integers(F) ! c) : c in Coefficients(R) ];

fh := hom< F -> F | -s, t >;
h := hom< Q -> Q | fh, x >;
fminus := f @ h;
gnew := Q ! (x^7*Evaluate(g, (x - 1)*(1 + 2*t/x)));
gnew mod fminus;


U<u> := PolynomialRing(Rationals());
evalst := func< j, k | hom< Q -> U | C, u >
                where C is hom< F -> Rationals() | j, k > >;
for p in [ [1,2], [7,1], [6,-7], [5,4], [1,4], [19,5] ] do
  j, k := Explode(p);
  N1 := NumberField( evalst(j,k)(f) );
  fD := Factorization(Discriminant(Integers(N1)));
  h1 := ClassNumber(N1 : Bound := 300);
  N2 := NumberField( evalst(-j,k)(f) );
  h2 := ClassNumber(N2 : Bound := 300);
  print <p, Min([ h1/h2, h2/h1 ]), fD, Signature(N1)>;
end for;
SetEchoInput(ei);