How to interpolate vertex normals for Phong lighting when barycentric coordinates are negative?
I have figured out the solution to this, and want to share it. What I am doing is I calculate different barycentric coordinate for each of the step I need to do: z-buffer, interpolation of normal, etc. 

The correct way is to use the barycentric coordinate that is calculated when checking z-buffer to test whether a pixel inside a triangle. Because that barycentric coordinate is guaranteed to be inside the [0, 1] range (otherwise the pixel wouldnt have been drawn), it is the correct barycentric coordinate to use.