I interpret it the way that the unit vector of the Z axis (in local coordinates) is your forward pointing vector?
The quaternion a describes the translation of the local coordinate frame to the global coordinate frame. Thus a vector v in local coordinates has direction a*v*a' (a'=conjugate of a=inverse for unit quaternion) in the global frame. The direction of the Z axis thus changes according to ka=a*k*a'.
The Z direction of the second object is thus kb=b*k*b'. The cosine of the angle between them is the scalar product of ka and kb which is also the real part of ka'*kb=a*k'*a'*b*k*b'. The real part of any quaternion v remains unchanged in a'*v*a, so
scal(ka,kb) = real( ka * kb' )
= -real( k * a' * b * k * b' * a ) = -real( k * (a' * b) * k * (a'* b)' )
I do not know the specifics of the GL implementation, but the operations are to compute the product p1=a'*b, rotate p1 about the Z-axis by 180 degrees to p2, what in practice is just flipping the signs of the i and j coefficients, and forming the scalar product of p1 and p2. The short version is
scal(ka,kb) = p1.w*p1.w - p1.x*p1.x - p1.y*p1.y + p1.z*p1.z