Why does opengl perform clipping in clip coordinates? Is there any reason why OpenGL performs primitives' clipping in clip coordinates? Can it be at least theoretically done in normalized device coordinates?
The space after clip space is normalized device coordinate space, which is obtained by dividing clip.xyz by clip.w. If W is zero... oops.

Oh sure, you could put a conditional test to see if W is zero or very close to zero. Or you can just do the clipping math in clip space and never have to worry about it. No post-clipping vertex can have a W of zero, because the clip box for each vertex is based on being in the closed range (-W, W). And if W is 0, then the closed range is an empty set, and thus the vertex is not on that range.

Also, negative values of W are outside of the clipping space. The closed range (-W, W) is inverted where W is negative, but the meaning of the range itself is not. Consider a W of -1; the range becomes (1, -1). There are no numbers that are simultaneously greater than 1 and less than -1. Therefore that space is empty.

And yet, division by a negative W can still land NDC points in the [-1, 1] NDC-space range, where they will not be clipped if you didn't clip them beforehand.