In order to calculate the density of states of single particle in the simple harmonic potential, we would calculate that $$ D(\epsilon)=\sum_{n}\delta(\epsilon-\epsilon_n) $$ where $\epsilon_n=(n+1/2)\hbar\omega$. In the limit $\hbar\omega\ll1$,we find that $$ D(\epsilon)\approx\frac{1}{\hbar\omega}\theta(\epsilon). $$
But I want to know how to calculate the $D(\epsilon)$ exactly by means of some special function for example.
There is a formal manipulation that can answer your question using well-known formulas. So, let me write $$ \epsilon_n=(2n+1)\epsilon_0 $$ being $\epsilon_0=\hbar\omega/2$ and $$ \delta(\epsilon-\epsilon_n)=\int_{-\infty}^{+\infty}\frac{dt}{2\pi}e^{-i(\epsilon-\epsilon_n)t}. $$ Then, $$ D(\epsilon)=\int_{-\infty}^{+\infty}\frac{dt}{2\pi}e^{-i\epsilon t} \sum_{n=0}^\infty e^{i(2n+1)\epsilon_0 t} $$ where I have formally exchanged the sum with the integral (just one of my yet-to-be-justified mathematical steps). The sum is normally not converging. E.g. if we truncate the spectrum at $n=k$ one gets $$ \sum_{n=0}^k e^{i(2n+1)\epsilon_0 t}= \frac{e^{i(3+2k)\epsilon_0 t}-e^{i\epsilon_0 t}}{e^{i2\epsilon_0 t}-1} $$ that for the exponential becoming $1$ yields that the sum goes to infinity just like $k$. But physicists have a lot of resources to cope with these situations. We can resort to one of the techniques in the Hardy's book and introduce a converging factor into the series as $$ \sum_{n=0}^\infty e^{i(2n+1)\epsilon_0 t-\delta n} = \frac{e^{\delta+i\epsilon_0 t}}{-e^{2i\epsilon_0 t}+e^{\delta}} $$ and the limit $\delta\rightarrow 0$ yields the result we wanted. So, finally $$ D(\epsilon)=-\int_{-\infty}^{+\infty}\frac{dt}{4\pi i}e^{-i\epsilon t}\frac{1}{\sin\epsilon_0 t} $$ that is the final result provided we add a rule to circumvent all the poles arising due to the sine function at the denominator (see below). Just notice that for $\epsilon_0 t\ll 1$ one has $\sin\epsilon_0 t\approx \epsilon_0 t$. Now, you have to add a rule to circumvent the pole at $t=0$ and this is done in the standard way by adding a $i\epsilon$ at the denominator yielding $$ D(\epsilon)=-\frac{1}{\hbar\omega}\int_{-\infty}^{+\infty}\frac{dt}{2\pi i}e^{-i\epsilon t}\frac{1}{t+i\epsilon}. $$ This is exactly the definition of the $\theta$ function and so $$ D(\epsilon)\approx \frac{1}{\hbar\omega}\theta(\epsilon). $$
