@mrocklin Can you comment on the correctness of this method? I know I should be the one providing it but I find it hard to give a formal argument. 

Also I'm not sure what shall be the result if the function is not invert-able? For example if `f` is `x**2` then the `solns` above are `sqrt(x)` and `-sqrt(x)`. We can return ` Union(*[imageset(f, self.base_set.intersect(imageset(d, f_inv, other))) for f_inv in solns])` but I think only one solution is sufficient.