By shifting and subtracting the bump map, you're finding the directional derivative of the bump map's height function.
By rotating the light vector into tangent space, then using the X and Y components for the shift values, you're finding the component of the perturbed normal vector aligned with the light. In tangent space, the unperturbed normal is a unit vector along the Z axis. When the shifted values are non-zero, they represent the magnitude of the component of the perturbed normal in the direction of the light source. Since the perturbed normal component is parallel to the light source vector (in tangent space), the dot product of this component with the light reduces to a scale operation, which is what a texture map with the texture environment set to modulate does.
Since the perturbed normal is relative to the smooth surface normal, we take the smoothed normal contribution into account when we add in the Gouraud shaded polygon.
There is an assumption that the perturbed normal is not much different from the smoothed surface unit normal, so that the length of the perturbed normal is not much different from one. If this assumption wasn't true, we'd have to create and modulate in an extra texture that would renormalize the perturbed normal. This can be done, at the cost of an extra texturing pass, if more accuracy is needed.