Fix squarefree_decomposition failure over GF(2)

src/sage/rings/finite_rings/finite_field_base.pyx

@@ -576,11 +576,21 @@ cdef class FiniteField(Field):
             ....:             for j in range(len(F)):
             ....:                 if i == j: continue
             ....:                 assert gcd(F[i][0], F[j][0]) == 1
+
+        Was broken in Sage <= 9.8
+
+            sage: R.<x> = GF(2)[]
+            sage: (x^2 + 1).squarefree_decomposition()
+            (x + 1)^2
+
         """
         from sage.structure.factorization import Factorization
         if f.degree() == 0:
             return Factorization([], unit=f[0])
 
+        def pth_root(x):
+            return x if self.degree() == 1 else x.pth_root()
+
         factors = []
         p = self.characteristic()
         unit = f.leading_coefficient()
@@ -589,7 +599,7 @@ cdef class FiniteField(Field):
         if T0.degree() > 0:
             der = T0.derivative()
             while der.is_zero():
-                T0 = T0.parent()([T0[p*i].pth_root() for i in range(T0.degree()//p + 1)])
+                T0 = T0.parent()([pth_root(T0[p*i]) for i in range(T0.degree()//p + 1)])
                 if T0 == 1:
                     raise RuntimeError
                 der = T0.derivative()
@@ -611,11 +621,11 @@ cdef class FiniteField(Field):
                         if T.degree() == 0:
                             break
                         # T is of the form sum_{i=0}^n t_i X^{pi}
-                        T0 = T0.parent()([T[p*i].pth_root() for i in range(T.degree()//p + 1)])
+                        T0 = T0.parent()([pth_root(T[p*i]) for i in range(T.degree()//p + 1)])
                         der = T0.derivative()
                         e = p*e
                         while der.is_zero():
-                            T0 = T0.parent()([T0[p*i].pth_root() for i in range(T0.degree()//p + 1)])
+                            T0 = T0.parent()([pth_root(T0[p*i]) for i in range(T0.degree()//p + 1)])
                             der = T0.derivative()
                             e = p*e
                         T = T0.gcd(der)