In one of the algorithms I am working on, there is a part asking for the square-free part of a multivariant polynomial. I can find it using a complicated way, namely,
I = ideal(f)
IRad = I.radical().groebner_basis()
But I think this must be a very inefficient way to do it. Is there a better way?