Computing $\sum\limits_{n=1}^q {1\over (kn-1)(kn)(kn+1)}$ as a telescoping series, when $k\geqslant3$?

One knows that $$S_1(q)=\sum_{n=1}^q {1\over (n-1)(n)(n+1)} = \sum_{n=1}^q\frac12\left({\frac{1}{(n-1)n}-\frac{1}{n(n+1)}}\right)$$ and the RHS can be easily telescoped. The same approach works for $$S_2(q)=\sum_{n=1}^q {1\over (2n-1)(2n)(2n+1)}$$ However, for $$S_3(q)=\sum_{n=1}^q {1\over (3n-1)(3n)(3n+1)}$$ it is impossible to telescope using the same method than in the two cases above. So:

How should $$S_k(q)=\sum_{n=1}^q {1\over (kn-1)(kn)(kn+1)}$$ where $k\geqslant3$ is an integer, be telescoped?

As you have noted, in the case of $k=1$, we can use Partial Fractions: $$ \begin{align} \frac1{(n-1)n(n+1)} &=\frac{1/2}{n-1}-\frac1n+\frac{1/2}{n+1}\\ &=\frac12\left(\frac1{n-1}-\frac1n\right)-\frac12\left(\frac1n-\frac1{n+1}\right)\tag{1} \end{align} $$

The case for $k=2$ is not as easy $$ \frac1{(2n-1)2n(2n+1)} =\frac{1/2}{2n-1}-\frac1{2n}+\frac{1/2}{2n+1}\tag{2} $$ The partial sums are hard to evaluate, though the infinite sum can be computed as an alternating harmonic series: $$ \begin{align} \sum_{n=1}^\infty\frac1{(2n-1)2n(2n+1)} &=\sum_{n=1}^\infty\left(\frac{1/2}{2n-1}-\frac1{2n}+\frac{1/2}{2n+1}\right)\\ &=\sum_{n=1}^\infty\left(\frac{1/2}{2n-1}-\frac1{2n}\right)+\sum_{n=2}^\infty\left(\frac{1/2}{2n-1}\right)\\ &=\underbrace{\sum_{n=1}^\infty\left(\frac1{2n-1}-\frac1{2n}\right)}_{\text{alternating harmonic series}}-\frac12\\ &=\log(2)-\frac12\tag{3} \end{align} $$

For $k\gt2$, even the infinite sums are hard because the terms for different $k$ do not overlap. Something other than telescoping series needs to be used. $$ \begin{align} \sum_{n=1}^\infty\frac1{(kn-1)kn(kn+1)} &=\sum_{n=1}^\infty\left(\frac{1/2}{kn-1}-\frac1{kn}+\frac{1/2}{kn+1}\right)\\ &=-\frac1{2k}\sum_{n=1}^\infty\left(\frac1n-\frac1{n-1/k}\right)-\frac1{2k}\sum_{n=1}^\infty\left(\frac1n-\frac1{n+1/k}\right)\\ &=-\frac1{2k}\left(H_{-1/k}+H_{1/k}\right)\tag{4} \end{align} $$ where $H_x$ is the extended Harmonic number as described in this answer. That answer shows how to compute certain values, but more can be computed.

Extended Harmonic Numbers with Rational Arguments

In extension of the values derived in this answer, we can compute the value of $H_{p/q}$ for any rational $p/q$.

We will use that for $x\in[-2\pi,2\pi]$ $$ \log\left(1-e^{ix}\right) =\log\left(2\sin(x/2)\right)+ix/2-\pi i/2\operatorname{sgn}(\sin(x/2))\tag{5} $$ and $$ \frac1q\sum_{k=0}^{q-1}e^{2\pi ijk/q}=[j\equiv0\pmod{q}]\tag{6} $$ where $[\cdots]$ are Iverson brackets.

