Ok, first thing is to split out the channels with the separate RGB node. We'll put down a few re-routes too so we don't go nuts and mix channels up:

Next, put in a Combine XYZ node to recompile the channels together. (note that you can use a Combine RGB node here instead of you like, they are identical aside from labels and the color of the output socket. Color and Vector and both a generic vector3 format inside Cycles, they are only differentiated for UI clarity reasons).

Now, currently we have Z=0. That's not going to work. As you noted in the OP, we need to derive Z from X and Y. The reason for this is that X and Y have already been normalized so the full XYZ vector is in a 0-1 range. Thus we can't just shove a 1.0 in there, we need to actually compute Z using the forumla you have above. (note that this setup is exactly what the Unreal 4 engine DeriveNormalZ does under the hood).
First the multiplications within the parentheses, X * X and Y * Y:

Then sum the result:

Subtract that from 1:

Now we need the square root of that result. That's the same as raising to the power of 0.5, so one more math node:

Finally, the particular normal map you provided appears to be DirectX formatted, so we need to flip the G(Y) channel too:

And it works!
