This function returns the primary decomposition of the ideal I.
Currently it is implemented ONLY for squarefree monomial ideals
using the Alexander dual technique.
See FrbPrimaryDecomposition for monomial ideals.
Example
/**/ Use R ::= QQ[x,y,z];
/**/ PrimaryDecomposition(***Ideal(xy, yz, zx)***);
[ideal(y, z), ideal(x, z), ideal(x, y)]