Reducing Polynomial Functions over a Finite Field

Question

Over $\mathbb{F}_{4}$, for example, I am looking for a function that will, for example, reduce the polynomial function (not an element of the finite field itself) $$x^5 + 6x^4 +x^3 + 1 \rightarrow x^2 + x^3 + 1$$

using the identity $x^q = x$ over a finite field of order $q$, and the fact that $p=0$ for the characteristic of the field $p$ (2 in the case above). Again, I am not looking for a way to reduce something modulo some irreducible polynomial. I am just looking for something that will reduce the powers (e.g., $x^5$ to $x^2$ in the example above) and convert any zero coefficients to zero (e.g., $6x^4$ in the example above)

I would also be interested in any potential ways to write such a function. I would need it to work over multiple variables (e.g., $x,y,z,...$) also. Maybe I could use an if statement to check if the power was greater than $q-1$, although with multiple summed terms and variables I wouldn't be completely sure how to do that.

Any advice would be appreciated... Thanks

EDIT: This function seems to work (modified from cvgmt's reply) as well as the one from Carl Woll's reply:

GFReduce[f_, p_, q_, var_] := PolynomialMod[PolynomialMod[f, p] //. var^k_ :> var^Mod[k, q - 1, 1], p];

Answer 1

I think you can use PolynomialReduce for this:

PolynomialReduce[
    x^5 + 6 x^4 + x^3 + 1,
    x^4 - x,
    x,
    Modulus->2
] //Last

1 + x^2 + x^3

Answer 2

Maybe

p=2;
q=4;
PolynomialMod[x^5 + 6 x^4 + x^3 + 1, p] //. x^k_ :> x^Mod[k, q, 1]

1 + x^2 + x^3.

Answer 3

There is a built-in function for this PolynomialRemainder.

Examples with constructed remainder:

Block[{f, d, rem, p = 7, q = 4},
 d = x^q - x;
 rem = x^(q - 1) + 1;
 f = (x^2 + 2 x + 5) d + p (x^4 + 5 x + 2) + rem;
 {PolynomialRemainder[f, d, x, Modulus -> p],
  GFReduce[f, p, q, x],
  rem}
 ]
(*  {1 + x^3, 1 + x^3, 1 + x^3}  *)

Block[{f, d, rem, p = 7, q = 4},
 d = x^q - x;
 rem = PolynomialMod[(x + p - 1)^(q - 2) // Expand, p];
 f = (x^2 + 2 x + 5) d + p (x^4 + 5 x + 2) + rem;
 {PolynomialRemainder[f, d, x, Modulus -> p],
  GFReduce[f, p, q, x],
  rem}
 ]
(*  {1 + 5 x + x^2, 1 + 5 x + x^2, 1 + 5 x + x^2}  *)