Number of $n$-element subsets of $\{1, 2, \dotsc, 3n\}$ with sum divisible by $n$ 
Is there any recurrence relations/formula/algorithm to count the number of $n$-element subsets of the set $\{1, 2, \dotsc, 3n\}$ with sum divisible by $n$? How about replacing $3n$ with $kn$?

I currently only know that if $k=2$, then the number of $n$-element subsets divisible by $n$ is
$$     
  \frac{(-1)^n}{n}\sum_{d \mid n} (-1)^d\phi({n\over d})\binom{2d}{d}
$$
from http://oeis.org/A169888, but a similar search for $k=3$ doesn't return any results. 
Is there any similar formula for $k \geq 3$?
Thanks.
Some bruteforced values for $k=3$:
$$
  \begin{array}{c|ccccccccc}
      n
    & 1
    & 2
    & 3
    & 4
    & 5
    & 6
    & 7
    & 8
    & 9
    \\
    \hline
      \text{$n$-element subsets}
    & 3
    & 6
    & 30
    & 126
    & 603
    & 3084
    & 16614
    & 91998
    & 520779
    \\
  \end{array}
$$
 A: We ask about the probability that a  set of size $n$ drawn from $[kn]$
has sum  divisible by $n$. The  exponential formula tells us  that the
cycle index $Z(P_n)$  of the unlabeled set  operator
$$\def\textsc#1{\dosc#1\csod}
\def\dosc#1#2\csod{{\rm #1{\small #2}}}\textsc{SET}$$ 
on $n$ slots has OGF
$$Z(P_n) = [w^n]
\exp\left(\sum_{l\ge 1} (-1)^{l-1} a_l \frac{w^l}{l}\right).$$
The desired probability is given by
$${kn\choose n}^{-1} \frac{1}{n}
\sum_{p=0}^{n-1} 
\left. Z
\left(P_n; \sum_{q=1}^{kn} z^q\right)\right|
_{z=\exp(2\pi ip/n)}.$$
This is
$${kn\choose n}^{-1} \frac{1}{n} 
\sum_{p=0}^{n-1} \left. [w^n]
\exp
\left(\sum_{l\ge 1} (-1)^{l-1} \left(\sum_{q=1}^{kn} z^{ql}\right)
\frac{w^l}{l}\right)\right|_{z=\exp(2\pi ip/n)}.$$
Evaluating the contribution for $p=0$ first we get
$${kn\choose n}^{-1} \frac{1}{n} 
[w^n] \exp\left(\sum_{l\ge 1} (-1)^{l-1} kn
\frac{w^l}{l}\right)
\\ = {kn\choose n}^{-1} \frac{1}{n} 
[w^n] \exp\left(kn \log(1+w)\right)
\\ = {kn\choose n}^{-1} \frac{1}{n} 
[w^n] (1+w)^{kn}
= {kn\choose n}^{-1} \frac{1}{n} {kn\choose n}
= \frac{1}{n}.$$
It remains to evaluate the contribution  from $1\le p\le n-1.$ Now for
these $p$ if $l$ is a multiple of $m = n/\gcd(p, n)$ we have
$$\sum_{q=1}^{kn} \exp(2\pi ip/n)^{ql} = kn.$$
We get zero otherwise. This yields for the remaining terms without the
scalar in front
$$\sum_{p=1}^{n-1} [w^n]
\exp\left(\sum_{l\ge 1} (-1)^{ml-1} kn
\frac{w^{ml}}{ml} \right)
= \sum_{p=1}^{n-1} [w^n]
\exp\left(-\frac{kn}{m} \sum_{l\ge 1} 
\frac{(-w)^{ml}}{l} \right)
\\ = \sum_{p=1}^{n-1} [w^n]
\exp\left(-\frac{kn}{m}\log\frac{1}{1-(-w)^m}\right)
= \sum_{p=1}^{n-1} [w^n] (1-(-w)^{n/\gcd(p, n)})^{k\gcd(p, n)}
\\ = \sum_{p=1}^{n-1} [w^n] (1+(-1)^{1+n/\gcd(p, n)}
w^{n/\gcd(p, n)})^{k\gcd(p, n)}.$$
This is
$$\sum_{p=1}^{n-1} {k\gcd(p,n) \choose \gcd(p,n)}
(-1)^{(1+n/\gcd(p, n)) \gcd(p,n)}.$$
Putting it all together we thus obtain
$$\bbox[5px,border:2px solid #00A000]{
\frac{1}{n} + (-1)^n
{kn\choose n}^{-1} \frac{1}{n} 
\sum_{p=1}^{n-1} {k\gcd(p,n) \choose \gcd(p,n)}
(-1)^{\gcd(p,n)}.}$$
While this formula will produce results it may perhaps be simplified.
Write
$$\frac{1}{n} = 
(-1)^n {kn\choose n}^{-1} \frac{1}{n} 
{k\gcd(n,n)\choose \gcd(n,n)} (-1)^{\gcd(n,n)}$$
to get
$$(-1)^n
{kn\choose n}^{-1} \frac{1}{n} 
\sum_{p=1}^{n} 
{k\gcd(p, n)\choose \gcd(p,n)} 
(-1)^{\gcd(p,n)}
\\ = (-1)^n
{kn\choose n}^{-1} \frac{1}{n} 
\sum_{d|n} \sum_{\gcd(p,n)=d}
{kd\choose d} 
(-1)^{d}
\\ = (-1)^n
{kn\choose n}^{-1} \frac{1}{n} 
\sum_{d|n} {kd\choose d} 
(-1)^{d} \sum_{\gcd(q,n/d)=1} 1.$$
We find the alternate closed form
$$\bbox[5px,border:2px solid #00A000]{
(-1)^n
{kn\choose n}^{-1} \frac{1}{n} 
\sum_{d|n} {kd\choose d} 
(-1)^d \varphi(n/d).}$$
These two formulae were verified  by simple enumeration for $1\le n\le
7$ and $1\le k\le 6,$ see below.

with(combinat);
with(numtheory);

ENUM :=
proc(k, n)
option remember;
local recurse, admit, total;

    admit := 0; total := 0;

    recurse :=
    proc(pos, selcount, sumsofar)
        if selcount = n then
            if sumsofar mod n = 0 then
                admit := admit + 1;
            fi;

            total := total + 1;

            return;
        fi;

        if pos > k*n then return fi;

        recurse(pos+1, selcount, sumsofar);
        recurse(pos+1, selcount + 1,
                sumsofar + pos);
    end;

    recurse(1, 0, 0);
    admit/total;
end;

X := (k, n) ->
1/n+(-1)^n*binomial(k*n,n)^(-1)*1/n
* add(binomial(k*gcd(p,n),gcd(p,n))*(-1)^gcd(p,n),
      p=1..n-1);

XX := (k, n) ->
(-1)^n*binomial(k*n,n)^(-1)*1/n
* add(binomial(k*d,d)*(-1)^d*phi(n/d),
      d in divisors(n));

Remark. This  computation is based  on material from  this MSE
link.
