1^1 undefined

edit

Consider the following code snippet:

# FACTORIZATION IN F_q[x]
def my_mul(x,y):
    if x == None: return y
    if y == None: return x
    return x*y

# Squarefree Decomposition in F_q[x]
def squarefree_decomposition(GF, f):
    p = GF.characteristic()
    # base case
    if f == 1: return [1]

    # recursive case
    df = f.diff()
    i_factors = [1]
    g = gcd(f,df)
    w = GF(f / g) 
    # collect squarefree factors such that p !| i
    while(w != 1):
        w_aux = gcd(g,w)
        g = GF(g / w_aux)
        i_factor = GF(w / w_aux)
        i_factors.append(i_factor)
        w = w_aux
    # here, g = prod[f_(j)^j for p | j]
    # collect squarefree factors such that p | j
    g_root = GF(g.coefficients(sparse=False)[::p]) # pth root of g
    g_root_decomposition = squarefree_decomposition(GF, g_root)
    j_factors = [g_root_decomposition[i/p] if i%p==0 else 1 for i in xrange(p*(len(g_root_decomposition)-1)+1)]

    # combine results
    res = map(my_mul, i_factors, j_factors)

    return res

# Example in F_q[x]
GF125X = GF(5^3)[x]
f = GF125X((x^5 + x^2 + x^1 + 1)^2*x^5)
decom = squarefree_decomposition(GF125X, f)
print("{} decomposes as {}".format(f, decom))
recom = [f_i^i for f_i,i in enumerate(decom)]
recom

When I execute this, I get the following output:

x^15 + 2*x^12 + 2*x^11 + 2*x^10 + x^9 + 2*x^8 + 3*x^7 + 2*x^6 + x^5 decomposes as [1, 1, x^5 + x^2 + x + 1, 1, 1, x]
Error in lines 28-28
Traceback (most recent call last):
  File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 1188, in execute
    flags=compile_flags) in namespace, locals
  File "", line 1, in <module>
  File "sage/rings/polynomial/polynomial_template.pxi", line 605, in sage.rings.polynomial.polynomial_zz_pex.Polynomial_template.__pow__ (build/cythonized/sage/rings/polynomial/polynomial_zz_pex.cpp:12318)
    raise NotImplementedError("%s^%s not defined."%(ee,self))
NotImplementedError: 1^1 not defined.

I notice I am confused. AFAIK, there is only one sensible result for 1^1, namely 1. Why do I get this error?

If I try creating a 1 polynomial in a finite field and taking its 1th power it works:

id_ = GF125X(1)
id_^1

Produces as output 1, as expected.