Here's something off the top of my head, trying to minimize expensive operations, without resorting to approximation hacks (ie. with infinite precision real numbers, the algorithm below is exactly correct, although in practice finite precision will introduce numerical errors). Recommendations/edits to improve performance are welcome.
Let s be the start point of the ray, and d a unit vector in the direction of the ray.
Let c be the center point of the sphere, and r its radius.
For simplicity, I'll assume that you only want points on the ray where it enters or kisses the sphere, forward from the start point. Intersections behind the start point, or exiting the sphere, are ignored (this means a ray originating inside the sphere detects no collision).
// Calculate ray start's offset from the sphere center
float3 p = s - c;
float rSquared = r * r;
float p_d = dot(p, d);
// The sphere is behind or surrounding the start point.
if(p_d > 0 || dot(p, p) < rSquared)
// Flatten p into the plane passing through c perpendicular to the ray.
// This gives the closest approach of the ray to the center.
float3 a = p - p_d * d;
float aSquared = dot(a, a);
// Closest approach is outside the sphere.
if(aSquared > rSquared)
// Calculate distance from plane where ray enters/exits the sphere.
float h = sqrt(rSquared - aSquared);
// Calculate intersection point relative to sphere center.
float3 i = a - h * d;
float3 intersection = c + i;
float3 normal = i/r;
// We've taken a shortcut here to avoid a second square root.
// Note numerical errors can make the normal have length slightly different from 1.
// If you need higher precision, you may need to perform a conventional normalization.
return (intersection, normal);