### 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 ~~multivariant ~~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!