confusion about gluProject and OpenGL clipping

Consider the fixed transformation pipeline of OpenGL, with the following parameters:

GL_MODELVIEW_MATRIX
0,175,303.109,0,688.503,-2741.84,1583,0,29.3148,5.52094,-3.18752,0,-87.4871,731.309,-1576.92,1
GL_PROJECTION_MATRIX
2.09928,0,0,0,0,3.73205,0,0,0,0,-1.00658,-1,0,0,-43.9314,0
GL_VIEWPORT
0,0,1920,1080

When I draw the faces of the unit cube I get the following:

By looking at the picture, I would expect half of the vertices to have pixel y-coordinate above 1080, and the other half to have a negative y-coordinate. Instead, with gluProject, all vertices have y > 1080:

model coordinate 0 0 0 -> screen coordinate 848.191 1474.61 0.989359
model coordinate 1 0 0 -> screen coordinate 821.586 1973.88 0.986045
model coordinate 0 1 0 -> screen coordinate -198317 667165 4.61719
model coordinate 1 1 0 -> screen coordinate -2957.48 12504.1 1.07433
model coordinate 0 0 1 -> screen coordinate 885.806 1479.77 0.989388
model coordinate 1 0 1 -> screen coordinate 868.195 1979.01 0.986088
model coordinate 0 1 1 -> screen coordinate -438501 1.39841e+06 8.60228
model coordinate 1 1 1 -> screen coordinate -3191.35 12592.4 1.07507

I could successfully reproduce the gluProject() results with my "custom" calculations.

Why the y-coordinate of all vertices is above 1080?

P.S. To draw the cube I rely on:

glBegin(GL_QUADS);
for(int d = 0; d < 3; ++d)
    for(int s = 0; s < 2; ++s)
        for(int v = 0; v < 4; ++v)
        {
            const int a = (v & 1) ^ (v >> 1 & 1);
            const int b = v >> 1 & 1;

            const int d0 = d;
            const int d1 = (d + 1) % 3;
            const int d2 = (d + 2) % 3;

            double p[3];

            p[d] = s;
            p[d1] = a;
            p[d2] = b;

            glColor3dv(p);
            glVertex3dv(p);
        }
glEnd();

I found the answer, in part thanks to this post.

The explanation is that the 4 vertices that have y < 0 in screen space, are also behind the camera, and so have w_clip < 0.

Perspective division (y_clip/w_clip) produces in turn a positive value in device independent coordinates and screen space.