nilrad:= function(R,I)
local rad, radgen, e,i,j,n;
   e:= Elements(R);
   n:= Size(R);
   rad:= I;
   radgen:= ShallowCopy(GeneratorsOfIdeal(I));
   for i in e do
      if not i in rad then
        if ForAny( [2..n], j -> i^j in I) then
           Add(radgen,i);
           rad:= Ideal(R, radgen);
        fi;
      fi;
   od;
   return rad;
 end;