Computing square free part of a multivariate polynomial
In one of the algorithms I am working on, there is a part asking for the square-free part of a multivariate 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? Thank you!