A glossy surface has "almost specular" reflections that look slightly blurred or soft.

Here's the usual situation for purely specular reflection, where $V$ points to the viewer and $R_V$ is $V$'s ideal reflection direction around the normal, $N$. $I_\textrm{in}$ is the light intensity coming in, while $I_\textrm{out}$ is the intensity going out.

If a surface is bumpy at the microscopic scale, photons from direction $R_V$ will not necessarily leave in the ideal reflection direction, $V$. The photons will scatter around direction $V$:

Symmetrically, if we want only to measure photons leaving in exactly the direction $V$, these photons could have arrived from several directions around $R_V$:

This is glossy reflection, where the light leaving a surface in a particular direction, $V$, has contributions arriving from many directions around $R_V$.

To simulate this glossy surface, we can cast a number of rays around the direction $R_V$ and recursively compute the incoming intensity from each direction. Each photon is modelled as undergoing a purely specular reflection off of a flat microfacet of the microscopically bumpy surface, so

$I_{\textrm{out},i} = k_s\ I_{\textrm{in},i}$

for the $i^\textrm{th}$ photon. Then the final outgoing intensity toward $V$ due to glossy reflection is the average of the individual intensities (plus the emissive and ambient Phong terms):

$I_\textrm{out} = \left( \frac1N \sum_{i=1}^N I_{\textrm{out},i} \right) + I_\textrm{ambient} + I_\textrm{emissive}$

Using more rays will get a less noisy image, at the cost of more raytracing time.

The glossiness, $g$, determines how spread out the rays are. A $g = 1$ corresponds to purely specular reflection, while $g = 0$ corresponds to purely diffuse reflection.

This $g$ can be considered as the cosine of the half-angle, $\Psi$, of a cone around $R_V$ through which photons arrive:

$g = \cos \Psi$

To generate one random ray that is inside this cone, place a disc of radius 1 centred on $R_V$ and perpendicular to $R_V$. This disc must be at distance $\ell = {1 \over \tan\Psi}$ for it to fit in the cone with half-angle $\Psi$:

Let $u$ and $v$ be two unit vectors that are perpendicular to $R_V$ and to each other (see above).

Let $a$ and $b$ be in the range $\mathbf{[-1,1]}$ with $a^2 + b^2 \leq 1$. Then $au + bv$ is directed from the disc centre to somewhere inside the disc and the point

$p + \ell\ R_V + a\ u + b\ v$

is on the disc.

Finally ... to generate a random ray inside the cone, pick $a$ and $b$ randomly in $[0,1]$ such that $a^2 + b^2 \leq 1$, then raytrace from $p$ through the point above.

The linalg.cpp package has functions perp1() and perp2() that produce the necessary perpendicular vectors, $u$ and $v$. You can use rejection sampling to find an $a$ and $b$ that satisfy the constraint above.