Optimal method for computing the Z-Buffer

Question

I am working on a simple 3D software rednerer and I'm just implementing the Z-Buffer for the purpose of hidden surface removal.

I've found many articles about Z-Buffer, but none shows the actual implementation of the Z value calculation, just theory. So I tried to implement the plane equation method as it seems more elegant than the barycentric linear interpolation. These two methods are described here.

So for every face I need to calculate corresponding plane equation in a world space, then convert the coordinates to a clip space and a screen space, and during a rendering phase I convert the screen space coordinate for every pixel back to the world space and calculate the Z value with the plane equation.

The result is pretty good and the hidden surface removal works almost perfect, except I still get some minor Z-fighting sometimes, but that's not the main issure right now.

The main problem is that my depth test creates artefacts on edges with certain angle as shown on the picture below:

I use my own triangle renderer and when I draw two adjacent triangles, there is always one pixel overlay on the edge, because they share two coordinates and draw the same line - that's a standard behaviour, right? Without the depth buffer, the face to be rendered later "wins" the overlay and everything looks fine. With depth testing, it seems that not every pixel on the edge has the same Z value and is fighting to be rendered.

And I have two questions:

Is this method the most optimal regarding time complexity?

And if so, do you have any idea what could cause these artefacts?

Thank you.