I think this can be seen as a matter of how you define the divergence angle $\theta$. Different authors seem to define it slightly differently, however in practice this is rarely an issue since $\theta$ is typically small enough that $2\tan(\theta/2) \approx \tan(\theta) \approx \theta$.
The Gaussian beam as a solution to the paraxial Helmholtz equation has the intensity
$$I(\rho, z) = I_0 \left[\frac{W_0}{W(z)}\right]^2 \exp \left[-\frac{2\rho^2}{W^2(z)}\right],$$
$$W(z) = W_0 \sqrt{1 + (z/z_0)^2},$$
$$W_0 = \sqrt{\frac{\lambda z_0}{\pi}},$$
where $\rho$ is the distance from the beam axis, $z$ is the position along the beam axis, $W_0$ is called the waist radius and $z_0$ the Rayleigh range (from Saleh & Teich). The fraction of the power of the beam contained by a circle of radius $\rho$ about the beam that is a distance $z$ away from the beam waist is
$$1 - \exp \left[-\frac{2\rho^2}{W^2(z)}\right].$$
If you define the beam diameter (as is usually done) as the diameter containing $1 - e^{-2}$ of the power, then the beam diameter is
$$D(z) = 2W(z) = 2W_0 \sqrt{1 + (z/z_0)^2} = D_0 \sqrt{1 + (z/z_0)^2}$$
with $D_0 = 2W_0$. Far away from the beam waist ($z \gg z_0$), we have
$$D(z) \approx D_0 \frac{z}{z_0}.$$
This diameter is defined by a cone along the beam axis with the vertex at the center of the beam. The half-angle of this cone is
$$ \theta_0 = \tan^{-1}\left(\frac{W_0}{z_0}\right) $$
Now, if we define the divergence to be the full vertex angle of this cone, i.e. $\theta = 2\theta_0$,
$$ D(z) = \sqrt{D_0^2 + \left(2\frac{W_0}{z_0}z\right)^2} = \sqrt{D_0^2 + \left(2z\tan\theta_0\right)^2} = \sqrt{D_0^2 + \left[2z\tan(\theta/2)\right]^2} $$
If we define the divergence as $\theta = \tan^{-1}(D_0/z_0)$ (although I don't see an obvious reason for doing so),
$$ D(z) = \sqrt{D_0^2 + \left(\frac{D_0}{z_0}z\right)^2} = \sqrt{D_0^2 + \left(z\tan\theta\right)^2} $$
If we define the divergence as $\theta = D_0/z_0$ so that $D(z) \approx \theta z$ for $z\gg z_0$ (this is twice the definition in Saleh & Teich),
$$ D(z) = \sqrt{D_0^2 + \left(\frac{D_0}{z_0}z\right)^2} = \sqrt{D_0^2 + \left(\theta z\right)^2}. $$
I believe the formulas you gave are approximations for these expressions that are valid for $z \approx 0$ and $z \gg z_0$.