How find this sum $\sum\limits_{n=0}^{\infty}\frac{1}{(3n+1)(3n+2)(3n+3)}$ Find this sum

$$I=\sum_{n=0}^{\infty}\dfrac{1}{(3n+1)(3n+2)(3n+3)}$$

My try: let

$$f(x)=\sum_{n=0}^{\infty}\dfrac{x^{3n+3}}{(3n+1)(3n+2)(3n+3)},|x|\le 1$$

then we have
$$f^{(3)}(x)=\sum_{n=0}^{\infty}x^{3n}=\dfrac{1}{1-x^3}$$
then we find the $f(x)$,Following is very ugly(can you someone can post your follow solution,) have other simple methods? Thank you very much.
 A: OK, so let's start with 
$$\frac{d^3}{dx^3} f(x) = \frac{1}{1-x^3}$$
where $f(0)=f'(0)=f''(0)=0$.  Note that we may use partial fractions and a little rearranging to find that
$$\frac{1}{1-x^3} = \frac13 \left [\frac{1}{1-x} + \frac{x+\frac12}{\left (x+\frac12\right)^2+\frac{3}{4}}+\frac{3}{2} \frac{1}{\left (x+\frac12\right)^2+\frac{3}{4}}\right ]$$
Integrating and using $f''(0)=0$ we get that
$$f''(x) = -\frac13 \log{(1-x)} + \frac16 \log{(1+x+x^2)} + \frac{1}{\sqrt{3}} \arctan{\frac{2 x+1}{\sqrt{3}}} - \frac{\pi}{6 \sqrt{3}}$$
From here the challenge will be managing the various terms rather than individual integrations, which are easily done through parts.  Note that, for example,
$$\int dx \, \log{(1+x+x^2)} = x \log{(1+x+x^2)} - \int dx \frac{2 x^2+x}{x^2+x+1}$$
You may use partial fractions on the latter integral and determine that
$$\int dx \frac{2 x^2+x}{x^2+x+1} = 2 x - \frac12 \log{(1+x+x^2)} - \sqrt{3} \arctan{\frac{2 x+1}{\sqrt{3}}} +C$$
Using a substitution ($u=(2 x+1)/\sqrt{3}$) and integration by parts, we similarly find that
$$\int dx \, \arctan{\frac{2 x+1}{\sqrt{3}}} = \frac{\sqrt{3}}{2}\left [\frac{2 x+1}{\sqrt{3}} \arctan{\frac{2 x+1}{\sqrt{3}}}  - \frac12 \log{(1+x+x^2)}\right ] + C$$
Omitting a lot of algebra, I find that, using $f'(0)=0$:
$$f'(x) = \frac13 (1-x) \log{(1-x)} - \frac16 (1-x) \log{(1+x+x^2)} + \frac{1+x}{\sqrt{3}}  \arctan{\frac{2 x+1}{\sqrt{3}}} - \frac{\pi}{6 \sqrt{3}} (1+x)$$
Now we integrate once more.  Again, the challenge is bookkeeping, as all the integrals may be done by parts.  Using $f(0)=0$, I get that
$$f(x) = -\frac16(1-x)^2 \log{(1-x)} + \frac{1}{12} (x^2-2 x-2) \log{(1+x+x^2)} + \frac{(x^2+2 x)}{2 \sqrt{3}} \arctan{\frac{2 x+1}{\sqrt{3}}}  - \frac{\pi}{12 \sqrt{3}} (x^2+2 x)$$
The sum is then 
$$f(1) = \frac{\sqrt{3}}{12} \pi - \frac14 \log{3} \approx 0.178$$
which checks out with, e.g., WA.
A: Your try is excellent... Note that $f(0)=f'(0)=f''(0)=0$ hence
$$
f(1)=\int_0^1f'(x)\,\mathrm dx=\int_0^1\int_0^xf''(y)\,\mathrm dy\,\mathrm dx=\int_0^1\int_0^x\int_0^yf'''(z)\,\mathrm dz\,\mathrm dy\,\mathrm dx,
$$
that is,
$$
2f(1)=\int_0^1(1-z)^2f'''(z)\,\mathrm dz=\int_0^1\frac{1-z}{1+z+z^2}\,\mathrm dz.
$$
The rest is routine. The change of variable $2z+1=\sqrt3t$ yields
$$
2f(1)=\int_{1/\sqrt3}^\sqrt3\frac{\sqrt3-t}{1+t^2}\,\mathrm dt=\left[\sqrt3\arctan t-\frac12\log(1+t^2)\right]_{1/\sqrt3}^\sqrt3.
$$
Note that $\arctan\sqrt3=\pi/3$ and $\arctan1/\sqrt3=\pi/6$, hence
$$
2f(1)=\sqrt3\cdot\left(\frac\pi3-\frac\pi6\right)-\frac12\log4+\frac12\log\frac43,
$$
that is,
$$
f(1)=\frac14\left[\frac\pi{\sqrt3}-\log3\right]\approx0.1788.
$$
Second method: The rational fraction is such that
$$
\frac2{(3n+1)(3n+2)(3n+3)}=\frac1{3n+1}-\frac2{3n+2}+\frac1{3n+3},
$$
hence
$$
2f(x)=x^2g_1(x)-2xg_2(x)+g_3(x),\qquad g_k(x)=\sum_{n\geqslant0}\frac{x^{3n+k}}{3n+k}.
$$
Thus, for each $k$,
$$
g'_k(x)=\sum_{n\geqslant0}x^{3n+k-1}=\frac{x^{k-1}}{1-x^3}.
$$
Since $g'_k(0)=0$ for every $k\geqslant1$, this yields
$$
2f(x)=x^2\int_0^x\frac{1}{1-t^3}\mathrm dt-2x\int_0^x\frac{t}{1-t^3}\mathrm dt+\int_0^x\frac{t^2}{1-t^3}\mathrm dt.
$$
The change of variable $t=xu$ yields
$$
2f(x)=x^3\int_0^1\frac{1}{1-x^3u^3}(1-2u+u^2)\mathrm du,
$$
that is,
$$
2f(x)=x^3\int_0^1\frac{1-u}{1-xu}\frac{1-u}{1+xu+x^2u^2}\mathrm du.
$$
When $x\to1$, one obtains once again
$$
2f(1)=\int_0^x\frac{1-u}{1+u+u^2}\mathrm du.
$$
More generally, for every integer $k\geqslant2$,
$$
\sum_{n\geqslant0}\frac{(k-1)!}{(kn+1)(kn+2)\cdots(kn+k)}=\int_0^1\frac{(1-u)^{k-2}}{1+u+\cdots+u^{k-1}}\,\mathrm du.
$$
A: I undeleted my answer and fixed it up since one of the other answers refers to it.

In the comments, Lucian mentioned the generalization $$\sum_{n=0}^{\infty} \frac{1}{(n+a)(n+b)(n+c)} = \frac{(b-c) \, \psi(a)+ (c-a) \, \psi(b)+(a-b) \, \psi(c)}{(a-b)(b-c)(c-a)} \tag{1}$$ where $\psi(x)$ is the digamma function;  $a,b,c >0$; $a \ne b$; $b \ne c$; and $a \ne c$.
We can prove $(1)$ by using the series representation $$\psi(x) = - \gamma + \sum_{n=0}^{\infty} \left(\frac{1}{n+1}-\frac{1}{n+x} \right) \, , \quad x>0.$$
Starting from the right side of $(1)$, we get
$$ \begin{align} &\frac{(b-c) \, \psi(a)+ (c-a) \, \psi(b)+(a-b) \, \psi(c)}{(a-b)(b-c)(c-a)} \\   &= \small\left[ (b-c)\left(-\gamma + \sum_{n=0}^{\infty} \left( \frac{1}{n+1} - \frac{1}{n+a} \right) \right)+ (c-a)\left(-  \gamma + \sum_{n=0}^{\infty} \left( \frac{1}{n+1} - \frac{1}{n+b} \right) \right)   + (a-b) \left( -\gamma +\sum_{n=0}^{\infty} \left( \frac{1}{n+1} - \frac{1}{n+c} \right) \right) \right] /{\left((a-b)(b-c)(c-a) \right)} \\  &= \frac{1}{(a-b)(b-c)(c-a)}\sum_{n=0}^{\infty} \left( \frac{c-b}{n+a}+\frac{a-c}{n+b} + \frac{b-a}{n+c} \right) \\ &=\frac{1}{(a-b)(b-c)(c-a)} \sum_{n=0}^{\infty} \ \frac{bc^{2}-b^{2}c+a^{2}c-ac^{2}+ab^{2}-a^{2}b}{(n+a)(n+b)(n+c)} \\ &= \frac{1}{(a-b)(b-c)(c-a)}\sum_{n=0}^{\infty} \frac{(a-b)(ac-c^{2}-ab+bc)}{(n+a)(n+b)(n+c)} \\ &=\frac{1}{(a-b)(b-c)(c-a)}\sum_{n=0}^{\infty} \frac{(a-b)(b-c)(c-a)}{(n+a)(n+b)(n+c)} \\ &= \sum_{n=0}^{\infty} \frac{1}{(n+a)(n+b)(n+c)} \end{align}$$
A: We all know that sums of the form $\displaystyle{\sum\frac1{(n+a)(n+b)}}$ or $\displaystyle{\sum\frac1{(n+a)(n+b)(n+c)}}$ etc., are telescopic in nature, and solving them is trivial: for integer values of a, b, and c, that is ! But what if a, b, and c are not integers? What then? Suddenly, things aren't so simple and trivial anymore, and the formerly banal and deceitfully tame-looking problem takes on whole new twist, gaining entirely unexpected dimensions of depth, meaning, and insight. The key, when dealing with such surprising turns of events, is to just take a step back, and try to parse or rephrase the older, well-worn solution in terms which might prove useful or relevant to circumventing the obstacles raised by the newly-encountered situation.

For instance, it is by no means difficult to show that $\displaystyle{\sum_{n=1}^\infty\frac1{(n+a)(n+b)}}=\frac{H_a-H_b}{a-b}$ , or that $$\sum_{n=1}^\infty\frac1{(n+a)(n+b)(n+c)}=\frac{(a-b)\cdot H_c\ +\ (b-c)\cdot H_a\ +\ (c-a)\cdot H_b}{(a-b)(b-c)(c-a)}$$ for natural values of a, b, c. Indeed, the user Random Variable has already proved it on this thread, though I  now realize that my initial formula is slightly mistaken, in the sense that either the sum should begin at $0$, with the $\gamma$ present, or at $1$, with the $\gamma$ absent. Now, the whole question becomes how to redefine $H_m$ , so as to be able to extend its meaning to non-natural arguments. For natural arguments, we have the well-known formula $H_m=\displaystyle\sum_{k=1}^m\frac1k$ . Now, let us take the simple function $f_k(x)=\displaystyle\frac{x^k}k$ , and notice that $f_k'(x)=x^{k-1}$ . And since we know that $\displaystyle{\sum_{k=1}^mx^{k-1}=\frac{1-x^m}{1-x}}$ , this finally allows us to conclude that $\displaystyle{H_m=\int_0^1\frac{1-x^m}{1-x}dx}$ , which, unlike the previous formula, can easily be extended to non-natural arguments as well. This has already been done by Euler two and a half centuries ago, so it's hardly new territory. Now, in our case, $\{a,b,c\}=\left\{\frac13,\frac23,1\right\}$ , and the value of $H_1$ is $1$, so all that's left to do is to compute $H_\frac13$ and $H_\frac23$ using the above formula, since $$\sum_{n=0}^\infty\frac1{(3n+1)(3n+2)(3n+3)}=\frac16+\frac1{27}\cdot\sum_{n=1}^\infty\frac1{\left(n+\frac13\right)\left(n+\frac23\right)(n+1)}$$

$$H_\frac13=\int_0^1\frac{1-\sqrt[3]x}{1-x}dx=\int_0^1\frac{1-t}{1-t^3}d\left(t^3\right)=\int_0^1\frac1{1+t+t^2}d\left(t^3\right)=\int_{\sqrt[3]0}^{\sqrt[3]1}\frac{3\,t^2}{1+t+t^2}dt=$$
$$=3\int_0^1\left(1-\frac{1+t}{1+t+t^2}\right)dt=3\,\bigg[\int_0^11\cdot dt-\tfrac12\int_0^1\frac{2+2t}{1+t+t^2}dt\bigg]=$$
$$=3\,\bigg[t|_0^1-\tfrac12\bigg(\int_0^1\frac1{1+t+t^2}dt+\int_0^1\frac{1+2t}{1+t+t^2}dt\bigg)\bigg]=$$
$$=3\,\bigg[1-\tfrac12\int_0^1\frac1{\left(t+\frac12\right)^2+\frac34}dt-\tfrac12\cdot\ln\left(1+t+t^2\right)_0^1\bigg]=$$
$$=3\,\bigg[1-\tfrac12\cdot\frac1{\sqrt3/2}\cdot\arctan\bigg(\frac{t+\frac12}{\sqrt3\big/2}\bigg)_0^1-\frac{\ln3}2\bigg]=3-\frac\pi{2\sqrt3}-\tfrac32\ln3.$$

Similarly for $\displaystyle{H_{\frac23}=\tfrac32+\frac\pi{2\sqrt3}-\tfrac32\ln3}$ . Then, by substituting these values back into the original formula above, we finally arrive at the desired result $\displaystyle{I=\frac{\pi\sqrt3-3\ln3}{12}}$ .
