Getting tight binding density of states more accurately I calculated numerically the density of states (DoS) for the 3-D tightbinding dispersion $\epsilon(k_x,k_y,k_z)=-2t\,(\cos k_x + \cos k_y + \cos k_z)$ and obtained the following plot [$t=1$ has been chosen]. 

What I did is summing over $k$-points of the lattice Green's function, $$G(k_x,k_y,k_z,\omega)=\frac{1}{\omega-\epsilon(k_x,k_y,k_z)+i0^+}$$ and finding the DoS from its imaginary part : $D(\omega)=-\frac{1}{\pi}\text{Im} \sum_{k_x}\sum_{k_y}\sum_{k_y} G(k_x,k_y,k_z,\omega)$.
One can easily notice that there are noises at low energies. Is there any alternative way to get better result? Like the one shown in a figure from a paper [Ref: arXiv:1207.4014] :

Can there be some mathematical standard expression that can be calculated through Mathematica or Matlab? 
Related bonus question : Can the same method be applied to an asymmetric triangular lattice having dispersion $\epsilon(k_x,k_y)=-2t\,(\cos k_x + \cos k_y)-2t'\,\cos(k_x+k_y)$ ?
 A: TL;DR
With Matlab or Mathematika I cannot help, but
in Python there is an implementation available: sc_dos

Here $D$ is the half-bandwidth $D = 6t$.
import numpy as np
import gftool as gt

eps = np.linspace(-1.2, 1.2, num=6001)
dos = gt.sc_dos(eps, half_bandwidth=1)

The evaluation of the DOS takes for me ~100 ms.

You already gave the correct expressions. We have the Green's function
$$
G(z) = \frac{1}{N} \sum_{\boldsymbol{k}} \frac{1}{z-\epsilon_{\boldsymbol{k}}}
$$
and the (normalized) density of states (DOS)
$$
D(\epsilon)
 = \frac{1}{N} \sum_{\boldsymbol{k}} \delta(\epsilon - \epsilon_{\boldsymbol{k}})
= -\frac{1}{\pi} \Im G(\epsilon+i0^+),
$$
where $\epsilon$ is a real energy variable.
The second equality is Sokhotski–Plemelj.
A naive sum over the points is extremely demanding, as a tremendous number of $\boldsymbol{k}$ points is necessary in 3 Dimensions.
To smoothen the function, we can evaluate the Green's function on a contour parallel to the real axis shifted by a finite $\eta>0$ into to upper complex half-plane:
$$
D(\epsilon) \approx -\frac{1}{\pi} \Im G(\epsilon+i\eta).
$$
The bigger we choose $\eta$ the smoother the function becomes, but on the other hand we loose features.
As we are only interested in the thermodynamic limit $N\rightarrow \infty$,
a smarter approach than just sampling $\boldsymbol{k}$, is to replace the sum by the integral. For integrals, we have more or less efficient algorithms.
So let's calculate
$$
G(z) = \int \frac{d^3 k}{{2\pi}^3} \frac{1}{z-\epsilon_{\boldsymbol{k}}}
$$
instead. The $\epsilon_\boldsymbol{k}$, is symmetric for all $k_{x_i}$: $\epsilon(k_{x_i}) = \epsilon(-k_{x_i})$, thus it is enough to integrate over a eighth of the Brillouin zone.
And finally we can use analytic results for the integrals. We note that we can express the 3D Green's function in terms of known results of the 1D and 2D Green's function as we have
$$
\epsilon_{\boldsymbol{k}} = \epsilon_{k_x, k_y, k_z}
= \epsilon^{2D}_{k_x, k_y} - 2t \cos(k_z)
= \epsilon^{1D}_{k_x} - 2t[\cos(k_y) + \cos(k_z)]
$$
and therefore
$$
G(z)
 =  \int \frac{dk_z}{2\pi} G^{2D} (z - 2t\cos(k_z))
 = \int \frac{dk_y}{2\pi} \int \frac{dk_z}{2\pi} G^{1D}(z - 2t[\cos(k_y) + \cos(k_z)]).
$$
The one-dimensional Green's function $G^{1D}(z)$ can be easily evaluated, the two-dimensional Green's function $G^{2D}(z)$ can be expressed in terms of the complete elliptic integral of first kind (which can be found in standard text-books). Using $G^{2D}(z)$ is basically the result given by bRost03.
A very smart guy named Joyce even found an expression for $G^{3D}(z)$ in 1973.
The equations are a bit lengthy and complicated, so I will avoid copying them here. But we implemented them in a Python module gftool>=0.8.0, see sc_dos. You will also find the relevant references there.
A: Sorry for the late response but hopefully this can be useful for someone else!
You can reduce the noise using an elliptic integral.
$$D(\varepsilon)=\frac{1}{4 \pi ^3 t}\int_{-\pi }^{\pi }d\phi K\left( \sqrt{1-\left(\frac{\varepsilon +2 t \cos\phi }{4 t}\right)^2}\right)$$ 
Where K is the Complete Elliptic Integral of the First Kind: http://mathworld.wolfram.com/CompleteEllipticIntegraloftheFirstKind.html.
It's not trivial to get here. And even from this expression the integral needs to be done numerically with care (it has singularities for many values of $\varepsilon$) but it should give better results. Running for six seconds in Mathematica gives me (with $t=\frac{1}{2\sqrt{3}}$):

A: I had the same trouble. I used the formula $\rho(\epsilon)=\sum_{\vec{k}}\delta(\epsilon - \epsilon(\vec{k}))$ to numerically calculate the density of states. I did a summation over $100$ $k-$values for each component, and used a Gaussian distribution with $\sigma =0.1$ for the delta function to get the following diagram. Using larger $\sigma$ ends up smoothening out the singularity in the derivative at around $\epsilon=\pm 2$.
The code was written in C++, and ran for about $60$ seconds. $y$ axis is $\rho(\epsilon)$, $x$ axis is $\epsilon$ and $t=1$.


PS: I did the summation over half the Brillouin zone, which was what I needed for my application.
