I'm trying to unproject the mouse position to get the position on the X-Z plane of a ray cast from the mouse. The camera is fully controllable by the user. Right now, the algorithm I'm using is...

Unproject the mouse into the camera to get the ray:

Vector3 p1 = Vector3.Unproject(new Vector3(x, y, 0), 0, 0, width, height, nearPlane, farPlane, viewProj;
Vector3 p2 = Vector3.Unproject(new Vector3(x, y, 1), 0, 0, width, height, nearPlane, farPlane, viewProj);
Vector3 dir = p2 - p1;
dir.Normalize();
Ray ray = Ray(p1, dir);

Then get the Y-intercept by using algebra:

float t = -ray.Position.Y / ray.Direction.Y;
Vector3 p = ray.Position + t * ray.Direction;

The problem is that the projected position is "jumpy". As I make small adjustments to the mouse position, the projected point moves in strange ways. For example, if I move the mouse one pixel up, it will sometimes move the projected position down, but when I move it a second pixel, the project position will jump back to the mouse's location. The projected location is always close to where it should be, but it does not smoothly follow a moving mouse. The problem intensifies as I zoom the camera out. I believe the problem is caused by numeric instability.

I can make minor improvements to this by doing some computations at double precision, and possibly abusing the fact that floating point calculations are done at 80-bit precision on x86, however before I start micro-optimizing this and getting deep into how the CLR handles floating point, I was wondering if there's an algorithmic change I can do to improve this?

EDIT: A little snooping around in .NET Reflector on SlimDX.dll:

public static Vector3 Unproject(Vector3 vector, float x, float y, float width, float height, float minZ, float maxZ, Matrix worldViewProjection)
{
    Vector3 coordinate = new Vector3();
    Matrix result = new Matrix();
    Matrix.Invert(ref worldViewProjection, out result);
    coordinate.X = (float) ((((vector.X - x) / ((double) width)) * 2.0) - 1.0);
    coordinate.Y = (float) -((((vector.Y - y) / ((double) height)) * 2.0) - 1.0);
    coordinate.Z = (vector.Z - minZ) / (maxZ - minZ);
    TransformCoordinate(ref coordinate, ref result, out coordinate);
    return coordinate;
}

// ...

public static void TransformCoordinate(ref Vector3 coordinate, ref Matrix transformation, out Vector3 result)
{
    Vector3 vector;
    Vector4 vector2 = new Vector4 {
        X = (((coordinate.Y * transformation.M21) + (coordinate.X * transformation.M11)) + (coordinate.Z * transformation.M31)) + transformation.M41,
        Y = (((coordinate.Y * transformation.M22) + (coordinate.X * transformation.M12)) + (coordinate.Z * transformation.M32)) + transformation.M42,
        Z = (((coordinate.Y * transformation.M23) + (coordinate.X * transformation.M13)) + (coordinate.Z * transformation.M33)) + transformation.M43
    };
    float num = (float) (1.0 / ((((transformation.M24 * coordinate.Y) + (transformation.M14 * coordinate.X)) + (coordinate.Z * transformation.M34)) + transformation.M44));
    vector2.W = num;
    vector.X = vector2.X * num;
    vector.Y = vector2.Y * num;
    vector.Z = vector2.Z * num;
    result = vector;
}

...which seems to be a pretty standard method of unprojecting a point from a projection matrix, however this serves to introduce another point of possible instability. Still, I'd like to stick with the SlimDX Unproject routine rather than writing my own unless it's really necessary.