// Conjecture A: every C2-extension of Q has class number 1 
//   false for Q(sqrt -5)
// Conjecture B: every C3-extension of Q has class number 1
//   this is the one we are going to test today
// Conjecture C: there are infinitely many number fields of class number 1 
//   still open actually

// How do we generate C3 extensions of Q, i.e. cubic abelian extensions of Q. 
// Most cubic polynomials have Galois group S3 rather than C3, so this needs 
// a bit of thought. For example, we can find them inside Q(zeta_p) because 
// Gal(Q(zeta_p)/Q)=Z/(p-1)Z is cyclic, and has a C3 quotient if (and only if) 
// p=1 mod 3.

// Version 1: Using built in Subfields function in Magma

R<x>:=PolynomialRing(Rationals());  // print polynomials nicer (x instead of $.1)
for p in PrimesUpTo(100) do 
  if p mod 3 ne 1 then continue; end if;  
  K<z>:=CyclotomicField(p);
  S:=Subfields(K,3)[1,1];  
  h:=ClassNumber(S);
  p,h,MinimalPolynomial(S.1);
end for;

// This looks like it takes a while for large primes, 
// and we can time whether it is Subfields or ClassNumber:

R<x>:=PolynomialRing(Rationals());
for p in PrimesUpTo(100) do 
  if p mod 3 ne 1 then continue; end if;  
  K<z>:=CyclotomicField(p);
  time S:=Subfields(K,3)[1,1];  
  time h:=ClassNumber(S);
  p,h,MinimalPolynomial(S.1);
end for;

// Version 2: Find K directly by looking at elements of Q(zeta_p) 
// that are Gal(Q(zeta_p)/K)-invariant

R<x>:=PolynomialRing(Rationals());

K<z>:=CyclotomicField(7);
f:=MinimalPolynomial(z+z^6);    // The maximal totally real subfield 
#GaloisGroup(NumberField(f));   // of Q(zeta_7) is cubic over Q

K<z>:=CyclotomicField(37);
[a: a in [1..36] | IsPower(GF(37)!a,3)];
&+[z^a: a in [1..36] | IsPower(GF(37)!a,3)];
MinimalPolynomial(&+[z^a: a in [1..36] | IsPower(GF(37)!a,3)]);

for p in PrimesUpTo(500) do    // find those p for which Conjecture B fails
  if p mod 3 ne 1 then continue; end if;  
  K<z>:=CyclotomicField(p);
  f:=MinimalPolynomial(&+[z^a: a in [1..p-1] | IsPower(GF(p)!a,3)]);
  S:=NumberField(f);
  h:=ClassNumber(S);
  if h eq 1 then continue; end if;
  p,h,MinimalPolynomial(S.1);
end for;

// Version 3:
// This is still on a slow side because constructing Q(zeta_p) for large p and 
// working in it involves linear algebra in a Q(zeta_p) viewed as a Q-vector 
// space of dimension p-1, which is quite bad for large primes. 
// Let us do it in C instead, and instead of MinimalPolynomial recognize complex 
// numbers and algebraic numbers using LinearRelation (which uses LLL)

R<x>:=PolynomialRing(Rationals());
C<i>:=ComplexField(100);    // C, i
pi:=Pi(C);                  // pi
for p in PrimesUpTo(500) do 
  if p mod 3 ne 1 then continue; end if;  
  z:=Exp(2*pi*i/p);         // zeta_p as a complex number
  r:=&+[z^a: a in [1..p-1] | IsPower(GF(p)!a,3)];
  f:=PolynomialRing(Q)!LinearRelation([1,r,r^2,r^3]);
  if LeadingCoefficient(f) eq -1 then f:=-f; end if;
  assert IsMonic(f);
  p,f;
end for;


// Version 4: Now let us tidy it up by making it into a function, 
// with checks and everything


function AbelianCubic(p)
  assert (Type(p) eq RngIntElt) and (p ge 2) and IsPrime(p) and (p mod 3 eq 1);
  C<i>:=ComplexField(100);
  pi:=Pi(C);
  z:=Exp(2*pi*i/p);
  r:=&+[z^a: a in [1..p-1] | IsPower(GF(p)!a,3)];
  f:=PolynomialRing(Q)!LinearRelation([1,r,r^2,r^3]);
  if LeadingCoefficient(f) eq -1 then f:=-f; end if;
  assert IsMonic(f);
  return NumberField(f);  
end function;


// Print primes p for which the class number is not 1
R<x>:=PolynomialRing(Rationals());
for p in PrimesUpTo(2000) do 
  if p mod 3 ne 1 then continue; end if;  
  S:=AbelianCubic(p);
  h:=ClassNumber(S);
  if h eq 1 then continue; end if;
  p,h,MinimalPolynomial(S.1);
end for;