// Magma code for the T-nonsplit form
// This works in the online Magma calculator for p up to about 40.

Lower := 3; Upper := 14;
for p in [Lower..Upper] do if IsPrime(p) then
a := PrimitiveRoot(p);
printf("\n==== p is %o, primitive root is %o. ====\n"), p, a;
K<ep> := CyclotomicField(p-1);
Dp := Decomposition(K,p); // Dp; list of primes over p
printf("Divisor of a-ep is %o.\n"),
  [Valuation(a-ep,Dp[i][1]) : i in [1..#Dp]];
RR<x,S,t> := PolynomialRing(K, 3);
P := S*t^(p-1) + p;
F := x^p - S*&+[ (Binomial(p,i) div p)*t^(i-1)*x^i : i in [1..p-1]];
RRbar, pi := quo< RR | P, F >;
nu := func< xx | RR!pi(xx) >;
// nu((1+t*x)^p-1); // Sanity check
sx := ep^(-1)*((1+t*x)^a-1) div t; sx; // image si(x)
si := hom< RR -> RR | sx, S, ep*t >; // si(1+t*x);
Orbx := [x]; for i in [2..p-1] do Append(~Orbx, nu(si(Orbx[i-1]))); end for; // Orbx;
// nu(si(Orbx[p-1])) eq x; // Sanity check
X := [z div MonomialCoefficient(z, t^(p-2)*x^(p-1))] where z is &+Orbx;
printf("First invariant X[1] is: \n %o\n"), X[1];
MonomialCoefficient(X[1], t^(p-2)*x^(p-1));
  for i in [2..p-1] do i;
  z := nu(X[1]*X[i-1]); c := MonomialCoefficient(z, t^(p-i-1)*x^(p-1));
  c;
  [Valuation(c,Dp[i][1]) : i in [1..#Dp]];
  Append(~X, z div c); // X;
  end for;
nu(X[1]*X[p-1]) eq S*X[1];
RRR<y,z,S,t> := PolynomialRing(K, 4);
P := S*t^(p-1) + p;
iy := hom< RR -> RRR | y,S,t >;
Fy := iy(F);
iz := hom< RR -> RRR | z,S,t >;
Fz := iz(F);
de := hom< RR -> RRR | y+z+t*y*z, S, t >;
Y := [iy(m) : m in X];
Z := [iz(m) : m in X];
time A1,A2 := IsDivisibleBy(z,P)
 where z is (p-1)*(de(X[1]) - Y[1] - Z[1])
  -t^(p-1)*&+[Y[i]*Z[p-i] : i in [1..p-1]];
A1;
end if; end for;