Thanks for writing this down in details.
I think there are some confusions in your arguments on what is the function that is being differentiated.
In particular, the user function above includes the indexing part and that is what makes all the <0 gradients 0, because the gradient is then 0 * 1/x which is going to be 0 for all non zero x. And you can also see where the nan comes from when x = 0.
I would also add that the chain rule that we use here is only properly defined when each constituent functions are continuously differentiable at the point where you evaluate them.
So even if the full function has well defined behavior, if you use a constituent at a point where it is not continuously differentiable, anything can happen (even though we try to get something sensible). In this case, you do evaluate log at x <=0 so all the theoretical guarantees disappear 
A simple example of that is to consider the identity function (gradient 1 everywhere) and write it as relu(x) - relu(-x). You will see that this gives you a gradient of 0 at 0 and there is nothing we can do about it.

What would be the issue of returning nan ?

Speed mainly. Checking if the input is positive or negative would be 2x the work of the current formula. And setting nans for the negative would end up being 3x slower in total (also use extra memory for the temporary mask/indices).