In an extension field, how do you to iteratively multiply an irreducible polynomial to obtain its powers, simplifying and printing the results in sagemath?

For example:

1) Within `GF(59^2)['x']`

, the polynomial `x^2 + 2x + 13`

is irreducible.

2) Say we use the label `'g'`

, we can use this irreducible polynomial to express g squared.

```
`g^2 = -2g - 13
= 57g + 46
```

3) Subsequent powers of g can be obtained by iteratively multiplying the previous expression by g.

for example:

```
g^2 = 57g+46
g^3 = g(g^2) = g(57g+46) = 57g^2 + 46g = 57 (57g + 46) + 46g (mod 59)
= 50g + 26,
g^4 = g(g^3) = g(50g + 26) = 50g^2 + 26g = 50(57g+46) + 26g = (50*57)g + (50 * 46) + 26 g (mod 59)
= 44g + 58,
g^5 = g(g^4) = g(44g + 58) = 44g^2 + 58g = 44(57g+46) + 58g = (44*57)g + (44 * 46) + 58g = 30g + 18+ 58g
=29g + 18
(...)
and so on, up to cycle, or g^59.
```

The question is:

```
What would the sagemath code look like to do the following a) initialize g squared, subsequently multiply by g, and simplify mod 59 to print the following output?
```

Sample output:

```
g^2 = 57*g + 46
g^3 = 50*g + 26
g^4 = 46*g + 58
g^5 = 25*g + 51
(....)
g^59 = ....
```

Possible pseudocode:

```
#initialize the ring
F = GF(59^2)['x']
#initialize the starting point to our g^2
g_current = F('57g + 46')
#loop over an appropriate range
for i in range (1,59):
#accumulate the current power, multiplying the previous by g, and then simplifying
g_current = F.simplify(g * g_current)
#print an appropriate line summary
print("g^" + str(i+2) + "=" + str(g_current))
```

Your help is appreciated.