Blender Cycle's z-coordinate vs. camera-centered z-coordinate (orthogonal distance)
The aforementioned formula is the correct, only that Blender provides the depth at $(x+0.5, y+0.5)$ at pixel $(x, y)$. The general solution is 

$$z_c=\frac{fz_s}{\sqrt{f^2+((c_x-x)^2+(c_y-y)^2)d^2}}$$

respectively with the offset:

$$z_c=\frac{fz_s}{\sqrt{f^2+((c_x-x-0.5)^2d_x^2+(c_y-y-0.5)^2)d^2}}$$

if density in x and y differs:

$$z_c=\frac{fz_s}{\sqrt{f^2+(c_x-x-0.5)^2d_x^2+(c_y-y-0.5)^2d_y^2}}$$