Ray tracing missing cosine

Question

I'm looking at the rendering equation and the source code for a classical ray tracer (so no Monte Carlo or anything yet included), included below in short pseudo code for clarity.

I have trouble understanding why the cosine from the rendering equation isn't included in the ray tracer implementation. I have included a comment on where I think it's missing.

I can see that with the BRDF for perfect specular reflection (f_r = (δ(cosθi − cosθr)δ(φi − φr ± π)) / cosθi) the cosine gets cancelled out. But what about it missing in combination with the phong reflection model as shown below?

for (loop through pixels) 
{        
  Ray ray(eye, getDirection(col+0.5,row+0.5));
  vec3 color = trace(ray, 0);      
}

trace(ray, recursion_depth)
{
  if (too deep, or ray hits nothing)
    return vec3(0,0,0);  // return black

  Hit hit(ray);
  return vec3(hit.Ed)           // light source color
    + classic_trace_lights(hit)         // direct light
    + classic_trace_reflection(hit, depth)     // specular reflection
    + classic_trace_refraction(hit, depth);    // specular refraction
}

classic_trace_reflection(hit, depth)
{
    ....
    vec3 R = ideal_reflected_direction(hit.vdir, hit.normal);
    Ray ray(hit.pos, R);
    return hit.Rs * classic_trace_recursive(ray, depth+1); // Rs = specular color
}

classic_trace_lights(...) { ... }

classic_trace_light(...)
{
    ....
    return phong(...) * light.Ed /* * COSINE MISSING HERE? */;
}
....