FiniteField
f = Function[x, x^2 + 1]
f = #^2+1 &
IrreduciblePolynomialQ[f[x], Modulus -> 7]
True
F = FiniteField[7, f]
Information[F, "FieldSize"]
49