
// The Magma script below tests condition (c) of Theorem 2 for
// K=\Q(\zeta_n)^+ with n<100 satisfying n \equiv 2 mod 4
// and 5 \nmid n. This computation is needed for the 
// proof of Corollary 1.1.



// units is a sequence of units in a number field K
// that generates a subgroup of the unit group of finite index.
// Returns unitsNew, which is a sequence of independent units
// that generates the same group modulo torsion. 

basis:=function(units);
	R:=Parent(units[1]);
	OK:=MaximalOrder(R);
	K:=NumberField(OK);
	r:=#units;
	s,t:=Signature(K);
	A:=FreeAbelianGroup(r);
	amb:=StandardLattice(r);	
	ker:=amb;
	lmt:=50;
	repeat
		prms:=PrimesUpTo(lmt,K);
		lmt:=2*lmt;
	until #prms ge 3*r; 
	for P in prms do
		RP,phiP:=ResidueClassField(P);
		RPs,piP:=MultiplicativeGroup(RP);
		imP:=[phiP(u) : u in units];
		assert &and[a ne RP!0 : a in imP];
		imP:=[a@@piP : a in imP];
		psi:=hom<A->RPs | imP>;
		ker:=ker meet sub<amb | [amb!Eltseq(A!k) : k in Generators(Kernel(psi))]>;
		ker:=sub<amb | [amb!k : k in Basis(ker)]>;
	end for;
	kermat:=Rows(LLLBasisMatrix(ker));
	kermat:=[amb!k : k in kermat];
	ker:=[Eltseq(k) : k in kermat | Length(k) lt 100];
	ker:=[k : k in ker | &*[OK | units[i]^k[i] : i in [1..r]] eq OK!1];
	if #ker ge 1 then
		M:=Matrix(ker);
		N,U,V:=SmithForm(M);
		Vin:=V^(-1);
		units:=[&*[OK | units[i]^v[i] : i in [1..r]] : v in Rows(Vin)];
		D:=Diagonal(N);
		j:=Max([i : i in [1..#D] | D[i] ne 0]);
		newUnits:=[units[i] : i in [(j+1)..r]];
		assert #newUnits eq s+t-1;
		return newUnits;
	else
		assert #units eq s+t-1;
		return units;
	end if;		
end function;


for n in [7..99] do
	if n mod 4 ne 2 and n mod 5 ne 0 then

	print "+++++++++++++++";
	print n;

	L<zet>:=CyclotomicField(n);
	OL:=MaximalOrder(L);
	K:=sub<L| zet+1/zet>;  // K=Q(\zeta_n)^+.
	OK:=MaximalOrder(K);
	if IsOdd(n) then // zd generates the roots of unity in L
			// and m is the order of zd.
		zd:=-zet;
		m:=2*n;
	else
		zd:=zet;
		m:=n;
	end if;
	// We construct the group of cyclotomic units.
	V:=[zd] cat [ 1-zet^a : a in [1..(n-1)]];
	// The cyclotomic units in L are the subgroups of units
	// in the multiplicative group generated by V.
	A:=AbelianGroup([m] cat [0 : i in [1..(n-1)]]); // Think of A as the subgroup of L^*
	// with basis V.

	// The units inside this subgroup generated by V are the elements of valuation 0
	// at all primes. The only primes that possible divide the elements of V are those
	// that divide $n$.
	ker:=A;
	facts:=Factorisation(n*OL);
	facts:=[P[1] : P in facts];
	for P in facts do
		phi:=hom<A ->FreeAbelianGroup(1) | [Valuation(u,P) : u in V]>;
		ker:=ker meet Kernel(phi);
	end for;  // Now ker is a subgroup of A corresponds to the units inside the subgroup generated by V.
	
	lmt:=50;
	repeat
		prms:=PrimesUpTo(lmt,L);
		lmt:=2*lmt;
	until #prms ge 3*#V; 

	for P in prms do
		if Valuation(n*OL,P) eq 0 then
			RP,phiP:=ResidueClassField(P);
			RPs,piP:=MultiplicativeGroup(RP);
			mu:=hom<A -> RPs | [ (phiP(v)/phiP(ComplexConjugate(v)))@@piP : v in V  ] >;
			ker:=ker meet Kernel(mu);
		end if;
	end for;	

	kerbas:=[Eltseq(A!k) : k in OrderedGenerators(ker)];
	cycBas:=[K!(&*[ V[i]^k[i] : i in [1..n]]) : k in kerbas]; 
	// The fact that this coercion works means that the elements
	// inside ker are all real.
	// cycBas generates the group of real cyclotomic units.

	cycBas:=basis(cycBas); // Independent system of 
	// cyclotomic units modulo torsion.

	cycBas:=[-1] cat cycBas;

	// We want to compute generators for the totally positive subgroup 
	// of this.

	B:=AbelianGroup([2] cat [0 : i in [1..#cycBas-1]]);
	C:=AbelianGroup([2 : i in [1..Degree(K)]]);

	psi:=function(u);
        	e:=RealEmbeddings(u);
	        v:=[0 : i in [1..(Degree(K))]];
        	for i in [1..Degree(K)] do
                	if e[i] lt 0 then
                        	v[i]:=1;
           	     end if;
       		end for;
        	return C!v;
	end function;

	mu:=hom<B->C | [psi(u) : u in cycBas]>;
	Mbas:=[[Integers()!r : r in Eltseq(B!b)] : b in OrderedGenerators(Kernel(mu))];
	cycPos:=[ &*[cycBas[i]^w[i] : i in [1..#cycBas]] : w in Mbas   ]; // These are
								// generators
								// for the 
							// totally positive 
							// cyclotomic units. 
	
	// In fact for the claim in the proof of Corollary 1.1 all 
	// we need is that the elements of cycPos are totally positive
	// units in K and we now double-check that.
	assert &and[IsTotallyPositive(m) : m in cycPos];
	assert &and[Norm(m) in {1,-1} : m in cycPos];
	assert &and[m in K : m in cycPos];

	facts:=Factorisation(5*OK);
	S5:={P[1] : P in facts};

	// Given a set S contained in S_5,
	// returns the product in (3) of the paper.
	norm:=function(u,S);
		assert S subset S5;
        	N:=1;
	        for P in S do
        	        F,pi:=ResidueClassField(P);
                	N:=N*Norm(pi(u));
	        end for;
        	return N;
	end function;

	if #S5 eq 1 then
        	print "true"; // In this case there are no non-empty proper subsets of S5
				// and so condition (c) of Theorem 2 is vacuously true.
	else
		tf:=true;
	        for i in [1..(#S5 -1)] do
        	        for S in Subsets(S5,i) do // S is a non-empty proper subset of S5 as the number of elements is
						// is between 1 and #S5-1.
                        	tf:= tf and  &or[ norm(u,S) ne 1 : u in cycPos  ];
           	     end for;
        	end for;
		print tf; // True means that condition (c) of Theorem 2 is true.
	end if;
	end if;
end for;


