Attach("RelAlgKTheory.m");
Attach("INB.m");

ZX<x> := PolynomialRing( Integers() );
load "D5Mipos.m";
load "D4Mipos.m";
load "D3Mipos.m";
load "Q8pols.m";
load "Q12pols.m";

/* We compute lists of tamely ramified quaternion extensions of order 8 and 12 */
Q8TamePols := [];
for f in Q8pols do
    if IsTamelyRamified(NumberField(f)) then Append(~Q8TamePols, f); end if;
end for;

Q12TamePols := [];
for f in Q12pols do
    if IsTamelyRamified(NumberField(f)) then Append(~Q12TamePols, f); end if;
end for;

/* Given a definig polynomial f, you can run the algorithm as in the following example. */
f := Q12pols[1];
X, isfree := INB( NumberField(f));
L := X[1]; theta := X[2]; Atheta := X[3]; AssOrd := X[4]; nb := X[5]; h := X[6]; QG := X[7];
print "isfree = ", isfree;

/* The following sequence of commands computes  nb as an element of OL  
   and checks if nb really generates OL over AssOrd */
lambda := QG ! ElementToSequence(nb);
delta := QGAction(lambda, theta, h);  /* this is the generating element as an element of OL */
delta;
IsIntegral(delta);
OL := MaximalOrder(L);
delta := OL ! delta;
conj := [ElementToSequence( OL ! QGAction(QG!e, OL ! delta, h)) : e in RowSequence(AssOrd`hnf) ];
Z := Matrix(conj);
Determinant(Z);           /* this determinant must be +-1 */


/* Some tame Q8 examples */
for i:=1 to 4 do
    K := NumberField( Q8TamePols[i] );
    X, isfree := INB(K);
    print "isfree = ", isfree;
end for;

/* Some tame Q12 examples */
for i:=1 to 4 do
    K := NumberField( Q12TamePols[i] );
    X, isfree := INB(K);
    print "isfree = ", isfree;
end for;

/* Some tame D5 examples */
for i:=1 to 4 do
    K := NumberField( D5_pols[i,1] );
    X, isfree := INB(K);
end for;

/* The following are two C_2 x Q_8 extensions */
F := QuadraticField(17);
FX := PolynomialRing(F);
g := Q8TamePols[1];
K := NumberField(FX ! g);
K := AbsoluteField( NumberField(FX ! g) );
IsTamelyRamified(K);
X, isfree := INB( K);  /* stably free and there exists a integral normal basis */

F := QuadraticField(17);
FX := PolynomialRing(F);
g := Q8TamePols[3];
K := NumberField(FX ! g);
K := AbsoluteField( NumberField(FX ! g) );
IsTamelyRamified(K);
X, isfree := INB( K);  /* not stably free  */

/* These are tame Q_16 polynomials */
f := x^16 - 7*x^15 + 120*x^14 - 495*x^13 + 4926*x^12 - 12075*x^11 +
113126*x^10 - 166765*x^9 +
1995351*x^8 - 923310*x^7 + 23515245*x^6 +17244952*x^5 + 150039566*x^4 +
161552045*x^3 + 59691900*x^2 - 741290510*x + 1133509705;

f := x^16 - 5*x^15 + 527*x^14 - 2604*x^13 + 111294*x^12 - 124501*x^11 +
9712704*x^10 + 57807612*x^9 + 310494907*x^8 + 6720634508*x^7 +
30076082564*x^6 +
158189254963*x^5 + 2200512547326*x^4 + 2483004315790*x^3 +
50325786433535*x^2 + 25312978481145*x + 336657798012895;

f := x^16 - 3*x^15 - 809*x^14 + 2281*x^13 + 235556*x^12 - 270219*x^11 -
32780609*x^10 -
48808294*x^9 + 2184241481*x^8 + 10403032006*x^7 - 40205410431*x^6 -
470493093343*x^5 - 1570153883576*x^4 - 2342919383186*x^3 - 1273123854386*x^2 +
163773591455*x + 84937276871;

/* Running the algorithm for Q_16 extensions takes a while, but it works */
X, isfree := INB( NumberField(f));
L := X[1]; theta := X[2]; Atheta := X[3]; AssOrd := X[4]; nb := X[5]; h := X[6]; QG := X[7];
print "isfree = ", isfree;


/* The coefficients for the polynomials used by Cougnard. */
CC := [
[1, -1, 62126, -565081, 1060385071, -16366741325, 465279400700, 7092550941085, 160472449673155],
[1, -1, 62126, -3117631, 465640921, -8688670925, 227961169550, -8096607143015, 78615823872205],
[1, -1, 62126, 455939, 1403958301, 12074281285, 13695081726440, 96273652090225, 45097657035837025],
[1, -1, 62126, -2096611, 913358191, -20639199515, 1748113943690, 13349986875835, 938956719025945],
[1, -1, 62126, -1586101, 1278883351, -30555856265, 10693466710040, -157578154508165, 31010272581373705],
[1, -1, 62126, -2096611, 1084889551, -67412125715, 4827203957690, -123415137229775, 950259835680775],
[1, -1, 62126, -565081, 1219664191, -19648299605, 8370234738980, -262621669734815, 4861757903277895],
[1, -1, 62126, -54571, 813298231, 24859493725, 1797896836850, 130052349272455, 2392810236402205]
                ];

CougnardPols := [ Polynomial(Reverse(p)) : p in CC];

/* In the following we check the 28 extensions given by Cougnard. Each of these
   extensions is stably free, but only the extensions for the pairs 
      (i,j) = (1,2), (1,8), (2,8), (3,5), (3,6), (5,6)
   are free.
   This may take several hours !!!
*/
		
Y := [**];
for i:=1 to #CougnardPols do
    for j:=i+1 to #CougnardPols do
        print "i = ", i, "  j = ", j;
        N1 := NumberField(CougnardPols[i]);
        N2 := NumberField(CougnardPols[j]);
	K := Compositum(N1, N2);
	X, isfree := INB(K);
	Append(~Y, [* i, j, isfree, X *]);
    end for;
end for;




/* The following checks locally free cancellation for all quaternion groups
of order < 32.

L := WithSkewFieldComp(31);
L := [l : l in L | not IsIsomorphic(QuaternionGroup(Integers()!(l[1]/4)), SmallGroup(l[1], l[2]))];

/* L contains the list of small groups which are not generalized
quaternion groups but have a skew field component in the Wedderburn decomposition of QG
*/
L;

C2Q8 := DirectProduct( CyclicGroup(2), QuaternionGroup(2) );
IsIsomorphic(SmallGroup(16,12), C2Q8);

C3Q8 := DirectProduct( CyclicGroup(3), QuaternionGroup(2) );
IsIsomorphic(SmallGroup(24,11), C3Q8);

C2Q12 := DirectProduct( CyclicGroup(2), QuaternionGroup(3) );
IsIsomorphic(SmallGroup(24,7), C2Q12);

/* For Sm(16, 4) the skew field component comes from a projection G ->> Q8 */
G := SmallGroup(16, 4);
Q8 := QuaternionGroup(2);
S := [s`subgroup : s in NormalSubgroups(G) | s`order eq 2];
[IsIsomorphic(quo<G | s>, Q8) : s in S];

/* For Sm(24, 1) the skew field component comes from a projection G ->> Q12 */
G := SmallGroup(24, 1);
Q12 := QuaternionGroup(3);
S := [s`subgroup : s in NormalSubgroups(G) | s`order eq 2];
[IsIsomorphic(quo<G | s>, Q12) : s in S];

/* Sm(24, 3) is the binary tetrahedral group T24 of order 24. For the description of
T24 see Vigneras, page 17, Theorem 3.7.
*/ 
G := SmallGroup(24,3);
D<i,j> := QuaternionAlgebra<Rationals() | -1, -1>;
Delta := MaximalOrder(D);
UDelta, iota := UnitGroup(Delta);
IsIsomorphic(UDelta, G);



/* The following checks the norm condition for all groups of order < 40 */
CheckNormConditionForQuaternionGroups(2,9);

L := WithSkewFieldComp(39);
L;
LL := [];
for v in L do
    T := CharacterTable(SmallGroup(v[1], v[2]));
    if Degree(CharacterField(T[v[3]])) ge 2 then
    Append(~LL, v);
    end if;
end for;
LL;


/* LL is a list of groups of order < 40 which have a skew field component
   which then is necessarily totally definite quaternion and for which
   the character field is not equal to Q. For each of these groups we must
   show that the corresponding Wedderburn components come from a quaternion
   order with the coefficient field unchanged.

   n = 16 Only the quaternion group
   n = 20 Only the quaternion group
   n = 24 Only the quaternion group
   n = 28 Only the quaternion group
   n = 32 There are 4 groups: S[32, 10] --> Q_16
                              S[32, 14] --> Q_16
			      S[32, 20] = Q_32
			      S[32, 41] = C_2 x Q_16
   n = 36 Only the quaternion group
*/

G := SmallGroup(32, 10);
Q16 := QuaternionGroup(4);
S := [s`subgroup : s in NormalSubgroups(G) | s`order eq 2];
[IsIsomorphic(quo<G | s>, Q16) : s in S];

G := SmallGroup(32, 14);
Q16 := QuaternionGroup(4);
S := [s`subgroup : s in NormalSubgroups(G) | s`order eq 2];
[IsIsomorphic(quo<G | s>, Q16) : s in S];

G := SmallGroup(32, 20);
Q32:= QuaternionGroup(8);
IsIsomorphic(G, Q32);

G := SmallGroup(32, 41);
C2Q16:= DirectProduct(CyclicGroup(2), QuaternionGroup(4));
IsIsomorphic(G, C2Q16);

/* It is remarkable that for n = 40 the quaternion group does not satisfy 
   the norm condition for all maximal orders. To be more precise, there
   are 3 maximal orders, 2 of them satisfy the norm condition, the third
   does not.
*/
CheckNormConditionForQuaternionGroups(10,10);


/* The follow sequence shows how one can find fields W for components
which occur in the Wedderburn decomposition of generalized quaternion algebras.
For |Q_n| <= 48 one can always take F(sqrt(-1)).
*/

/* Compute the list of quaternion algebras */
QuatW := ComputeListOfTotDefQuatAlg(2, 12) ;

for i:=1 to #QuatW do
    print "i = ", i;
    D := QuatW[i,1]; Delta := QuatW[i, 3];
    FinRamPlaces, InfRamPlaces := RamifiedPlaces(D);
    F := CoefficientField(D);
    print "Degree(F) = ", Degree(F);
    FX<y> := PolynomialRing(F);
    E := NumberField(y^2 + 1);
    OE := MaximalOrder(E);

    /* Check that the finite ramified places are unramified in E/F and not split */
    for P in FinRamPlaces do
        if Degree(F) eq 1 then
            Factorization(P);
        else
            Factorization(ideal<OE | Basis(P)>);
        end if; 
    end for;

    /* Finally compute an embedding */
    Embed(E, D);
    print " ********************************************************** ";
end for;