Numerical unstability of calculating KL(1-P||1-Q)

Given two raw logit vectors p and q, I would like to calculate -(P-Q)*log((1-P)/(1-Q)), which is equals to KL(1-P||1-Q) + KL(1-Q||1-P), where P/Q is the probability distribution of p/q after softmax().

Here is my code:

-(F.softmax(p, dim=1) - F.softmax(q, dim=1)) * torch.log((1.0001 - F.softmax(p, dim=1)) / (1.0001 - F.softmax(q, dim=1)))

but it turns out to be numerical unstable. How can I fix it?