Intersecting a ray with spheres: GLSL

To the point, here is a dump of some very small code to intersect a ray with N spheres. I used it for my wada rendering so its hardwired to four spheres. It is a function so can be dropped into anyones code. Its dumb and just does a linear search over the spheres. However, you can get reasonable speeds up to 20 spheres or so depending on your card. It is of course ps3.0.

p: the point of origin of the ray.
rd: ray direction.
si: OUTPUT, the closest sphere hit.
t is returned, the distance to the closest intersection point.

Assumption: spheres are passed in as vec4s into a uniform, in the normal GLSL manner. Each vec4 is (x,y,z,r2): centre and radius squared.

uniform vec4 sph[4];  // in this case 4 spheres

float isect (in vec3 p, in vec3 rd, out vec4 si){
float t=999.9,tnow,b,disc;
  for (int i=0; i<4; i++) {  //each sphere
    tnow = 9999.9;           //temporary closest hit is far away
    //next 4 lines are the intersect routine for a sphere
    vec3 sd=sph[i].xyz-p;    
    b = dot ( rd,sd );
    disc = b*b + sph[i].w - dot ( sd,sd );
    if (disc>0.0) tnow = b - sqrt(disc);
    // hit, so compare and store if this is the closest
    if ((tnow>0.0001)&&(t>tnow)) {t=tnow; si=sph[i];}
  }
  return t;
}

A couple of explanations:

disc = b*b + sph[i].w - dot ( sd,sd );

sph[i].w is already a radius squared so saves some code here. Also dot(sd,sd) is short for length(sd) - less bytes.

if ((tnow>0.0001)&&(t>tnow))...

tnow is the current intersection parameter, the distance along rd of the current intersection. If tnow is very small, we may be intersection with the sphere we are already on and must ignore it. 0.0001 is arbitrary.
Of course to get the actual intersection you need to do this in the main routine:

if (t!=999.9) inter=p+t*rd;

or something similar.