Go-GL "Project" method giving unexpected results

Go-GL's Project method is giving me unexpectedly large screen coordinates.

To summarize:

// Screen is 800x600.
projection := mgl32.Perspective(
    mgl32.DegToRad(45), // Field of view (45 degrees).
    800.0 / 600.0,      // Aspect ratio.
    0.1,                // Near Z at 0.1.
    10)                 // Far Z at 10.
camera := mgl32.LookAtV(
    mgl32.Vec3{0, 0.1, 10}, // Camera out on Z and slightly above.
    mgl32.Vec3{0, 0, 0},    // Looking at the origin.
    mgl32.Vec3{0, 1, 0}     // Up is positive Y.
model := mgl32.Ident4()     // Simple model matrix, to avoid confusion.
modelView := camera.Mul4(model) // The model-view matrix (== camera, here).

// Okay, so what does the origin translate to? Expect center-of-screen, with arbitrary-seeming depth.
origin := mgl32.Vec3{0, 0, 0}
screenOrigin := mgl32.Project(origin, modelView, projection, 0, 0, 800, 600)
fmt.Printf("Origin: (%v, %v, %v)", screenOrigin[0], screenOrigin[1], screenOrigin[2])

// What about the point 5 to the right of the origin?
// Expect increased X, but still less than screenWidth.
// Y and Z should be the same as for the origin.
// In my actual program, I drew (-1,-1,-1)-(1,1,1) cube at (5,0,0) and
// it is completely visible in the window.
test := mgl32.Vec3{5, 0, 0}
screenTest := mgl32.Project(test, modelView, projection, 0, 0, 800, 600)
fmt.Printf("Test:   (%v, %v, %v)", screenTest[0], screenTest[1], screenTest[2])

The results?

Origin: (400, 300, 5.500255)
Test:   (4021.32, 300, 5.500255)

4021.32? That's WAY off the screen!

I've pretty much eliminated all the variables I can think of except for the cube rendering I'm using as a hint basicall. My code is based on this, but I moved the camera and the cube: https://github.com/go-gl/examples/tree/master/gl41core-cube

I get the same results if I use (0, 0, 0) as the input vector and translate Ident4() by 5 units in +X (which makes sense).

So what am I doing wrong? Based on the position of the cube I'm drawing, I would expect something like X = 625.

Out of curiosity I looked into the source code of mgl32. It seems that the project method is simply wrong:

obj4 := obj.Vec4(1)

vpp := projection.Mul4(modelview).Mul4x1(obj4)
win[0] = float32(initialX) + (float32(width)*(vpp[0]+1))/2
win[1] = float32(initialY) + (float32(height)*(vpp[1]+1))/2
win[2] = (vpp[2] + 1) / 2

As one might notice, the perspective divide is missing after the matrix multiplication, thus the function will not work whenever w != 1.0f.

In the first example, one cannot notice the error since x and y are zero after the matrix multiplication, thus dividing by w does not make a difference (except for z). In the second example, the x coordinate is not zero anymore and the bug gets visible.