gistfile1.txt

So if your projection matrix P =
 m11 m12 m13 m14
 m21 m22 m23 m24
 m31 m32 m33 m34
 m41 m42 m43 m44
 
 Then
 p_clip = (x_clip,y_clip,z_clip,w_clip)^T = P * (x,y,z,1)^T
 
 We don't actually care about z_clip here, so drop the third row:
 
 P2 = 
  m11 m12 m13 m14
  m21 m22 m23 m24
  m41 m42 m43 m44
  
  The fourth column is now useless (it will contain all zeroes... it was only needed to compute z_clip), so skip that too:
    
 P3 = 
  m11 m12 m14
  m21 m22 m24
  m41 m42 m44
 
Now:
 p_clip = (x_clip, y_clip, w_clip)^T = P * (x,y,z) 
 
 Note that we have clip space xy and w but no z. 
 
 To inverse, do 3x3 matrix inverse and: 
 view_pos = P3^-1 * (x_clip,y_clip,1)
 
 top_left_dir = normalize(P3^1 * (-1,1,1)^T)
 bottom_left_dir = normalize(P3^1 * (-1,-1,1)^T)
 bottom_right_dir = normalize(P3^1 * (1,-1,1)^T)
 top_right_dir = normalize(P3^1 * (1,1,1)^T)
 
 Note that view_pos = (x*w, y*w, z*w) which is a view *ray*. Normally you'd divide by w to get the actual position along the ray, but
 in our case we *want* the ray. So just ignore the divide by w (which we can't do anyway since we don't have the w anymore due to 
 dropping those rows/columns), and just normalize (or don't, if you don't need unit vectors). 
 If you need these things to point to the actual corners on the near plane then scale the vector by near_z/z.
 
 So these are all direction vectors in view space. The origin is obviously (0,0,0) in view space. 
 Rotate direction vector to world space (using inverse view matrix) and use camera position as origin, if you need these in world space.
 
 If this doesn't work, I'd start by validating that your matrix routines are correct (e.g. verify that M*M^-1 = I etc.).