Well, hello then!

I need to find intersection between line and sphere in GLSL. Right now my solution is based on Paul Bourke's page and was ported to GLSL this way:

// The line passes through p1 and p2:
vec3 p1 = (...);
vec3 p2 = (...);

// Sphere center is p3, radius is r:
vec3 p3 = (...);
float r = ...;

float x1 = p1.x; float y1 = p1.y; float z1 = p1.z;
float x2 = p2.x; float y2 = p2.y; float z2 = p2.z;
float x3 = p3.x; float y3 = p3.y; float z3 = p3.z;

float dx = x2 - x1;
float dy = y2 - y1;
float dz = z2 - z1;

float a = dx*dx + dy*dy + dz*dz;
float b = 2.0 * (dx * (x1 - x3) + dy * (y1 - y3) + dz * (z1 - z3));
float c = x3*x3 + y3*y3 + z3*z3 + x1*x1 + y1*y1 + z1*z1 - 2.0 * (x3*x1 + y3*y1 + z3*z1) - r*r;

float test = b*b - 4.0*a*c;

if (test >= 0.0) {
  // Hit (according to Treebeard, "a fine hit").
  float u = (-b - sqrt(test)) / (2.0 * a);
  vec3 hitp = p1 + u * (p2 - p1);
  // Now use hitp.
}

It works perfectly! But it seems slow... I'm new at GLSL. You can answer this questions in two ways:

1. Tell me there is no solution, showing some proof or strong evidence.
2. Tell me about GLSL features (vector APIs, primitive operations) that makes the above algorithm faster, showing some example.

Thanks a lot!