Hello. I have a finite-field k and its degree n extension E

```
k = GF(4, 'a')
n = 64
E = k.extension(n, 'X')
```

where I perform exponentiation

```
v = E.random_element()
res = v^123456789
```

My goal is to track an element in k^n, mapped to E, as the exponentiation in E occurs. Ideally, I would write it as

```
alphas = var(["alpha_%d" % i for i in range(n)])
v = 0
for i in range(n):
v += alphas[i]*E.gen()^i
res = v^123456789
```

But this leads to "TypeError: positive characteristic not allowed in symbolic computations".

I noticed that polynomials could be used instead of symbolic variables (e.g., https://ask.sagemath.org/question/59697), so another idea is to write it as

```
k = GF(4, 'a')
n = 8
R = PolynomialRing(k, ["alpha_%d" % i for i in range(n)])
P = PolynomialRing(k, 'x')
E = R.extension(P.irreducible_element(n), 'X')
v = 0
for i in range(n):
v += R.gen(i)*E.gen()^i
res = v^65536
```

But this leads to "OverflowError: exponent overflow (65536)". It's also slow for larger values of n.

The next idea is to reduce $ \alpha_i^q = \alpha_i $ using an ideal

```
k = GF(4, 'a')
n = 4
R = PolynomialRing(k, n, 'x')
I = R.ideal([g^(k.order() - 1) - 1 for g in R.gens()])
QR = R.quotient(I, 'alpha')
P = PolynomialRing(k, 'x')
BR = PolynomialRing(QR, 'Y')
E = BR.extension(P.irreducible_element(n), 'X')
v = 0
for i in range(n):
v += QR.gen(i)*E.gen()^i
res = v^65537
```

It seems to work for values of n = 4, but for anything larger it gets stuck on computing primary decomposition for the ideal. See https://ask.sagemath.org/question/59987 for details.

Is there a way to do this in SageMath or, perhaps, using one of the underlying systems, such as Singular, and then move back to SageMath?