// When are two degree 3 extensions isomorphic?

// Suppose that x1, x2, x3 are the roots of the first degree 3 polynomial
// and x4, x5, x6 are the roots of the second.

// There is a unique conjugacy class of subgroups H of
// S_3 x S_3 with the property that the stabilizer of 1 and the stabilizer
// of 4 are conjugate if and only if the subgroup is contained in H.
// One H that works is generated by (1,2)(4,5) and (1,3,2)(4,6,5).
// A polynomial that's fixed by H is x1x4 + x2x5 + x3x6.

R<x1,x2,x3,x4,x5,x6> := PolynomialRing(Rationals(),6);
theta1 := x1*x4 + x2*x5 + x3*x6;
theta2 := x1*x4 + x2*x6 + x3*x5;
theta3 := x1*x5 + x2*x4 + x3*x6;
theta4 := x1*x5 + x2*x6 + x3*x4;
theta5 := x1*x6 + x2*x4 + x3*x5;
theta6 := x1*x6 + x2*x5 + x3*x4;

F<t> := PolynomialRing(R);
pol1 := (t-x1)*(t-x2)*(t-x3);
pol2 := (t-x4)*(t-x5)*(t-x6);
pol3 := (t-theta1)*(t-theta2)*(t-theta3)*(t-theta4)*(t-theta5)*(t-theta6);

a1 := Coefficient(pol1,2);
a2 := Coefficient(pol1,1);
a3 := Coefficient(pol1,0);
b1 := Coefficient(pol2,2);
b2 := Coefficient(pol2,1);
b3 := Coefficient(pol2,0);

I := IdealWithFixedBasis([a1,a2,a3,b1,b2,b3]);

cof5 := Coefficient(pol3,5);
cof4 := Coefficient(pol3,4);
cof3 := Coefficient(pol3,3);
cof2 := Coefficient(pol3,2);
cof1 := Coefficient(pol3,1);
cof0 := Coefficient(pol3,0);

alist :=
[[a1 ],
[a2, a1^2],
[a3, a1*a2, a1^3],
[a1*a3, a2^2, a2*a1^2, a1^4],
[a3*a2, a3*a1^2, a2^2*a1, a2*a1^3, a1^5 ],
[a3^2, a3*a2*a1, a3*a1^3, a2^3, a2^2*a1^2, a2*a1^4, a1^6]];

AAAA<A1,A2,A3,B1,B2,B3> := PolynomialRing(Rationals(),6);

Alist := [
[A1 ],
[A2, A1^2],
[A3, A1*A2, A1^3],
[A1*A3, A2^2, A2*A1^2, A1^4],
[A3*A2, A3*A1^2, A2^2*A1, A2*A1^3, A1^5 ],
[A3^2, A3*A2*A1, A3*A1^3, A2^3, A2^2*A1^2, A2*A1^4, A1^6]];

blist := [
[b1 ],
[b2, b1^2],
[b3, b1*b2, b1^3],
[b1*b3, b2^2, b2*b1^2, b1^4],
[b3*b2, b3*b1^2, b2^2*b1, b2*b1^3, b1^5 ],
[b3^2, b3*b2*b1, b3*b1^3, b2^3, b2^2*b1^2, b2*b1^4, b1^6]];


Blist := [
[B1 ],
[B2, B1^2],
[B3, B1*B2, B1^3],
[B1*B3, B2^2, B2*B1^2, B1^4],
[B3*B2, B3*B1^2, B2^2*B1, B2*B1^3, B1^5 ],
[B3^2, B3*B2*B1, B3*B1^3, B2^3, B2^2*B1^2, B2*B1^4, B1^6]];

len := [1,2,3,4,5,7];

coflist := [ Coefficient(pol3,5), Coefficient(pol3,4), Coefficient(pol3,3),
Coefficient(pol3,2), Coefficient(pol3,1), Coefficient(pol3,0) ];