For $0\lt p\le q$, $$ \begin{align} H_{p/q} &=\sum_{n=1}^\infty\left(\frac1n-\frac1{n+p/q}\right)\\ &=q\sum_{n=1}^\infty\left(\frac1{qn}-\frac1{qn+p}\right)\\ &=\frac qp+q\sum_{j=1}^\infty\left(\vphantom{\frac1j}[j\equiv0\pmod{q}]-[j\equiv p\pmod{q}]\right)\frac1j\\ &=\frac qp+q\sum_{j=1}^\infty\frac1q\sum_{k=1}^{q-1}\left(e^{2\pi ijk/q}-e^{2\pi i(j-p)k/q}\right)\frac1j\\ &=\frac qp+q\sum_{j=1}^\infty\frac1q\sum_{k=1}^{q-1}e^{2\pi ijk/q}\left(1-e^{-2\pi ipk/q}\right)\frac1j\\ &=\frac qp-\sum_{k=1}^{q-1}\log\left(1-e^{2\pi ik/q}\right)\left(1-e^{-2\pi ikp/q}\right)\\ &=\frac qp-\sum_{k=1}^{q-1}\left[\scriptsize\log\left(2\sin\left(\frac{\pi k}q\right)\right)+i\left(\frac{\pi k}q-\frac\pi2\right)\right]\left[\scriptsize\left(1-\cos\left(\frac{2\pi kp}q\right)\right)+i\sin\left(\frac{2\pi kp}q\right)\right]\\ &=\bbox[5px,border:2px solid #C0A000]{\frac qp-\sum_{k=1}^{q-1}{\scriptsize\left[\log\left(2\sin\left(\frac{\pi k}q\right)\right)\left(1-\cos\left(\frac{2\pi kp}q\right)\right)+\left(\frac{\pi}2-\frac{\pi k}q\right)\sin\left(\frac{2\pi kp}q\right)\right]}}\tag{7} \end{align} $$ Because $(7)$ is valid for $0\lt p\le q$, we can use $$ H_{-1/k}=H_{(k-1)/k}-\frac{k}{k-1}\tag{8} $$ Thus, we get $$ \bbox[5px,border:2px solid #C0A000]{\sum_{n=1}^\infty\frac1{(kn-1)kn(kn+1)} =-\frac1{2k}\left(H_{(k-1)/k}+H_{1/k}-\frac{k}{k-1}\right)}\tag{9} $$ where each of the Harmonic numbers can be computed with $(7)$.

Here is a list of sums for the first $6$ values of $k$ $$ \begin{array}{c|c} k&\sum\limits_{n=1}^\infty\frac1{(kn-1)kn(kn+1)}\\\hline 1&\infty\text{, or $\frac14$ if we sum from $n=2$}\\ 2&-\frac12+\log(2)\\ 3&-\frac12+\frac12\log(3)\\ 4&-\frac12+\frac34\log(2)\\ 5&\scriptsize-\frac12-\frac12\log(2)+\frac{5+\sqrt5}{20}\log\left(5+\sqrt5\right)+\frac{5-\sqrt5}{20}\log\left(5-\sqrt5\right)\\ 6&-\frac12+\frac13\log(2)+\frac14\log(3) \end{array}\tag{10} $$ Mathematica implementation of $(7)$:

h[p_,q_] := q/p-Sum[(1-Cos[2Pi k p/q])Log[2Sin[Pi k/q]] + Sin[2Pi k p/q](1/2-k/q)Pi,{k,1,q-1}]

Partial fractions:

$$\frac1{(3n-1)3n(3n+1)}=\frac a{3n-1}+\frac b{3n}+\frac c{3n+1}\implies$$

$$1=3an(3n+1)+b(3n-1)(3n+1)+3cn(3n-1)$$

Now, attach values to $\;n\;$ to get the values $\;a,b,c\;$ :

$$\begin{align*}&n=0\implies&1=-b\\{}\\ &n=\frac13\implies&1=2a\implies a=\frac12\\{}\\ &n=-\frac13\implies&-1=-2c\implies c=\frac12\end{align*}$$

Thus:

$$\frac1{(3n-1)3n(3n+1)}=\frac12\left(\frac1{3n-1}-\frac2{3n}+\frac1{3n+1}\right)$$