K := NumberField(x^4+x+1);
epsilon := K.1^2+1;
K := Parent(epsilon);
Q<u1,u2,u3,u4> := PolynomialRing(K,4);
P  := PolynomialRing(Integers());
thugp := hom<K->P | P.1>;
T  := u1 + u2*K.1 + u3*K.1^2 + u4*K.1^3;
ps := epsilon*T^2;
R  := Monomials(ps); 
S  := Coefficients(ps);
A  := &+[Coefficient(thugp(S[i]),2) * R[i] : i in [1 .. #R]];
B  := &+[Coefficient(thugp(S[i]),2) * R[i] : i in [1 .. #R]];
thugp(S[1]);