for i in [1..6] do
  expolist := [ <x,y,z> : x in [0..i], y in [0..i], z in [0..i] | x+y+z eq i ];
  monlist := [ x1^a[1]*x2^a[2]*x3^a[3]*x4^b[1]*x5^b[2]*x6^b[3] : a, b in expolist ];
  M := ZeroMatrix(Rationals(),len[i]^2,#monlist);
  rowdata := [];
  for aa in [1..len[i]] do
    for bb in [1..len[i]] do
      rownum := aa+(bb-1)*len[i];
      term := alist[i][aa]*blist[i][bb];
      Append(~rowdata,<aa,bb>);
      for j in [1..#monlist] do
        M[rownum][j] := MonomialCoefficient(term,monlist[j]);
      end for;
    end for;
  end for;
  vec := ZeroMatrix(Rationals(),1,#monlist);
  for j in [1..#monlist] do
    vec[1][j] := MonomialCoefficient(coflist[i],monlist[j]);
  end for;
  chk, sol, nul := IsConsistent(M,vec);
  printf "For i = %o, consistent = %o. Dimension nul = %o. Solution = %o.\n",i,chk,Dimension(nul),sol;
  ans := &+[ AAAA!sol[1][j]*Alist[i][rowdata[j][1]]*Blist[i][rowdata[j][2]] : j in [1..NumberOfRows(M)]];
  printf "Formula = %o.\n",ans;
end for; 

// chk5 = -2*A1*B1
// chk4 = A1^2*B1^2 + 2*A1^2*B2 + 2*A2*B1^2 - 6*A2*B2
// chk3 = -2*A1^3*B1*B2 - 2*A1^3*B3 - 2*A1*A2*B1^3 + 5*A1*A2*B1*B2 + 9*A1*A2*B3 
- 2*A3*B1^3 + 9*A3*B1*B2 - 27*A3*B3
// chk2 = 2*A1^4*B1*B3 + A1^4*B2^2 + 3*A1^2*A2*B1^2*B2 - 9*A1^2*A2*B1*B3 - 
6*A1^2*A2*B2^2 + 2*A1*A3*B1^4 - 9*A1*A3*B1^2*B2 + 27*A1*A3*B1*B3 + A2^2*B1^4 - 
6*A2^2*B1^2*B2 + 9*A2^2*B2^2
// chk1 = -2*A1^5*B2*B3 - 2*A1^3*A2*B1^2*B3 - A1^3*A2*B1*B2^2 + 15*A1^3*A2*B2*B3
- 2*A1^2*A3*B1^3*B2 + 9*A1^2*A3*B1*B2^2 - 27*A1^2*A3*B2*B3 - A1*A2^2*B1^3*B2 + 
9*A1*A2^2*B1^2*B3 + 3*A1*A2^2*B1*B2^2 - 27*A1*A2^2*B2*B3 - 2*A2*A3*B1^5 + 
15*A2*A3*B1^3*B2 - 27*A2*A3*B1^2*B3 - 27*A2*A3*B1*B2^2 + 81*A2*A3*B2*B3
// chk0 = A1^6*B3^2 + A1^4*A2*B1*B2*B3 - 9*A1^4*A2*B3^2 - 2*A1^3*A3*B1^3*B3 + 
A1^3*A3*B1^2*B2^2 + 9*A1^3*A3*B1*B2*B3 - 4*A1^3*A3*B2^3 + A1^2*A2^2*B1^3*B3 - 
9*A1^2*A2^2*B1*B2*B3 + A1^2*A2^2*B2^3 + 27*A1^2*A2^2*B3^2 + A1*A2*A3*B1^4*B2 + 
9*A1*A2*A3*B1^3*B3 - 9*A1*A2*A3*B1^2*B2^2 - 27*A1*A2*A3*B1*B2*B3 + 
18*A1*A2*A3*B2^3 - 4*A2^3*B1^3*B3 + A2^3*B1^2*B2^2 + 18*A2^3*B1*B2*B3 - 
4*A2^3*B2^3 - 27*A2^3*B3^2 + A3^2*B1^6 - 9*A3^2*B1^4*B2 + 27*A3^2*B1^2*B2^2 - 
27*A3^2*B2^3