using Nemo
R, x = PolynomialRing(ZZ, "x") S, y = PolynomialRing(R, "y") f = (x^2 + 2x + 1)*y^3 + (x - 1)*y + 4x - 1
f^2