Trouble with Phong Shading

I am writing a shader according to the Phong Model. I am trying to implement this equation:

where n is the normal, l is direction to light, v is the direction to the camera, and r is the light reflection. The equations are described in more detail in the Wikipedia article.

As of right now, I am only testing on directional light sources so there is no r^2 falloff. The ambient term is added outside the below function and it works well. The function maxDot3 returns 0 if the dot product is negative, as it usually done in the Phong model.

Here's my code implementing the above equation:

#include "PhongMaterial.h"

PhongMaterial::PhongMaterial(const Vec3f &diffuseColor, const Vec3f &specularColor, 
                            float exponent,const Vec3f &transparentColor, 
                             const Vec3f &reflectiveColor,float indexOfRefraction){

            _diffuseColor = diffuseColor;
            _specularColor = specularColor;
            _exponent = exponent;
            _reflectiveColor = reflectiveColor;
            _transparentColor = transparentColor;


}

Vec3f PhongMaterial::Shade(const Ray &ray, const Hit &hit, 
                    const Vec3f &dirToLight, const Vec3f &lightColor) const{

        Vec3f n,l,v,r;
        float nl;

        l = dirToLight;
        n = hit.getNormal();
        v = -1.0*(hit.getIntersectionPoint() - ray.getOrigin());

        l.Normalize();
        n.Normalize();
        v.Normalize();

        nl = n.maxDot3(l);
        r = 2*nl*(n-l);
        r.Normalize();

        return (_diffuseColor*nl + _specularColor*powf(v.maxDot3(r),_exponent))*lightColor;
}

Unfortunately, the specular term seems to disappear for some reason. My output:

Correct output:

The first sphere only has diffuse and ambient shading. It looks right. The rest have specular terms and produce incorrect results. What is wrong with my implementation?

This line looks wrong:

r = 2*nl*(n-l);

2*nl is a scalar, so this is in the direction of n - l, which is clearly the wrong direction (you also normalize the result, so multiplying by 2*nl does nothing). Consider when n and l point in the same direction. The result r should also be in the same direction but this formula produces the zero vector.

I think your parentheses are misplaced. I believe it should be:

r = (2*nl*n) - l;

We can check this formula on two boundaries easily. When n and l point in the same direction, nl is 1 so the result is also the same vector which is correct. When l is tangent to the surface, nl is zero and the result is -l which is also correct.