Krdytkyu

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)^dphi({nover 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}
$$

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

$$deftextsc#1{dosc#1csod}
defdosc#1#2csod{{rm #1{small #2}}}textsc{SET}$$

on $n$ slots has OGF

$$Z(P_n) = [w^n]
expleft(sum_{lge 1} (-1)^{l-1} a_l frac{w^l}{l}right).$$

The desired probability is given by

$${knchoose n}^{-1} frac{1}{n}
sum_{p=0}^{n-1}
left. Z
left(P_n; sum_{q=1}^{kn} z^qright)right|
_{z=exp(2pi ip/n)}.$$

This is

$${knchoose n}^{-1} frac{1}{n}
sum_{p=0}^{n-1} left. [w^n]
exp
left(sum_{lge 1} (-1)^{l-1} left(sum_{q=1}^{kn} z^{ql}right)
frac{w^l}{l}right)right|_{z=exp(2pi ip/n)}.$$

Evaluating the contribution for $p=0$ first we get

$${knchoose n}^{-1} frac{1}{n}
[w^n] expleft(sum_{lge 1} (-1)^{l-1} kn
frac{w^l}{l}right)
\ = {knchoose n}^{-1} frac{1}{n}
[w^n] expleft(kn log(1+w)right)
\ = {knchoose n}^{-1} frac{1}{n}
[w^n] (1+w)^{kn}
= {knchoose n}^{-1} frac{1}{n} {knchoose n}
= frac{1}{n}.$$

It remains to evaluate the contribution from $1le ple n-1.$ Now for
these $p$ if $l$ is a multiple of $m = n/gcd(p, n)$ we have

$$sum_{q=1}^{kn} exp(2pi 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]
expleft(sum_{lge 1} (-1)^{ml-1} kn
frac{w^{ml}}{ml} right)
= sum_{p=1}^{n-1} [w^n]
expleft(-frac{kn}{m} sum_{lge 1}
frac{(-w)^{ml}}{l} right)
\ = sum_{p=1}^{n-1} [w^n]
expleft(-frac{kn}{m}logfrac{1}{1-(-w)^m}right)
= sum_{p=1}^{n-1} [w^n] (1-(-w)^{n/gcd(p, n)})^{kgcd(p, n)}
\ = sum_{p=1}^{n-1} [w^n] (1+(-1)^{1+n/gcd(p, n)}
w^{n/gcd(p, n)})^{kgcd(p, n)}.$$

This is

$$sum_{p=1}^{n-1} {kgcd(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
{knchoose n}^{-1} frac{1}{n}
sum_{p=1}^{n-1} {kgcd(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 {knchoose n}^{-1} frac{1}{n}
{kgcd(n,n)choose gcd(n,n)} (-1)^{gcd(n,n)}$$

to get

$$(-1)^n
{knchoose n}^{-1} frac{1}{n}
sum_{p=1}^{n}
{kgcd(p, n)choose gcd(p,n)}
(-1)^{gcd(p,n)}
\ = (-1)^n
{knchoose n}^{-1} frac{1}{n}
sum_{d|n} sum_{gcd(p,n)=d}
{kdchoose d}
(-1)^{d}
\ = (-1)^n
{knchoose n}^{-1} frac{1}{n}
sum_{d|n} {kdchoose d}
(-1)^{d} sum_{gcd(q,n/d)=1} 1.$$

We find the alternate closed form

$$bbox[5px,border:2px solid #00A000]{
(-1)^n
{knchoose n}^{-1} frac{1}{n}
sum_{d|n} {kdchoose d}
(-1)^d varphi(n/d).}$$

These two formulae were verified by simple enumeration for $1le nle
7$ and $1le kle 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.