Number of $n$-element subsets of $\{1, 2, \dotsc, 3n\}$ with sum divisible by $n$

Question

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} $$

Answer 1

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.