`is_square` of power series ring fails for rings with zero divisors

### Steps To Reproduce

 ```
sage: R.<x> = PowerSeriesRing(Zmod(16))
sage: ((x+4)^2).is_square()
False
```

### Expected Behavior

returns True

We have a few ways around this. Either implement it, or make it raise `NotImplementedError` by default for non-integral domains. (We probably need to make it raise by default, maybe with `check=False` to let users pretend some `Zmod(p*q)` to be prime.)