I'm trying to make a software z-buffer implementation, however, after I generate the z-buffer and proceed with the vertex culling, I get pretty severe discrepancies between the vertex depth and the depth of the buffer at their projected coordinates on the screen (i.e. zbuffer[v.xp][v.yp] != v.z, where xp and yp are the projected x and y coordinates of the vertex v), sometimes by a small fraction of a unit and sometimes by 2 or 3 units. Here's what I think is happening:

Each triangle's data structure holds the plane's (that is defined by the triangle) coefficients (a, b, c, d) computed from its three vertices from their normal:

void computeNormal(Vertex *v1, Vertex *v2, Vertex *v3, double *a, double *b, double *c) {
  double a1 = v1 -> x - v2 -> x;
  double a2 = v1 -> y - v2 -> y;
  double a3 = v1 -> z - v2 -> z;

  double b1 = v3 -> x - v2 -> x;
  double b2 = v3 -> y - v2 -> y;
  double b3 = v3 -> z - v2 -> z;

  *a = a2*b3 - a3*b2;
  *b = -(a1*b3 - a3*b1);
  *c = a1*b2 - a2*b1;
}

void computePlane(Poly *p) {
  double x = p -> verts[0] -> x;
  double y = p -> verts[0] -> y;
  double z = p -> verts[0] -> z;

  computeNormal(p -> verts[0], p -> verts[1], p -> verts[2], &p -> a, &p -> b, &p -> c);

  p -> d = p -> a * x + p -> b * y + p -> c * z;
}

The z-buffer just holds the smallest depth at the respective xy coordinate by somewhat casting rays to the polygon (I haven't quite got interpolation right yet so I'm using this slower method until I do) and determining the z coordinate from the reversed perspective projection formulas (which I got from here:

double z = -(b*Ez*y + a*Ez*x - d*Ez)/(b*y + a*x + c*Ez - b*Ey - a*Ex);

Where x and y are the pixel's coordinates on the screen; a, b, c, and d are the planes coefficients; Ex, Ey, and Ez are the eye's (camera's) coordinates.

This last formula does not accurately give the exact vertices' z coordinate at their projected x and y coordinates on the screen, probably because of some floating point inaccuracy (i.e. I've seen it return something like 3.001 when the vertex's z-coordinate was actually 2.998).

Here is the portion of code that hides the vertices that shouldn't be visible:

for(i = 0; i < shape.nverts; ++i) {
  double dist = shape.verts[i].z;
  if(z_buffer[shape.verts[i].yp][shape.verts[i].xp].z < dist)
    shape.verts[i].visible = 0;
  else
    shape.verts[i].visible = 1;
}

How do I solve this issue?

EDIT

I've implemented the near and far planes of the frustum, with 24 bit accuracy, and now I have some questions:

1. Is this what I have to do this in order to resolve the flickering?
2. When I compare the z value of the vertex with the z value in the buffer, do I have to convert the z value of the vertex to z' using the formula, or do I convert the value in the buffer back to the original z, and how do I do that?
3. What are some decent values for near and far?

Thanks in advance.