Firstly, extracting from just the projection matrix is a good place to start debugging.
To answer your question:
In this case, your near and far plane values should be <0, 0, 1, -1> and <0, 0, -1, 1000> respectively. In practice, it might be off by a little due to floating point rounding errors.
Now to help you solve your larger problem:
You say you are getting a near and far plane of <0, 0, -1, -0.3328891> and <0, 0, -1, 1> respectively. These values comes directly from summing/subtracting rows or columns from your projection matrix, so writing out the projection matrix values on paper and eyeballing it might make it obvious.
Make sure your projection matrix is correct.
Make sure your plane normalization routine is correct (<a / |n|, b / |n|, c / |n|, d / |n|>, where |n| is the length of <a, b, c> ONLY).
Make sure you are not getting your rows and columns mixed up.
You are basically going to need to comb through your code step-by-step and make sure everything has been done correctly.
I recently did this (for Direct3D though) and had to go through the same steps.