hmmm...maybe sifted could use a transform option so one could do

```
sifted = sift(a, lambda arg: arg[1] in b1, lambda arg: arg[0]*arg[1]**S.Half)
p = _mexpand([Add._from_args(sifted[True])**2 - Add._from_args(sifted[False])**2)  # I think _mexpand will work like this
```

This also might be a good place to use `Add._from_args` since if there is only one term it will be correct and if there are more they will be reordered by the squaring and expanding.