MAINT: use scipy.special.xlogy

This should fix the edge cases of the domain of `entr` and `kl_div` so that it allows zero in some places that were not allowed before this PR, and so that it improves numerical stability near zero.  It still needs unit tests and it needs to return the correct values (-inf or inf or whatever) outside the domain but I hope these could be added by someone else after this PR is merged.